radio-astro-tools / spectral-cube

Library for reading and analyzing astrophysical spectral data cubes
http://spectral-cube.rtfd.org
BSD 3-Clause "New" or "Revised" License
98 stars 65 forks source link

Cube reprojection tutorial failing on learn.astropy.org #920

Open jeffjennings opened 1 month ago

jeffjennings commented 1 month ago

Hi @adamginsburg and @e-koch, I'm hoping for help with fixing the cube reprojection tutorial on learn.astropy.org. The Build tutorials workflow in astropy-tutorials has been failing at least since last year, and the cube reprojection tutorial notebook now errors when executing, presumably due to using a newer version of spectral-cube than when the notebook was written. The error occurs in the second cell of Step 6 in the notebook at

cube2vel_reproj = cube2vel_spatialspectralsmooth.reproject(cube1vel.header)

with traceback

ValueError                                Traceback (most recent call last)
Cell In[81], line 1
----> 1 cube2vel_reproj = cube2vel_spatialspectralsmooth.reproject(cube1vel.header)
      2 cube2vel_reproj

File /opt/miniconda3/envs/astropy-tutorials/lib/python3.12/site-packages/spectral_cube/utils.py:49, in warn_slow.<locals>.wrapper(self, *args, **kwargs)
     43 elif warn_how and not self._is_huge:
     44     # TODO: add check for whether cube has been loaded into memory
     45     warnings.warn(\"This function ({0}) requires loading the entire cube into \"
     46                   \"memory and may therefore be slow.\".format(str(function)),
     47                   PossiblySlowWarning
     48                  )
---> 49 return function(self, *args, **kwargs)

File /opt/miniconda3/envs/astropy-tutorials/lib/python3.12/site-packages/spectral_cube/spectral_cube.py:2705, in BaseSpectralCube.reproject(self, header, order, use_memmap, filled, **kwargs)
   2697 newcube, newcube_valid = reproject_interp((data,
   2698                                            self.header),
   2699                                           newwcs,
   (...)
   2702                                           order=order,
   2703                                           **reproj_kwargs)
   2704 if np.all(np.isnan(newcube)):
-> 2705     raise ValueError(\"All values in reprojected cube are nan.  This can be caused\"
   2706                      \" by an error in which coordinates do not 'round-trip'.  Try \"
   2707                      \"setting ``roundtrip_coords=False``.  You might also check \"
   2708                      \"whether the WCS transformation produces valid pixel->world \"
   2709                      \"and world->pixel coordinates in each axis.\"
   2710                      )
   2712 return self._new_cube_with(data=newcube,
   2713                            wcs=newwcs,
   2714                            mask=BooleanArrayMask(newcube_valid.astype('bool'),
   2715                                                  newwcs),
   2716                            meta=self.meta,
   2717                           )

ValueError: All values in reprojected cube are nan.  This can be caused by an error in which coordinates do not 'round-trip'.  Try setting ``roundtrip_coords=False``.  You might also check whether the WCS transformation produces valid pixel->world and world->pixel coordinates in each axis."
}

The input cube, cube2vel_spatialspectralsmooth, is not all nan before the projection. I tried adding roundtrip_coords=False to the call as suggested, and checked the header used in the call with

wcs.WCS(cube1vel.header)

which returns

WCS Keywords

Number of WCS axes: 3
CTYPE : 'RA---SIN' 'DEC--SIN' 'VRAD' 
CRVAL : 266.5442833334 -28.70493055556 -43509.079452344 
CRPIX : 126.0 126.0 1.0 
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 : -7.222222222222e-05 7.222222222222e-05 2002.6282784629 
NAXIS : 250  250  75

Also for reference,

wcs.WCS(cube2vel_spatialspectralsmooth.header)

returns

WCS Keywords

Number of WCS axes: 3
CTYPE : 'RA---SIN' 'DEC--SIN' 'VRAD' 
CRVAL : 266.5442833333 -28.704955 -43509.079452344 
CRPIX : 211.0 211.0 1.0 
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 : -3.055555555556e-05 3.055555555556e-05 2002.6282784629 
NAXIS : 420  420  75

and

cube2vel_spatialspectralsmooth

is

SpectralCube with shape=(75, 420, 420) and unit=Jy / beam:
 n_x:    420  type_x: RA---SIN  unit_x: deg    range:   266.537002 deg:  266.551600 deg
 n_y:    420  type_y: DEC--SIN  unit_y: deg    range:   -28.711371 deg:  -28.698569 deg
 n_s:     75  type_s: VRAD      unit_s: km / s  range:      -43.509 km / s:     104.685 km / s

Any advice or fixes you have would be helpful, thanks!

keflavich commented 1 month ago

Just to verify - changing roundtrip_coords=False did not fix the problem?

jeffjennings commented 1 month ago

That's right, I still get the same error when using roundtrip_coords=False.

e-koch commented 1 month ago

This is (I suspect) the same as problem as in #874 : SpectralCoord is reverting back to frequencies whenever RESTFRQ is defined in the header.

We could hack in something in spectral-cube that avoids the transformation back to freq; or an option upstream in astropy that avoids this by default

jeffjennings commented 1 month ago

Yes you're right; I tried the workaround in #874 - it worked and the reprojection looks correct. Thank you!

For the moment I'm just adding the workaround in the tutorial to get the notebook to build. But I can leave this issue open if you'd like and revert the workaround in the tutorial once there's a fix here or upstream.