nan in SpectralCube data range causes IndexError in subcube

[tvwenger]

Summary

I am attempting to read the HI4PI data cube with spectral_cube. When I do so, the reported data ranges are nan. Subsequently, any attempt to extract a sub-cube fails. I presume these issues are related.

Example

import spectral_cube
print(spectral_cube.__version__)
# 0.6.5
from spectral_cube import SpectralCube
cube = SpectralCube.read("CAR.fits")
print(cube)

Outputs:

SpectralCube with shape=(945, 2163, 4323) and unit=K:
 n_x:   4323  type_x: GLON-CAR  unit_x: deg    range:          nan deg:         nan deg
 n_y:   2163  type_y: GLAT-CAR  unit_y: deg    range:          nan deg:         nan deg
 n_s:    945  type_s: VRAD      unit_s: m / s  range:  -606648.293 m / s:  609426.638 m / s

Attempting to extract a sub-cube:

import astropy.units as u
cube.subcube(ylo=-1.0*u.deg, yhi=1.0*u.deg)

Raises:

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/spectral_cube.py:329, in BaseSpectralCube.__repr__(self)
    323 else:
    324     s += " and unit={0}:\n".format(self.unit)
    325 s += (" n_x: {0:6d}  type_x: {1:8s}  unit_x: {2:5s}"
    326       "  range: {3:12.6f}:{4:12.6f}\n".format(self.shape[2],
    327                                               self.wcs.wcs.ctype[0],
    328                                               self.wcs.wcs.cunit[0],
--> 329                                               self.longitude_extrema[0],
    330                                               self.longitude_extrema[1],))
    331 s += (" n_y: {0:6d}  type_y: {1:8s}  unit_y: {2:5s}"
    332       "  range: {3:12.6f}:{4:12.6f}\n".format(self.shape[1],
    333                                               self.wcs.wcs.ctype[1],
   (...)
    336                                               self.latitude_extrema[1],
    337                                              ))
    338 s += (" n_s: {0:6d}  type_s: {1:8s}  unit_s: {2:5s}"
    339       "  range: {3:12.3f}:{4:12.3f}".format(self.shape[0],
    340                                             self.wcs.wcs.ctype[2],
   (...)
    343                                             self.spectral_extrema[1],
    344                                            ))

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/utils.py:18, in cached.<locals>.wrapper(self, *args)
     14 @wraps(func)
     15 def 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)]

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/base_class.py:298, in SpatialCoordMixinClass.longitude_extrema(self)
    295 @property
    296 @cached
    297 def longitude_extrema(self):
--> 298     return self.world_extrema[0]

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/utils.py:18, in cached.<locals>.wrapper(self, *args)
     14 @wraps(func)
     15 def 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)]

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/base_class.py:281, in SpatialCoordMixinClass.world_extrema(self)
    279     latlon_corners = [self.world[y, x] for y,x in corners]
    280 else:
--> 281     latlon_corners = [self.world[0, y, x][1:] for y,x in corners]
    283 lon = u.Quantity([x for y,x in latlon_corners])
    284 lat = u.Quantity([y for y,x in latlon_corners])

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/base_class.py:281, in <listcomp>(.0)
    279     latlon_corners = [self.world[y, x] for y,x in corners]
    280 else:
--> 281     latlon_corners = [self.world[0, y, x][1:] for y,x in corners]
    283 lon = u.Quantity([x for y,x in latlon_corners])
    284 lat = u.Quantity([y for y,x in latlon_corners])

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/cube_utils.py:250, in SliceIndexer.__getitem__(self, view)
    249 def __getitem__(self, view):
--> 250     result = self._func(self._other, view)
    251     if isinstance(result, da.Array):
    252         result = result.compute()

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/base_class.py:219, in SpatialCoordMixinClass.world(self, view)
    217 inds = np.ogrid[[slice(0, s) for s in self.shape]]
    218 inds = np.broadcast_arrays(*inds)
--> 219 inds = [i[view] for i in inds[::-1]]  # numpy -> wcs order
    221 shp = inds[0].shape
    222 inds = np.column_stack([i.ravel() for i in inds])

File ~/miniconda3/envs/astro/lib/python3.11/site-packages/spectral_cube/base_class.py:219, in <listcomp>(.0)
    217 inds = np.ogrid[[slice(0, s) for s in self.shape]]
    218 inds = np.broadcast_arrays(*inds)
