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
95 stars 61 forks source link

Spatial reproject not working for cubes with different rest frequencies #874

Open SpacialTree opened 1 year ago

SpacialTree commented 1 year ago

I am attempting to spatially reproject one cube to another.

# cube17.header['RESTFRQ] = 230538000000.0
# cube_29.header['RESTFRQ] = 220398684200.0

cube17.reproject(cube_29.header)

I am getting the following error:

ValueError                                Traceback (most recent call last)
/scratch/local/60838376/ipykernel_89087/2137793806.py in <module>
----> 1 reproj_cube17 = interp_cube17.reproject(cube_29.header)

/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/utils.py in wrapper(self, *args, **kwargs)
     47                           PossiblySlowWarning
     48                          )
---> 49         return function(self, *args, **kwargs)
     50     return wrapper
     51 

/blue/adamginsburg/adamginsburg/repos/spectral-cube/spectral_cube/spectral_cube.py in reproject(self, header, order, use_memmap, filled, **kwargs)
   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 "

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 issue seems to be that reproject is trying to take the frequency space of the cubes into account when reprojecting, making the reprojected cube be full of NaNs.

A temporary fix is to make a copy of the target header and set its rest frequency to that of the cube you are trying to reproject.

tgthdr = cube_29.header
tgthdr['RESTFRQ'] = cube_17.header['RESTFRQ'], 'HACKHACKHACK'
reproj_cube17 = cube_17.reproject(tgthdr)
keflavich commented 1 year ago

this appears to be a recent change in how spectral coordinate transforms are handled by astropy. I haven't found the origin of the issue.

keflavich commented 1 month ago

I'm officially endorsing the hackaround. @astrofrog do you happen to have any idea what caused this change? It was something in astropy that's getting called by reproject. The WCS is trying to convert back to frequency even though we're trying to get world coordinates in velocity, so all the hard work we did to get into matched velocity coordinates is getting undone. It's not easy to track exactly where this happens, though.

astrofrog commented 1 month ago

It's probably related to the fact that reproject uses pixel_to_pixel in astropy which now uses SpectralCoord to do conversions in world coordinates for the spectral axis. I think it would be good to have the option of not converting back to frequency space though.