Problem using WCS.pixel_to_world with WCS tab file

[tiagopereira]

Description

I am trying to use astropy to read a FITS file produced by the SSTRED pipeline from data taken at the Swedish 1-m Solar Telescope. The data can be read ok, but I have some trouble with the WCS object when trying to run operations such as pixel_to_world.

Expected behavior

When running w.pixel_to_world(), I expected to get an array with the world coordinates. Instead, I get an error.

Steps to Reproduce

Here is a minimal WCS header extracted from the actual object (original file is too large):

WCSAXES =                    5 / Number of coordinate axes
CRPIX1  =                  0.0 / Pixel coordinate of reference point
CRPIX2  =                  0.0 / Pixel coordinate of reference point
CRPIX3  =                  0.0 / Pixel coordinate of reference point
CRPIX4  =                  1.0 / Pixel coordinate of reference point
CRPIX5  =                  0.0 / Pixel coordinate of reference point
CDELT1  =                  1.0 / [arcsec] Coordinate increment at reference poin
CDELT2  =                  1.0 / [arcsec] Coordinate increment at reference poin
CDELT3  =                  1.0 / [nm] Coordinate increment at reference point
CDELT4  =                  1.0 / Coordinate increment at reference point
CDELT5  =                  1.0 / [s] Coordinate increment at reference point
CUNIT1  = 'arcsec'             / Units of coordinate increment and value
CUNIT2  = 'arcsec'             / Units of coordinate increment and value
CUNIT3  = 'nm'                 / Units of coordinate increment and value
CUNIT5  = 's'                  / Units of coordinate increment and value
CTYPE1  = 'HPLN-TAB'           / Coordinate type codeundefined projection
CTYPE2  = 'HPLT-TAB'           / Coordinate type codeundefined projection
CTYPE3  = 'WAVE-TAB'           / Vacuum wavelength
CTYPE4  = 'STOKES'             / Coordinate type code
CTYPE5  = 'UTC--TAB'           / Coordinate type code
CRVAL1  =                  0.0 / [arcsec] Coordinate value at reference point
CRVAL2  =                  0.0 / [arcsec] Coordinate value at reference point
CRVAL3  =                  0.0 / [nm] Coordinate value at reference point
CRVAL4  =                  1.0 / Coordinate value at reference point
CRVAL5  =                  0.0 / [s] Coordinate value at reference point
PV1_3   =                  1.0 / [deg] alias for LONPOLE (has precedence)
PV2_3   =                  2.0 /     projection parameter
PV3_3   =                  3.0 / Coordinate transformation parameter
PV5_3   =                  4.0 / Coordinate transformation parameter
PS1_0   = 'WCS-TAB '           / Coordinate transformation parameter
PS1_1   = 'HPLN+HPLT+WAVE+TIME' / Coordinate transformation parameter
PS1_2   = 'HPLN-INDEX'         / Coordinate transformation parameter
PS2_0   = 'WCS-TAB '           / Coordinate transformation parameter
PS2_1   = 'HPLN+HPLT+WAVE+TIME' / Coordinate transformation parameter
PS2_2   = 'HPLT-INDEX'         / Coordinate transformation parameter
PS3_0   = 'WCS-TAB '           / Coordinate transformation parameter
PS3_1   = 'HPLN+HPLT+WAVE+TIME' / Coordinate transformation parameter
PS5_0   = 'WCS-TAB '           / Coordinate transformation parameter
PS5_1   = 'HPLN+HPLT+WAVE+TIME' / Coordinate transformation parameter
LATPOLE =                 90.0 / [deg] Native latitude of celestial pole
CNAME1  = 'Spatial X'          / Axis name for labelling purposes
CNAME2  = 'Spatial Y'          / Axis name for labelling purposes
CNAME3  = 'Wavelength'         / Axis name for labelling purposes
CNAME5  = 'Time since DATEREF, increases with dim. 3 and 5' / Axis name for labe
CSYER1  =                 60.0 / [arcsec] Systematic error in coordinate
CSYER2  =                 60.0 / [arcsec] Systematic error in coordinate
TIMESYS = 'UTC'                / Time scale
DATEREF = '2021-06-22T00:00:00.000000' / ISO-8601 fiducial time
MJDREF  =              59387.0 / [d] MJD of fiducial time
DATE-OBS= '2021-06-22T16:10:31' / ISO-8601 time of observation
MJD-OBS =      59387.673969907 / [d] MJD of observation
DATE-BEG= '2021-06-22T16:10:33.61340' / ISO-8601 time at start of observation
MJD-BEG =      59387.674000155 / [d] MJD at start of observation
DATE-AVG= '2021-06-22T16:18:47.91310' / ISO-8601 time at midpoint of observation
MJD-AVG =      59387.679721216 / [d] MJD at midpoint of observation
DATE-END= '2021-06-22T16:27:36.28459' / ISO-8601 time at end of observation
MJD-END =      59387.685836627 / [d] MJD at end of observation
XPOSURE =             0.180015 / [s] Exposure (integration) time
OBSGEO-X=            5327403.0 / [m] observatory X-coordinate
OBSGEO-Y=           -1718726.0 / [m] observatory Y-coordinate
OBSGEO-Z=            3051730.0 / [m] observatory Z-coordinate

NOTE: with only this header the WCS will be missing the -TAB data.

If I put that into an WCS object called w, the first try is:

w.pixel_to_world(1, 1, 1, 1, 1)

This fails with an error:

NotImplementedError: SPECSYS= not yet supported

