3D WCS world coordinates are not SkyCoord objects

[dstansby]

I am trying to work with a 3D WCS, and extract SkyCoord coordinates from it. This is the code I'm using:

import astropy.wcs
import numpy as np
import astropy.units as u

data = np.random.rand(4, 3, 5)
# Define WCS transformations in an astropy WCS object.
wcs = astropy.wcs.WCS(naxis=3)
wcs.wcs.ctype = 'DSUN', 'HGLT-TAN', 'HGLN-TAN'
wcs.wcs.cunit = 'm', 'deg', 'deg'
wcs.wcs.cdelt = 1000, 0.5, 0.4
wcs.wcs.crpix = 0, 2, 2
wcs.wcs.crval = 10, 0.5, 1

print(wcs.pixel_to_world(1, 1, 1))

but this prints [<Quantity 2010. m>, <Quantity 0.5 deg>, <Quantity 1. deg>] instead of a SkyCoord object. This might be an upstream astropy issue, or an issue with our registration of coordinate systems with astropy.

[Cadair]

This has some layers to it...

1. Your example is wrong, you need to import sunpy.coordinates to register the solar coordinate frames and WCS converters. If you do this it looks like this:

import astropy.wcs
import numpy as np
import astropy.units as u

import sunpy.coordinates

data = np.random.rand(4, 3, 5)
# Define WCS transformations in an astropy WCS object.
wcs = astropy.wcs.WCS(naxis=3)
wcs.wcs.ctype = 'DSUN', 'HGLT-TAN', 'HGLN-TAN'
wcs.wcs.cunit = 'm', 'deg', 'deg'
wcs.wcs.cdelt = 1000, 0.5, 0.4
wcs.wcs.crpix = 0, 2, 2
wcs.wcs.crval = 10, 0.5, 1

print(wcs.pixel_to_world(1, 1, 1))

[<Quantity 2010. m>, <SkyCoord (HeliographicStonyhurst: obstime=None, rsun=695700.0 km): (lon, lat) in deg
    (1., 0.5)>]

closer.

2. To actually make this correctly fully 3D would take some work. The way the code is written in astropy just isn't built for it at all. The relevant section is here: https://github.com/astropy/astropy/blob/main/astropy/wcs/wcsapi/fitswcs.py#L382-L399

I am not sure how we could teach astropy that DSUN is a valid 3rd argument to the SkyCoord constructor without subclassing WCS (something I have strived long and hard to avoid). I think we probably need to open an issue upstream to discuss this and see what people say.

[dstansby]

Thanks for first fix. On the second point, I think it might even be modifications to fitswcs that are needed. The motivation for doing this isn't tied to FITS, so I shall go and take a look at other options (I guess gwcs?) for what I want.

[dstansby]

I managed to make a gwcs 🎉 , so going to close this.

[nabobalis]

Worth an example maybe?

[dstansby]

Maybe... this code I'm writing will at least end up in pfsspy, so once that's sorted I'll have a think about where best to put an example.

[Cadair]

I think this is worth keeping open, if you read Thompson and the FITS specs this should work, so it would be good to be able to support it.