From bf1b98b305de176c724052bb8cf4223b392e695f Mon Sep 17 00:00:00 2001 From: Stuart Mumford Date: Wed, 24 Aug 2022 15:13:10 +0100 Subject: [PATCH] Fix a whole bunch of bugs todo with -TAB in wcsapi/fitswcs --- astropy/wcs/wcsapi/fitswcs.py | 12 ++++--- astropy/wcs/wcsapi/tests/example_4d_tab.fits | 1 + astropy/wcs/wcsapi/tests/test_fitswcs.py | 35 ++++++++++++++++--- .../wcsapi/wrappers/tests/test_sliced_wcs.py | 22 ++++++------ docs/changes/wcs/13571.bugfix.rst | 2 ++ setup.cfg | 2 +- 6 files changed, 53 insertions(+), 21 deletions(-) create mode 100644 astropy/wcs/wcsapi/tests/example_4d_tab.fits create mode 100644 docs/changes/wcs/13571.bugfix.rst diff --git a/astropy/wcs/wcsapi/fitswcs.py b/astropy/wcs/wcsapi/fitswcs.py index 30c87577ce8c..1f3b5b6ea43b 100644 --- a/astropy/wcs/wcsapi/fitswcs.py +++ b/astropy/wcs/wcsapi/fitswcs.py @@ -416,7 +416,11 @@ def _get_components_and_classes(self): else: kwargs = {} kwargs["frame"] = celestial_frame - kwargs["unit"] = u.deg + # Very occasionally (i.e. with TAB) wcs does not convert the units to degrees + kwargs["unit"] = ( + u.Unit(self.wcs.cunit[self.wcs.lng]), + u.Unit(self.wcs.cunit[self.wcs.lat]), + ) classes["celestial"] = (SkyCoord, (), kwargs) @@ -445,7 +449,7 @@ def _get_components_and_classes(self): earth_location = EarthLocation(*self.wcs.obsgeo[:3], unit=u.m) # Get the time scale from TIMESYS or fall back to 'utc' - tscale = self.wcs.timesys or "utc" + tscale = self.wcs.timesys.lower() or "utc" if np.isnan(self.wcs.mjdavg): obstime = Time( @@ -695,8 +699,8 @@ def value_from_spectralcoord(spectralcoord): # Initialize delta reference_time_delta = None - # Extract time scale - scale = self.wcs.ctype[i].lower() + # Extract time scale, and remove any algorithm code + scale = self.wcs.ctype[i].split("-")[0].lower() if scale == "time": if self.wcs.timesys: diff --git a/astropy/wcs/wcsapi/tests/example_4d_tab.fits b/astropy/wcs/wcsapi/tests/example_4d_tab.fits new file mode 100644 index 000000000000..4776304227d2 --- /dev/null +++ b/astropy/wcs/wcsapi/tests/example_4d_tab.fits @@ -0,0 +1 @@ +SIMPLE = T / conforms to FITS standard BITPIX = 8 / array data type NAXIS = 0 / number of array dimensions EXTEND = T WCSAXES = 4 WCSNAME = 'CompositeFrame' CTYPE1 = 'GLON-TAB' CUNIT1 = 'arcsec ' PS1_0 = 'WCS-TABLE' PV1_1 = 1 PS1_1 = 'coordinates' PV1_3 = 1 CRVAL1 = 1 CRPIX1 = 1.0 PC1_1 = 1.0 CDELT1 = 1.0 CTYPE2 = 'GLAT-TAB' CUNIT2 = 'arcsec ' PS2_0 = 'WCS-TABLE' PV2_1 = 1 PS2_1 = 'coordinates' PV2_3 = 2 CRVAL2 = 1 CRPIX2 = 1.0 PC2_2 = 1.0 CDELT2 = 1.0 CTYPE3 = 'FREQ-TAB' CUNIT3 = 'Hz ' PS3_0 = 'WCS-TABLE' PV3_1 = 1 PS3_1 = 'coordinates' PV3_3 = 3 CRVAL3 = 1 CRPIX3 = 1.0 PC3_3 = 1.0 CDELT3 = 1.0 CTYPE4 = 'TIME-TAB' CUNIT4 = '' PS4_0 = 'WCS-TABLE' PV4_1 = 1 PS4_1 = 'coordinates' PV4_3 = 4 CRVAL4 = 1 CRPIX4 = 1.0 PC4_4 = 1.0 CDELT4 = 1.0 TIMESYS = 'UTC ' END XTENSION= 'BINTABLE' / binary table extension BITPIX = 8 / array data type NAXIS = 2 / number of array dimensions NAXIS1 = 3840 / length of dimension 1 NAXIS2 = 1 / length of dimension 2 PCOUNT = 0 / number of group parameters GCOUNT = 1 / number of groups TFIELDS = 1 / number of table fields TTYPE1 = 'coordinates' TFORM1 = '480D ' TDIM1 = '(4,5,4,3,2)' EXTNAME = 'WCS-TABLE' / extension name EXTVER = 1 / extension value END @m#Ȏ ,U]@8\(@<@j؉=V \@8\(@<@fV @8\(@<@bJ'vGHV \@8\(@<@_nYU]@8\(@<@nm]U"@8\(@<@l@+V%~ކ@8\(@<@fV@@8\(@<@`#2V%~ޅ@8\(@<@]$2WD!U"@8\(@<@oUԽ@8\(@<@nmqmhV8rQ/@8\(@<@fV`@8\(@<@]$Fe%$V8rQ/@8\(@<@ZHWUԼ@8\(@<@pV@8\(@<@pV@@8\(@<@FV@8\(@<@VV@@8\(@<@VV@8\(@<@m#Ȏ ,U]@8Q@<@j؉=V \@8Q@<@fV @8Q@<@bJ'vGHV \@8Q@<@_nYU]@8Q@<@nm]U"@8Q@<@l@+V%~ކ@8Q@<@fV@@8Q@<@`#2V%~ޅ@8Q@<@]$2WD!U"@8Q@<@oUԽ@8Q@<@nmqmhV8rQ/@8Q@<@fV`@8Q@<@]$Fe%$V8rQ/@8Q@<@ZHWUԼ@8Q@<@pV@8Q@<@pV@@8Q@<@FV@8Q@<@VV@@8Q@<@VV@8Q@<@m#Ȏ ,U]@8G{@<@j؉=V \@8G{@<@fV @8G{@<@bJ'vGHV \@8G{@<@_nYU]@8G{@<@nm]U"@8G{@<@l@+V%~ކ@8G{@<@fV@@8G{@<@`#2V%~ޅ@8G{@<@]$2WD!U"@8G{@<@oUԽ@8G{@<@nmqmhV8rQ/@8G{@<@fV`@8G{@<@]$Fe%$V8rQ/@8G{@<@ZHWUԼ@8G{@<@pV@8G{@<@pV@@8G{@<@FV@8G{@<@VV@@8G{@<@VV@8G{@<@m#Ȏ ,U]@8\(@<@j؉=V \@8\(@<@fV @8\(@<@bJ'vGHV \@8\(@<@_nYU]@8\(@<@nm]U"@8\(@<@l@+V%~ކ@8\(@<@fV@@8\(@<@`#2V%~ޅ@8\(@<@]$2WD!U"@8\(@<@oUԽ@8\(@<@nmqmhV8rQ/@8\(@<@fV`@8\(@<@]$Fe%$V8rQ/@8\(@<@ZHWUԼ@8\(@<@pV@8\(@<@pV@@8\(@<@FV@8\(@<@VV@@8\(@<@VV@8\(@<@m#Ȏ ,U]@8Q@<@j؉=V \@8Q@<@fV @8Q@<@bJ'vGHV \@8Q@<@_nYU]@8Q@<@nm]U"@8Q@<@l@+V%~ކ@8Q@<@fV@@8Q@<@`#2V%~ޅ@8Q@<@]$2WD!U"@8Q@<@oUԽ@8Q@<@nmqmhV8rQ/@8Q@<@fV`@8Q@<@]$Fe%$V8rQ/@8Q@<@ZHWUԼ@8Q@<@pV@8Q@<@pV@@8Q@<@FV@8Q@<@VV@@8Q@<@VV@8Q@<@m#Ȏ ,U]@8G{@<@j؉=V \@8G{@<@fV @8G{@<@bJ'vGHV \@8G{@<@_nYU]@8G{@<@nm]U"@8G{@<@l@+V%~ކ@8G{@<@fV@@8G{@<@`#2V%~ޅ@8G{@<@]$2WD!U"@8G{@<@oUԽ@8G{@<@nmqmhV8rQ/@8G{@<@fV`@8G{@<@]$Fe%$V8rQ/@8G{@<@ZHWUԼ@8G{@<@pV@8G{@<@pV@@8G{@<@FV@8G{@<@VV@@8G{@<@VV@8G{@< \ No newline at end of file diff --git a/astropy/wcs/wcsapi/tests/test_fitswcs.py b/astropy/wcs/wcsapi/tests/test_fitswcs.py index e0c37644faf0..36c46911a868 100644 --- a/astropy/wcs/wcsapi/tests/test_fitswcs.py +++ b/astropy/wcs/wcsapi/tests/test_fitswcs.py @@ -4,6 +4,7 @@ import warnings from itertools import product +from pathlib import Path import numpy as np import pytest @@ -20,6 +21,7 @@ SkyCoord, SpectralCoord, ) +from astropy.io import fits from astropy.io.fits import Header from astropy.io.fits.verify import VerifyWarning from astropy.tests.helper import assert_quantity_allclose @@ -145,7 +147,7 @@ def test_simple_celestial(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], ICRS) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert_allclose(wcs.pixel_to_world_values(29, 39), (10, 20)) assert_allclose(wcs.array_index_to_world_values(39, 29), (10, 20)) @@ -275,7 +277,7 @@ def test_spectral_cube(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -395,7 +397,7 @@ def test_spectral_cube_nonaligned(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -491,7 +493,7 @@ def test_time_cube(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], ICRS) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["time"][0] is Time assert wcs.world_axis_object_classes["time"][1] == () @@ -905,7 +907,7 @@ def test_caching_components_and_classes(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], ICRS) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) wcs.wcs.radesys = "FK5" @@ -1348,3 +1350,26 @@ def check_wcs(header): # Check fall back to scale='utc' del hdr["TIMESYS"] check_wcs(hdr) + + +def test_fits_tab(): + """ + This test is a regression test for https://github.com/astropy/astropy/issues/12095 + + It checks the following: + - If a spatial WCS isn't converted to units of deg by wcslib it still works. + - If TIMESYS is upper case we parse it correctly + - If a TIME CTYPE key has an algorithm code (in this case -TAB) it still works. + + The file used here was generated by gWCS and then edited to add the TIMESYS key. + """ + with fits.open(Path(__file__).parent / "example_4d_tab.fits") as hdul: + with pytest.warns(FITSFixedWarning): + w = WCS(header=hdul[0].header, fobj=hdul) + world = w.pixel_to_world(0, 0, 0, 0) + assert isinstance(world[0], SkyCoord) + assert world[0].data.lat.unit is u.arcsec + assert world[0].data.lon.unit is u.arcsec + assert isinstance(world[1], u.Quantity) + assert isinstance(world[2], Time) + assert world[2].scale == "utc" diff --git a/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py b/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py index 623cd2e5c721..861dbef8b3e6 100644 --- a/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py +++ b/astropy/wcs/wcsapi/wrappers/tests/test_sliced_wcs.py @@ -146,7 +146,7 @@ def test_ellipsis(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -227,7 +227,7 @@ def test_spectral_slice(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert_allclose(wcs.pixel_to_world_values(29, 44), (10, 25)) assert_allclose(wcs.array_index_to_world_values(44, 29), (10, 25)) @@ -305,7 +305,7 @@ def test_spectral_range(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -385,7 +385,7 @@ def test_celestial_slice(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -467,7 +467,7 @@ def test_celestial_range(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -558,7 +558,7 @@ def test_celestial_range_rot(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -661,7 +661,7 @@ def test_no_array_shape(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wcs.world_axis_object_classes["spectral"][0] is Quantity assert wcs.world_axis_object_classes["spectral"][1] == () @@ -758,7 +758,7 @@ def test_ellipsis_none_types(): assert wcs.world_axis_object_classes["celestial"][0] is SkyCoord assert wcs.world_axis_object_classes["celestial"][1] == () assert isinstance(wcs.world_axis_object_classes["celestial"][2]["frame"], Galactic) - assert wcs.world_axis_object_classes["celestial"][2]["unit"] is u.deg + assert wcs.world_axis_object_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert_allclose(wcs.pixel_to_world_values(29, 39, 44), (10, 20, 25)) assert_allclose(wcs.array_index_to_world_values(44, 39, 29), (10, 20, 25)) @@ -951,7 +951,7 @@ def test_dropped_dimensions(): assert wao_classes["celestial"][0] is SkyCoord assert wao_classes["celestial"][1] == () assert isinstance(wao_classes["celestial"][2]["frame"], Galactic) - assert wao_classes["celestial"][2]["unit"] is u.deg + assert wao_classes["celestial"][2]["unit"] == (u.deg, u.deg) sub = SlicedLowLevelWCS(wcs, np.s_[5, :5, 12]) @@ -975,7 +975,7 @@ def test_dropped_dimensions(): assert wao_classes["celestial"][0] is SkyCoord assert wao_classes["celestial"][1] == () assert isinstance(wao_classes["celestial"][2]["frame"], Galactic) - assert wao_classes["celestial"][2]["unit"] is u.deg + assert wao_classes["celestial"][2]["unit"] == (u.deg, u.deg) def test_dropped_dimensions_4d(cube_4d_fitswcs): @@ -999,7 +999,7 @@ def test_dropped_dimensions_4d(cube_4d_fitswcs): assert wao_classes["celestial"][0] is SkyCoord assert wao_classes["celestial"][1] == () assert isinstance(wao_classes["celestial"][2]["frame"], ICRS) - assert wao_classes["celestial"][2]["unit"] is u.deg + assert wao_classes["celestial"][2]["unit"] == (u.deg, u.deg) assert wao_classes["spectral"][0:3] == (u.Quantity, (), {}) assert wao_components[0] == ("celestial", 0, "spherical.lon.degree") diff --git a/docs/changes/wcs/13571.bugfix.rst b/docs/changes/wcs/13571.bugfix.rst new file mode 100644 index 000000000000..26225953bfe1 --- /dev/null +++ b/docs/changes/wcs/13571.bugfix.rst @@ -0,0 +1,2 @@ +Fix bugs with high level WCS API on ``wcs.WCS`` object when using ``-TAB`` +coordinates. diff --git a/setup.cfg b/setup.cfg index a4b4a7e5b1ee..d220cda056fe 100644 --- a/setup.cfg +++ b/setup.cfg @@ -112,7 +112,7 @@ astropy = astropy.cfg, CITATION astropy.cosmology = *.ecsv astropy.utils.tests = data/.hidden_file.txt astropy.tests.figures = *.json -astropy.wcs = include/*/*.h +astropy.wcs = include/*/*.h, wcsapi/tests/example_4d_tab.fits astropy.wcs.tests = extension/*.c [tool:pytest]