This seems to be an error with the original file, it is missing the SPECSYS keyword adding it manually, I can proceed a little longer:

w.wcs.specsys = 'TOPOCENT'
w.pixel_to_world(1, 1, 1, 1, 1)

Fails with:

ValueError: Unrecognized time CTYPE=UTC--TAB

I thought maybe the error was the double dash, so I did:

w.wcs.ctype[-1] = 'UTC-TAB'
w.pixel_to_world(1, 1, 1, 1, 1)

and then I get:

InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3074 of file cextern/wcslib/C/wcs.c:
Table parameters set for non-table axis type.

System Details

Linux-3.10.0-1127.19.1.el7.x86_64-x86_64-with-glibc2.10
Python 3.8.5 (default, Sep  4 2020, 07:30:14)
[GCC 7.3.0]
Numpy 1.20.2
astropy 4.3.1
Scipy 1.6.2
Matplotlib 3.4.1

[dhomeier]

Ah, trying to build the WCS from the header only raises the error

  File "/opt/lib/python3.9/site-packages/astropy/wcs/wcs.py", line 448, in __init__
    tmp_wcsprm = _wcs.Wcsprm(header=tmp_header_bytes, key=key,
ValueError: HDUList is required to retrieve -TAB coordinates and/or indices.

so you'll have to extract a minimal set of the data section.

Or could you try and insert in your astropy installation after the lines

astropy/astropy/wcs/wcsapi/fitswcs.py

Lines 587 to 589 in ed434a3

	# Extract time scale
	scale = self.wcs.ctype[i].lower()

this code

                     # Strip algorithm code
                     if '-' in scale:
                         scale, algorithm = scale.split('-')[0], scale.split('-')[-1]
                         warnings.warn(f'Ignoring unsupported algorithm '
                                       f'{algorithm.upper()} for scale {scale.upper()}',
                                       UserWarning)

[tiagopereira]

It took me a while but I finally managed to get one of these files in a much smaller size. I've placed them on this google drive folder.. The error should appear with both files (they have the same data, only one is transposed in the spectral/temporal dimension). In these example files the HPLN and HPLT have only 20 points each.

[Cadair]

This seems to be an issue with the HighLevelWCS layer not the wcslib layer as this works:

In [5]: w.pixel_to_world_values(0,0,0,0,0)
Out[5]: 
(array(244.35426644),
 array(377.38281263),
 array(485.92380808),
 array(1.),
 array(36157.38193))

[Cadair]

So after a little investigation, I think the issue here is that this check:

astropy/astropy/wcs/wcsapi/fitswcs.py

Lines 659 to 660 in b6352e1

	elif scale not in Time.SCALES:
	raise ValueError(f'Unrecognized time CTYPE={self.wcs.ctype[i]}')

is not agnostic to the existence of --TAB on the end of the CTYPE. I would need to go and check the standard to figure out what the correct approach here, but I think just stripping anything after - might work?

[pllim]

@mcara , have you memorized the standards by now? 😆

[Cadair]

it works 🥳

In [4]: w.pixel_to_world(0,0,0,0,0)
Out[4]: 
[<SkyCoord (Helioprojective: obstime=2021-08-04T09:56:50.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2021-08-04T09:56:50.000, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
     (0.00174495, 6.02251189, 1.51772714e+11)>): (Tx, Ty) in arcsec
     (244.35426644, 377.38281263)>,
 <SpectralCoord 
    (observer: <ITRS Coordinate (obstime=59430.4188856662, location=(0., 0., 0.) km): (x, y, z) in m
                   (5327403., -1718726., 3051730.)
                (v_x, v_y, v_z) in km / s
                   (0., 0., 0.)>
     target: <Helioprojective Coordinate (obstime=2021-08-04T09:56:50.000, rsun=695700.0 km, observer=<HeliographicStonyhurst Coordinate (obstime=2021-08-04T09:56:50.000, rsun=695700.0 km): (lon, lat, radius) in (deg, deg, m)
                 (0.00174495, 6.02251189, 1.51772714e+11)>): (Tx, Ty, distance) in (arcsec, arcsec, kpc)
                 (0., 0., 1000.)
              (d_Tx, d_Ty, d_distance) in (arcsec / s, arcsec / s, km / s)
                 (0., 0., 0.)>
     observer to target (computed from above):
       radial_velocity=-11980.853083640337 km / s
       redshift=-0.03919626429118017)
   485.92380808 nm>,
 <Quantity 1.>,
 <Time object: scale='utc' format='mjd' value=59430.41848821678>]

[mcara]

I also would like to add that currently time axis does not work well with -TAB due to a bug in WCSLIB. I am not sure what are all the consequences of the bug but selecting temporal axis using sub may not work as expected.

This was reported upstream a few months ago.

[tiagopereira]

@mcara Using astropy 5.1.1, with the update above, still gives an error with the example files, but this time it is different:

w.pixel_to_world(1, 1, 1, 1, 1)
(...)
ScaleValueError: Scale value 'UTC' not in allowed values ('tai', 'tcb', 'tcg', 'tdb', 'tt', 'ut1', 'utc', 'local')

Seems like it the value "UTC" should have been lowercased. If I try to set it to lowercase manually, I get yet a different error:

w.wcs.ctype[-1] = 'utc--tab'
w.pixel_to_world(1, 1, 1, 1, 1)
(...)
InconsistentAxisTypesError: ERROR 4 in wcs_types() at line 3021 of file cextern/wcslib/C/wcs.c:
Unrecognized projection code (tab in CTYPE5).

[Cadair]

I think I fixed that in #13571 but it needs attention.