Fix indexing for propagating PV and PS cards.

[low-sky]

PVi_j cards are being lost when reducing dimensions. This appears to come from wcs.get_pv() returning tuples with the coordinates given in FITS 1-based axis indexing whereas the list of what axes to include (inds) is 0-based. Included +1 to shift axes so they line up.

[keflavich]

The blame on this points to me, but damned if I know where that logic comes from, and the code style doesn't look like mine, so there's a good chance this code snippet exists somewhere else. I want to hunt it down in case this also exists as a failure mode in astropy somewhere.

[keflavich]

Ah, here's the issue I was referencing:
astropy/astropy#3588 (comment)

[keflavich]

@low-sky you were encountering this issue when taking moments, right? We definitely want to add a regression test.

[keflavich]

This is a minimalist example I tried based on your earlier report:

from spectral_cube import SpectralCube
from astropy import wcs
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord
import astropy.units as u
s=cube = SpectralCube.read(('G008.67_B6_spw1_12M_sio.image.fits'))
pos = SkyCoord(271.58799991,
               -21.62130537,
               unit=(u.deg, u.deg), frame='fk5')
x0, y0, _ = s.wcs.wcs_world2pix(pos.ra.value, pos.dec.value, 0, 0)
box = 4
subcube = s[:3, int(y0-box):int(y0+box), int(x0-box):int(x0+box)]
moment = subcube.moment0()
x1, y1, _ = subcube.wcs.wcs_world2pix(pos.ra.value,
                                      pos.dec.value, 0, 0)
x2, y2 = subcube.wcs.celestial.wcs_world2pix(pos.ra.value,
                                             pos.dec.value, 0)
x3, y3 = moment.wcs.wcs_world2pix(pos.ra.value,
                                  pos.dec.value, 0)
print('Cube          :', x1, y1)
print('Cube.celestial:', x2, y2)
print('Moment        :', x3, y3)

the result is

Cube          : 4.999916332536277 4.000020250002777
Cube.celestial: 4.999916332536277 4.000020250002777
Moment        : 4.999916332536277 4.000020250002777

So this probably has to do with the type of projection - I'm guessing the Cygnus data you were working with is in Galactic?

[keflavich]

OK, yes, with a Galactic coordinate, I get... a closely related but different error:

IndexError                                Traceback (most recent call last)
<ipython-input-43-b1ae781b9bd3> in <module>
     11 box = 4
     12 subcube = s[:3, int(y0-box):int(y0+box), int(x0-box):int(x0+box)]
---> 13 moment = subcube.moment0()
     14 x1, y1, _ = subcube.wcs.wcs_world2pix(pos.ra.value,
     15                                       pos.dec.value, 0, 0)

/orange/adamginsburg/repos/spectral-cube/spectral_cube/spectral_cube.py in moment0(self, axis, how)
   1608            map [, wcs]
   1609            The moment map (numpy array) and, if wcs=True, the WCS object
-> 1610            describing the map
   1611 
   1612         Notes

/orange/adamginsburg/repos/spectral-cube/spectral_cube/spectral_cube.py in moment(self, order, axis, how)
   1571         channel and :math:`x` is the spectral coordinate:
   1572 
-> 1573         Moment 0:
   1574 
   1575         .. math:: M_0 \\int I dx

/orange/adamginsburg/repos/spectral-cube/spectral_cube/_moments.py in moment_auto(cube, order, axis)
    186     strategy = dict(cube=moment_cubewise, ray=moment_raywise,
    187                     slice=moment_slicewise)
--> 188     return strategy[iterator_strategy(cube, axis)](cube, order, axis)

/orange/adamginsburg/repos/spectral-cube/spectral_cube/_moments.py in moment_cubewise(cube, order, axis)
    159     """
    160 
--> 161     pix_cen = cube._pix_cen()[axis]
    162     data = cube._get_filled_data() * cube._pix_size_slice(axis)
    163 

/orange/adamginsburg/repos/spectral-cube/spectral_cube/utils.py in wrapper(self, *args)
     16         # The cache lives in the instance so that it gets garbage collected
     17         if (func, args) not in self._cache:
---> 18             self._cache[(func, args)] = func(self, *args)
     19         return self._cache[(func, args)]
     20 

/orange/adamginsburg/repos/spectral-cube/spectral_cube/spectral_cube.py in _pix_cen(self)
   1354             A rest wavelength or frequency with appropriate units.  Required if
   1355             output type is velocity.  The cube's WCS should include this
-> 1356             already if the *input* type is velocity, but the WCS's rest
   1357             wavelength/frequency can be overridden with this parameter.
   1358 

/orange/adamginsburg/repos/spectral-cube/spectral_cube/cube_utils.py in __getitem__(self, view)
    225         self._other = _other
    226 
--> 227     def __getitem__(self, view):
    228         result = self._func(self._other, view)
    229         if isinstance(result, da.Array):

/orange/adamginsburg/repos/spectral-cube/spectral_cube/base_class.py in world(self, view)
    187         inds = np.ogrid[[slice(0, s) for s in self.shape]]
    188         inds = np.broadcast_arrays(*inds)
--> 189         inds = [i[view] for i in inds[::-1]]  # numpy -> wcs order
    190 
    191         shp = inds[0].shape

/orange/adamginsburg/repos/spectral-cube/spectral_cube/base_class.py in <listcomp>(.0)
    187         inds = np.ogrid[[slice(0, s) for s in self.shape]]
    188         inds = np.broadcast_arrays(*inds)
--> 189         inds = [i[view] for i in inds[::-1]]  # numpy -> wcs order
    190 
    191         shp = inds[0].shape

IndexError: index 0 is out of bounds for axis 1 with size 0

which at least isn't silent!

[keflavich]

cube.wcs.celestial.wcs.get_pv() == cube.wcs.wcs.get_pv() is True for this cube.

i.e.:

>>> cube.wcs.celestial.wcs.get_pv()
[(1, 0, 1.0), (1, 1, 0.0), (1, 2, 0.0)]

[low-sky]

looking at the original astropy comment, should this just be replaced with a wcs.sub call? I’ll work try to work up a more minimal example.

[keflavich]

I think you're right that we should just use wcs.sub, as long as we can get the ordering to work. LMK if you have a minimal example that we can use as a regression test.

[keflavich]

We're still holding this back until we get a regression test