--> 219 inds = [i[view] for i in inds[::-1]]  # numpy -> wcs order
    221 shp = inds[0].shape
    222 inds = np.column_stack([i.ravel() for i in inds])

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

Expected Behavior

I did not expect to see a nan in the above data ranges for this data cube, and I expected to be able to successfully extract a sub-cube.

Context

Inspecting the header with astropy:

from astropy.io import fits
with fits.open("CAR.fits") as hdulist:
    print(hdulist[0].header)

Outputs:

SIMPLE  =                    T / conforms to FITS standard
BITPIX  =                  -32 / array data type
NAXIS   =                    3 / number of array dimensions
NAXIS1  =                 4323
NAXIS2  =                 2163
NAXIS3  =                  945
OBJECT  = 'HI4PI   '           / The HI 4-PI Survey
TELESCOP= 'Effelsberg 100m RT; ATNF Parkes 64-m' / Telescope names
ORIGIN  = 'AIfA/MPIfR Bonn; ATNF Sydney' / Organisations or Institutions
REFERENC= 'HI4PI Collaboration 2016' / A&A
RESTFRQ =        1420405751.77
RESTWAV =       0.211061140541
CDELT1  = -0.08333333330000001
CRPIX1  =               2162.0
CRVAL1  =                  0.0
CTYPE1  = 'GLON-CAR'
CUNIT1  = 'deg     '
CDELT2  =  0.08333333330000001
CRPIX2  =               1082.0
CRVAL2  =                  0.0
CTYPE2  = 'GLAT-CAR'
CUNIT2  = 'deg     '
WCSAXES =                    3
LONPOLE =                  0.0
LATPOLE =                 90.0
CRVAL3  =                  0.0
CRPIX3  =     471.921630003202
CDELT3  =     1288.21496912415
CUNIT3  = 'm/s     '
CTYPE3  = 'VRAD    '
CNAME3  = 'LSRK radio velocity'
SPECSYS = 'LSRK    '
SSYSOBS = 'LSRK    '
SSYSSRC = 'LSRK    '
VELOSYS =                  0.0
ZSOURCE =                  0.0
BUNIT   = 'K       '
BPA     =                  0.0
BMAJ    =               0.2706
BMIN    =               0.2706
EQUINOX =               2000.0
CHECKSUM= 'ZBR9f9O9ZAO9d7O9'   / HDU checksum updated 2016-09-18T12:29:42
DATASUM = '2870199985'         / data unit checksum updated 2016-09-18T12:29:42

astropy.WCS seems to parse the header ok:

from astropy.wcs import WCS
wcs = WCS(hdulist[0].header)
print(wcs)

Outputs:

WCS Keywords

Number of WCS axes: 3
CTYPE : 'GLON-CAR'  'GLAT-CAR'  'VRAD'
CRVAL : 0.0  0.0  0.0
CRPIX : 2162.0  1082.0  471.921630003202
PC1_1 PC1_2 PC1_3  : 1.0  0.0  0.0
PC2_1 PC2_2 PC2_3  : 0.0  1.0  0.0
PC3_1 PC3_2 PC3_3  : 0.0  0.0  1.0
CDELT : -0.0833333333  0.0833333333  1288.21496912415
NAXIS : 4323  2163  945

[keflavich]

Could you check whether cube[:, 1:-1, 1:-1] is any better behaved? I'll look more closely when I'm at my computer again, but one hypothesis is that we're hitting the edges of the sky in a rectangular WCS and that's breaking assumptions somewhere

[tvwenger]

Thanks for the suggestion @keflavich indeed it works:

SpectralCube with shape=(945, 2161, 4321) and unit=K:
 n_x:   4321  type_x: GLON-CAR  unit_x: deg    range:   180.000000 deg:  180.000000 deg
 n_y:   2161  type_y: GLAT-CAR  unit_y: deg    range:   -90.000000 deg:   90.000000 deg
 n_s:    945  type_s: VRAD      unit_s: m / s  range:  -606648.293 m / s:  609426.638 m / s

[keflavich]

Great, let's regard that as a useful workaround and leave this issue open until we figure out how to deal with ... closed edge cases?