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

RA not returned correctly for cubes with orthographic/synthesis projection #864

Open TimothyADavis opened 1 year ago

TimothyADavis commented 1 year ago

Dear spectral-cube team,

I have come across a fun issue where the RA coordinate grids returned by e.g. cube.world[:] are not correct when the FITS file has an orthographic/synthesis projection (e.g. RA---SIN and DEC--SIN).

For instance in my cube the RA coordinate predicted by cube.world for some given pixel is -164.004 deg, while cube.wcs.pixel_to_world returns the correct 195.819 deg.

Many thanks, Tim Davis

keflavich commented 1 year ago

Could you provide the header for said file?

TimothyADavis commented 1 year ago

Hi Adam - of course, find that attached. CARTA and astropy.wcs seem to cope with it fine.

component1.fits-Header.txt

low-sky commented 1 year ago

I think there's just a difference in convention here since the cube.world is returning results from astropy.wcs.all_pixel2world whereas you might be seeing the results from astropy.wcs.pixel_to_world. Since the header you supplied has a negative CRVAL1 the former is returning a negative value and the latter adds the 360° to put it inside the standard range. I am unable to replicate a small offset beyond this 360° shift that could be attributable to mishandling a projection since it's all astropy under the hood.

However, it might be good to move the SpectralCube.world attributes to using astropy.wcs.pixel_to_world since it comes with some units and astropy is catching up on spectral axis handling.

keflavich commented 1 year ago

We're using the lower-level routines rather than the Coordinates routines in part for performance reasons; I think there is overhead in creating a huge array of coordinates. However, we could check what the performance cost is and see if switching to pixel_to_world is acceptable.

@TimothyADavis Could you verify that there is a genuine offset, not just a wrap? (-164 = 196, but your coordinates gave a difference of ~0.17 deg)

TimothyADavis commented 1 year ago

After some investigation, it seems that it is just the wrap (I may have extracted the wrong pixel to compare before). Negative RAs just go on to break code further up my pipeline, which is why I started seeing the issue. I can always add the wrap myself- just thought I would bring it up as it seemed strange to get different results from these two routines. Cheers!

keflavich commented 1 year ago

I agree it is a bit weird. astropy coordinates have a wrap_at option to specify the convention (you want -180 < lon < 180 near lon=0, 0 < lon < 360 at lon ~ 180), but I'm not sure that that option exists for the lower-level pix2world routines.

e-koch commented 1 year ago

Thanks for raising this, @TimothyADavis.

I noticed in the header that the CRVALs are negative:

CRVAL1 = -1.641802500000E+02
CRVAL2 = -1.742302777778E+01

It'd be ideal if spectral-cube via astropy coords can handle the wrapping, but I'm wondering if anyone knows if the negative CRVALs are accepted under the FITS standard? This might be a bug upstream in WSClean and how it writes headers.

preshanth commented 1 year ago

@e-koch the fits standard only says CRVAL for a given coordinate system needs to be floating point no other specifications beyond that. Negative values are permissible but the question is does it make sense as a reference for the coordinate system?