Open ChrisBeaumont opened 10 years ago
This is a more general error with the WCS, and not just a saving issue:
In [21]: s.moment0(axis=1).wcs.to_header_string()
ERROR: InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 1973 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
[astropy.wcs.wcs]
---------------------------------------------------------------------------
InconsistentAxisTypesError Traceback (most recent call last)
<ipython-input-21-62aead4e4bfb> in <module>()
----> 1 s.moment0(axis=1).wcs.to_header_string()
/Users/beaumont/anaconda/lib/python2.7/site-packages/astropy/wcs/wcs.pyc in to_header_string(self, relax)
1679 header cards.
1680 """
-> 1681 return str(self.to_header(relax))
1682
1683 def footprint_to_file(self, filename=None, color='green', width=2):
/Users/beaumont/anaconda/lib/python2.7/site-packages/astropy/wcs/wcs.pyc in to_header(self, relax, key)
1663
1664 if self.wcs is not None:
-> 1665 header_string = self.wcs.to_header(relax)
1666 header = fits.Header.fromstring(header_string)
1667 else:
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 1973 of file cextern/wcslib/C/wcs.c:
Unmatched celestial axes.
Right... WCS doesn't like having one of its celestial axes dropped. That should probably be allowed, but I guess it may not be at a deep level?
I think we need @mdboom's input - if we want an RA - frequency image plane (and it is genuinely aligned with RA), should that be allowed by WCS?
@ChrisBeaumont - What are we doing for non-aligned axes? I think in pvextractor
we just declared the spatial axis "arbitrary/offset" to work around this issue (no longer a "celestial axis" but some other unspecified spatial axis)
For non-aligned axes, I suspect the WCS is wrong (that is, I think it just drops the integrated-over axis). A generic offset axis (or two) is probably a better idea in these cases
That's an interesting question. I think wcslib doesn't like having a single celestial axis, because all of the projections are written assuming two, of course. You'd need to (somehow) define a constant for the missing axis, and I don't believe wcslib has a way of doing that. We could get Mark Calabretta involved in the discussion -- if you can boil this down to the minimum raw astropy.wcs calls that produces the situation, I can further convert that to C for his consideration of the problem.
Here's a method I wrote for SunPy that adds a useless axis. It's (mostly) untested, and it's a hack that doesn't really address the underlying problem. But still:
from astropy.wcs import WCS
import numpy as np
from astropy.wcs._wcs import InconsistentAxisTypesError
import re
def add_celestial_axis(wcs):
'''
Creates a copy of the given wcs and returns it, with an extra meaningless
celestial axes to allow for certain operations. The given WCS must already
have an unmatched celestial axis.
Parameters
----------
wcs: astropy.wcs.WCS object
The world coordinate system to add an axis to.
'''
outwcs = WCS(naxis=wcs.naxis + 1)
wcs_params_to_preserve = ['cel_offset', 'dateavg', 'dateobs', 'equinox',
'latpole', 'lonpole', 'mjdavg', 'mjdobs', 'name',
'obsgeo', 'phi0', 'radesys', 'restfrq',
'restwav', 'specsys', 'ssysobs', 'ssyssrc',
'theta0', 'velangl', 'velosys', 'zsource']
for par in wcs_params_to_preserve:
setattr(outwcs.wcs, par, getattr(wcs.wcs, par))
new_wcs_axes_params = {'crpix': [0], 'cdelt': [1], 'crval': [0],
'cname': ['redundant axis'], 'ctype': ['HPLN-TAN'],
'crota': [0], 'cunit': ['deg']}
try:
naxis = wcs.naxis
oldpc = wcs.wcs.pc
newpc = np.eye(naxis + 1)
newpc[:naxis, :naxis] = oldpc
outwcs.wcs.pc = newpc
except AttributeError:
pass
for param in new_wcs_axes_params:
try:
oldattr = list(getattr(wcs.wcs, param))
newattr = oldattr + new_wcs_axes_params[param]
setattr(outwcs.wcs, param, newattr)
except AttributeError: # Some attributes may not be present. Ignore.
pass
# Change the projection if we have two redundant celestial axes.
try:
outwcs.get_axis_types()
except InconsistentAxisTypesError as err:
projection = re.findall(r'expected [^,]+', err.message)[0][9:]
outwcs.wcs.ctype[-1] = projection
return outwcs
The more we dig into slicing, projecting, etc., the more it appears we need to preserve the full WCS when slicing. A 1D slice should have a length-1 coordinate for each of the sliced-out coordinates, and so on. This is fine and pretty easy for slices, but nontrivial for projections (what coordinate do you use when taking the max
along an axis?).
One of the main reasons I'd avoided this solution previously was the inconvenience and unexpected behavior of getting an array with shape (100,100,1,1)
when loading a 2D image that has a 4D WCS. As a user who (previously) didn't care at all about the WCS for most objects, I found that irritating and frustrating. To avoid passing this difficulty on to other users, though, we have already written the software (this package) that can appropriately determine the dimensions of FITS files and return the appropriate objects.
The solution I'm proposing will solve the related issue: https://github.com/radio-astro-tools/spectral-cube/issues/367
Works:
Fails: