spacetelescope / gwcs

Generalized World Coordinate System: provides tools for managing WCS in a general way
https://gwcs.readthedocs.io/en/latest/
46 stars 48 forks source link

Bug in APE-14 description of WCS #495

Open astrofrog opened 6 months ago

astrofrog commented 6 months ago

I'm running into issues with a solar GWCS where the high-level objects are astropy coordinates with default values of arcseconds. To reproduce this you will need the following dataset:

https://github.com/astropy/reproject/blob/main/reproject/tests/data/aia_171_level1.asdf

The issue:

In [1]: import asdf

In [3]: aia = asdf.open("/Users/tom/Code/Astropy/reproject/reproject/tests/data/aia_171_level1.asdf")

In [5]: aia['wcs'].world_axis_object_components
Out[5]: [('celestial', 0, 'spherical.lon'), ('celestial', 1, 'spherical.lat')]

In [7]: aia['wcs'].world_axis_units
Out[7]: ('deg', 'deg')

In [8]: coord = aia['wcs'].pixel_to_world(1, 2)

In [9]: coord.spherical.lon
Out[9]: <Longitude -1203.0967792 arcsec>

Basically the issue is that world_axis_object_components is incorrect because if applied to the high level objects, it will not give values that could be passed directly to world_to_pixel_values (which it should). I think the following would be more correct:

In [5]: aia['wcs'].world_axis_object_components
Out[5]: [('celestial', 0, 'spherical.lon.deg'), ('celestial', 1, 'spherical.lat.deg')]

Is this an issue in GWCS, or in the ASDF file? (I don't understand where all the layers live).

nden commented 6 months ago

I'm not well versed in the API. Which method gives a wrong result and what do you expect to get? Also how would one apply the world_axis_object_components function to high level objects?

For a spherical system the method is defined here
Ping @Cadair as he contributed to this.

astrofrog commented 6 months ago

@nden the main thing that needs to be done is that world_axis_object_components should be more specific about what units to get the values in, e.g.

In [5]: aia['wcs'].world_axis_object_components
Out[5]: [('celestial', 0, 'spherical.lon.deg'), ('celestial', 1, 'spherical.lat.deg')]

or it could be that world_axis_units should be arcsec instead of deg. I'll let @Cadair chime in though.

Cadair commented 6 months ago

I think that it's world_axis_units that's wrong. We should support things in native arcsec.

I am not working today, but will try and look next week.

Cadair commented 6 months ago

ok, I have tracked this down.

In the asdf file included in the reproject test suite some muppet (me), didn't set the units of the CelestialFrame object to u.arcsec which I should have done.

However, this brings up an interesting point, it appears to be true that the units of CelestialFrame should match the default units of the SphericalRepresentation on the reference_frame (in this case arcsec), so why don't we automatically default to that?

Cadair commented 1 day ago

I have now fixed this in #512 and also #457