Invalid WCS when saving moment0 maps of non-spectral axes

[ChrisBeaumont]

Works:

cube.moment0(axis=0).write('test.fits')
cube.moment0(axis=1)

Fails:

In [13]: s.moment0(axis=1).write('test.fits')
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-13-2419fce08bed> in <module>()
----> 1 s.moment0(axis=1).write('test.fits')

/Users/beaumont/alma/spectral-cube/spectral_cube/spectral_cube.pyc in write(self, filename, format, overwrite)
    124             format = determine_format(filename)
    125         if format == 'fits':
--> 126             self.hdu.writeto(filename, clobber=overwrite)
    127         else:
    128             raise ValueError("Unknown format '{0}' - the only available "

/Users/beaumont/alma/spectral-cube/spectral_cube/spectral_cube.pyc in hdu(self)
    104             hdu = fits.PrimaryHDU(self.value)
    105         else:
--> 106             hdu = fits.PrimaryHDU(self.value, header=self.wcs.to_header())
    107         hdu.header['BUNIT'] = self.unit.to_string(format='fits')
    108         return hdu

/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.

[ChrisBeaumont]

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.

[keflavich]

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?

[keflavich]

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)

[ChrisBeaumont]

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

[mdboom]

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.

[mateoi]

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

[keflavich]

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: #367