File: C:/Users/fred/anaconda3/Lib/site-packages/astropy/wcs/wcsapi/tests/test_fitswcs.py
# Note that we test the main astropy.wcs.WCS class directly rather than testing
# the mix-in class on its own (since it's not functional without being used as
# a mix-in)
import warnings
from contextlib import nullcontext
from itertools import product
import numpy as np
import pytest
from numpy.testing import assert_allclose, assert_array_equal, assert_equal
from packaging.version import Version
from astropy import units as u
from astropy.coordinates import (
FK5,
ICRS,
ITRS,
EarthLocation,
Galactic,
SkyCoord,
SpectralCoord,
StokesCoord,
)
from astropy.io import fits
from astropy.io.fits import Header
from astropy.io.fits.verify import VerifyWarning
from astropy.tests.helper import PYTEST_LT_8_0, assert_quantity_allclose
from astropy.time import Time
from astropy.units import Quantity
from astropy.units.core import UnitsWarning
from astropy.utils import iers
from astropy.utils.data import get_pkg_data_filename
from astropy.utils.exceptions import AstropyUserWarning
from astropy.wcs._wcs import __version__ as wcsver
from astropy.wcs.wcs import WCS, FITSFixedWarning, NoConvergence, Sip
from astropy.wcs.wcsapi.fitswcs import VELOCITY_FRAMES, custom_ctype_to_ucd_mapping
###############################################################################
# The following example is the simplest WCS with default values
###############################################################################
WCS_EMPTY = WCS(naxis=1)
WCS_EMPTY.wcs.crpix = [1]
def test_empty():
wcs = WCS_EMPTY
# Low-level API
assert wcs.pixel_n_dim == 1
assert wcs.world_n_dim == 1
assert wcs.array_shape is None
assert wcs.pixel_shape is None
assert wcs.world_axis_physical_types == [None]
assert wcs.world_axis_units == [""]
assert wcs.pixel_axis_names == [""]
assert wcs.world_axis_names == [""]
assert_equal(wcs.axis_correlation_matrix, True)
assert wcs.world_axis_object_components == [("world", 0, "value")]
assert wcs.world_axis_object_classes["world"][0] is Quantity
assert wcs.world_axis_object_classes["world"][1] == ()
assert wcs.world_axis_object_classes["world"][2]["unit"] is u.one
assert_allclose(wcs.pixel_to_world_values(29), 29)
assert_allclose(wcs.array_index_to_world_values(29), 29)
assert np.ndim(wcs.pixel_to_world_values(29)) == 0
assert np.ndim(wcs.array_index_to_world_values(29)) == 0
assert_allclose(wcs.world_to_pixel_values(29), 29)
assert_equal(wcs.world_to_array_index_values(29), (29,))
assert np.ndim(wcs.world_to_pixel_values(29)) == 0
assert np.ndim(wcs.world_to_array_index_values(29)) == 0
# High-level API
coord = wcs.pixel_to_world(29)
assert_quantity_allclose(coord, 29 * u.one)
assert np.ndim(coord) == 0
coord = wcs.array_index_to_world(29)
assert_quantity_allclose(coord, 29 * u.one)
assert np.ndim(coord) == 0
coord = 15 * u.one
x = wcs.world_to_pixel(coord)
assert_allclose(x, 15.0)
assert np.ndim(x) == 0
i = wcs.world_to_array_index(coord)
assert_equal(i, 15)
assert np.ndim(i) == 0
###############################################################################
# The following example is a simple 2D image with celestial coordinates
###############################################################################
HEADER_SIMPLE_CELESTIAL = """
WCSAXES = 2
CTYPE1 = RA---TAN
CTYPE2 = DEC--TAN
CRVAL1 = 10
CRVAL2 = 20
CRPIX1 = 30
CRPIX2 = 40
CDELT1 = -0.1
CDELT2 = 0.1
CROTA2 = 0.
CUNIT1 = deg
CUNIT2 = deg
"""
with warnings.catch_warnings():
warnings.simplefilter("ignore", VerifyWarning)
WCS_SIMPLE_CELESTIAL = WCS(Header.fromstring(HEADER_SIMPLE_CELESTIAL, sep="\n"))
def test_simple_celestial():
wcs = WCS_SIMPLE_CELESTIAL
# Low-level API
assert wcs.pixel_n_dim == 2
assert wcs.world_n_dim == 2
assert wcs.array_shape is None
assert wcs.pixel_shape is None
assert wcs.world_axis_physical_types == ["pos.eq.ra", "pos.eq.dec"]
assert wcs.world_axis_units == ["deg", "deg"]
assert wcs.pixel_axis_names == ["", ""]
assert wcs.world_axis_names == ["", ""]
assert_equal(wcs.axis_correlation_matrix, True)
assert wcs.world_axis_object_components == [
("celestial", 0, "spherical.lon.degree"),
("celestial", 1, "spherical.lat.degree"),
]
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"] == (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))
assert_allclose(wcs.world_to_pixel_values(10, 20), (29.0, 39.0))
assert_equal(wcs.world_to_array_index_values(10, 20), (39, 29))
# High-level API
coord = wcs.pixel_to_world(29, 39)
assert isinstance(coord, SkyCoord)
assert isinstance(coord.frame, ICRS)
assert_allclose(coord.ra.deg, 10)
assert_allclose(coord.dec.deg, 20)
coord = wcs.array_index_to_world(39, 29)
assert isinstance(coord, SkyCoord)
assert isinstance(coord.frame, ICRS)
assert_allclose(coord.ra.deg, 10)
assert_allclose(coord.dec.deg, 20)
coord = SkyCoord(10, 20, unit="deg", frame="icrs")
x, y = wcs.world_to_pixel(coord)
assert_allclose(x, 29.0)
assert_allclose(y, 39.0)
i, j = wcs.world_to_array_index(coord)
assert_equal(i, 39)
assert_equal(j, 29)
# Check that if the coordinates are passed in a different frame things still
# work properly
coord_galactic = coord.galactic
x, y = wcs.world_to_pixel(coord_galactic)
assert_allclose(x, 29.0)
assert_allclose(y, 39.0)
i, j = wcs.world_to_array_index(coord_galactic)
assert_equal(i, 39)
assert_equal(j, 29)
# Check that we can actually index the array
data = np.arange(3600).reshape((60, 60))
coord = SkyCoord(10, 20, unit="deg", frame="icrs")
index = wcs.world_to_array_index(coord)
assert_equal(data[index], 2369)
coord = SkyCoord([10, 12], [20, 22], unit="deg", frame="icrs")
index = wcs.world_to_array_index(coord)
assert_equal(data[index], [2369, 3550])
###############################################################################
# The following example is a spectral cube with axes in an unusual order
###############################################################################
HEADER_SPECTRAL_CUBE = """
WCSAXES = 3
CTYPE1 = GLAT-CAR
CTYPE2 = FREQ
CTYPE3 = GLON-CAR
CNAME1 = Latitude
CNAME2 = Frequency
CNAME3 = Longitude
CRVAL1 = 10
CRVAL2 = 20
CRVAL3 = 25
CRPIX1 = 30
CRPIX2 = 40
CRPIX3 = 45
CDELT1 = -0.1
CDELT2 = 0.5
CDELT3 = 0.1
CUNIT1 = deg
CUNIT2 = Hz
CUNIT3 = deg
"""
with warnings.catch_warnings():
warnings.simplefilter("ignore", VerifyWarning)
WCS_SPECTRAL_CUBE = WCS(Header.fromstring(HEADER_SPECTRAL_CUBE, sep="\n"))
def test_spectral_cube():
# Spectral cube with a weird axis ordering
wcs = WCS_SPECTRAL_CUBE
# Low-level API
assert wcs.pixel_n_dim == 3
assert wcs.world_n_dim == 3
assert wcs.array_shape is None
assert wcs.pixel_shape is None
assert wcs.world_axis_physical_types == [
"pos.galactic.lat",
"em.freq",
"pos.galactic.lon",
]
assert wcs.world_axis_units == ["deg", "Hz", "deg"]
assert wcs.pixel_axis_names == ["", "", ""]
assert wcs.world_axis_names == ["Latitude", "Frequency", "Longitude"]
assert_equal(
wcs.axis_correlation_matrix,
[[True, False, True], [False, True, False], [True, False, True]],
)
assert len(wcs.world_axis_object_components) == 3
assert wcs.world_axis_object_components[0] == (
"celestial",
1,
"spherical.lat.degree",
)
assert wcs.world_axis_object_components[1][:2] == ("spectral", 0)
assert wcs.world_axis_object_components[2] == (
"celestial",
0,
"spherical.lon.degree",
)
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"] == (u.deg, u.deg)
assert wcs.world_axis_object_classes["spectral"][0] is Quantity
assert wcs.world_axis_object_classes["spectral"][1] == ()
assert wcs.world_axis_object_classes["spectral"][2] == {}
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))
assert_allclose(wcs.world_to_pixel_values(10, 20, 25), (29.0, 39.0, 44.0))
assert_equal(wcs.world_to_array_index_values(10, 20, 25), (44, 39, 29))
# High-level API
coord, spec = wcs.pixel_to_world(29, 39, 44)
assert isinstance(coord, SkyCoord)
assert isinstance(coord.frame, Galactic)
assert_allclose(coord.l.deg, 25)
assert_allclose(coord.b.deg, 10)
assert isinstance(spec, SpectralCoord)
assert_allclose(spec.to_value(u.Hz), 20)
coord, spec = wcs.array_index_to_world(44, 39, 29)
assert isinstance(coord, SkyCoord)
assert isinstance(coord.frame, Galactic)
assert_allclose(coord.l.deg, 25)
assert_allclose(coord.b.deg, 10)
assert isinstance(spec, SpectralCoord)
assert_allclose(spec.to_value(u.Hz), 20)
coord = SkyCoord(25, 10, unit="deg", frame="galactic")
spec = 20 * u.Hz
with pytest.warns(AstropyUserWarning, match="No observer defined on WCS"):
x, y, z = wcs.world_to_pixel(coord, spec)
assert_allclose(x, 29.0)
assert_allclose(y, 39.0)
assert_allclose(z, 44.0)
# Order of world coordinates shouldn't matter
with pytest.warns(AstropyUserWarning, match="No observer defined on WCS"):
x, y, z = wcs.world_to_pixel(spec, coord)
assert_allclose(x, 29.0)
assert_allclose(y, 39.0)
assert_allclose(z, 44.0)
with pytest.warns(AstropyUserWarning, match="No observer defined on WCS"):
i, j, k = wcs.world_to_array_index(coord, spec)
assert_equal(i, 44)
assert_equal(j, 39)
assert_equal(k, 29)
# Order of world coordinates shouldn't matter
with pytest.warns(AstropyUserWarning, match="No observer defined on WCS"):
i, j, k = wcs.world_to_array_index(spec, coord)
assert_equal(i, 44)
assert_equal(j, 39)
assert_equal(k, 29)
HEADER_SPECTRAL_CUBE_NONALIGNED = (
HEADER_SPECTRAL_CUBE.strip()
+ "\n"
+ """
PC2_3 = -0.5
PC3_2 = +0.5
"""
)
with warnings.catch_warnings():
warnings.simplefilter("ignore", VerifyWarning)
WCS_SPECTRAL_CUBE_NONALIGNED = WCS(
Header.fromstring(HEADER_SPECTRAL_CUBE_NONALIGNED, sep="\n")
)
def test_spectral_cube_nonaligned():
# Make sure that correlation matrix gets adjusted if there are non-identity
# CD matrix terms.
wcs = WCS_SPECTRAL_CUBE_NONALIGNED
assert wcs.world_axis_physical_types == [
"pos.galactic.lat",
"em.freq",
"pos.galactic.lon",
]
assert wcs.world_axis_units == ["deg", "Hz", "deg"]
assert wcs.pixel_axis_names == ["", "", ""]
assert wcs.world_axis_names == ["Latitude", "Frequency", "Longitude"]
assert_equal(
wcs.axis_correlation_matrix,
[
[True, True, True],
[False, True, True],
[True, True, True],
],
)
# NOTE: we check world_axis_object_components and world_axis_object_classes
# again here because in the past this failed when non-aligned axes were
# present, so this serves as a regression test.
assert len(wcs.world_axis_object_components) == 3
assert wcs.world_axis_object_components[0] == (
"celestial",
1,
"spherical.lat.degree",
)
assert wcs.world_axis_object_components[1][:2] == ("spectral", 0)
assert wcs.world_axis_object_components[2] == (
"celestial",
0,
"spherical.lon.degree",
)
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"] == (u.deg, u.deg)
assert wcs.world_axis_object_classes["spectral"][0] is Quantity
assert wcs.world_axis_object_classes["spectral"][1] == ()
assert wcs.world_axis_object_classes["spectral"][2] == {}
###############################################################################
# The following example is from Rots et al (2015), Table 5. It represents a
# cube with two spatial dimensions and one time dimension
###############################################################################
HEADER_TIME_CUBE = """
SIMPLE = T / Fits standard
BITPIX = -32 / Bits per pixel
NAXIS = 3 / Number of axes
NAXIS1 = 2048 / Axis length
NAXIS2 = 2048 / Axis length
NAXIS3 = 11 / Axis length
DATE = '2008-10-28T14:39:06' / Date FITS file was generated
OBJECT = '2008 TC3' / Name of the object observed
EXPTIME = 1.0011 / Integration time
MJD-OBS = 54746.02749237 / Obs start
DATE-OBS= '2008-10-07T00:39:35.3342' / Observing date
TELESCOP= 'VISTA' / ESO Telescope Name
INSTRUME= 'VIRCAM' / Instrument used.
TIMESYS = 'UTC' / From Observatory Time System
TREFPOS = 'TOPOCENT' / Topocentric
MJDREF = 54746.0 / Time reference point in MJD
RADESYS = 'ICRS' / Not equinoctal
CTYPE2 = 'RA---ZPN' / Zenithal Polynomial Projection
CRVAL2 = 2.01824372640628 / RA at ref pixel
CUNIT2 = 'deg' / Angles are degrees always
CRPIX2 = 2956.6 / Pixel coordinate at ref point
CTYPE1 = 'DEC--ZPN' / Zenithal Polynomial Projection
CRVAL1 = 14.8289418840003 / Dec at ref pixel
CUNIT1 = 'deg' / Angles are degrees always
CRPIX1 = -448.2 / Pixel coordinate at ref point
CTYPE3 = 'UTC' / linear time (UTC)
CRVAL3 = 2375.341 / Relative time of first frame
CUNIT3 = 's' / Time unit
CRPIX3 = 1.0 / Pixel coordinate at ref point
CTYPE3A = 'TT' / alternative linear time (TT)
CRVAL3A = 2440.525 / Relative time of first frame
CUNIT3A = 's' / Time unit
CRPIX3A = 1.0 / Pixel coordinate at ref point
OBSGEO-B= -24.6157 / [deg] Tel geodetic latitude (=North)+
OBSGEO-L= -70.3976 / [deg] Tel geodetic longitude (=East)+
OBSGEO-H= 2530.0000 / [m] Tel height above reference ellipsoid
CRDER3 = 0.0819 / random error in timings from fit
CSYER3 = 0.0100 / absolute time error
PC1_1 = 0.999999971570892 / WCS transform matrix element
PC1_2 = 0.000238449608932 / WCS transform matrix element
PC2_1 = -0.000621542859395 / WCS transform matrix element
PC2_2 = 0.999999806842218 / WCS transform matrix element
CDELT1 = -9.48575432499806E-5 / Axis scale at reference point
CDELT2 = 9.48683176211164E-5 / Axis scale at reference point
CDELT3 = 13.3629 / Axis scale at reference point
PV1_1 = 1. / ZPN linear term
PV1_3 = 42. / ZPN cubic term
"""
with warnings.catch_warnings():
warnings.simplefilter("ignore", (VerifyWarning, FITSFixedWarning))
WCS_TIME_CUBE = WCS(Header.fromstring(HEADER_TIME_CUBE, sep="\n"))
def test_time_cube():
# Spectral cube with a weird axis ordering
wcs = WCS_TIME_CUBE
assert wcs.pixel_n_dim == 3
assert wcs.world_n_dim == 3
assert wcs.array_shape == (11, 2048, 2048)
assert wcs.pixel_shape == (2048, 2048, 11)
assert wcs.world_axis_physical_types == ["pos.eq.dec", "pos.eq.ra", "time"]
assert wcs.world_axis_units == ["deg", "deg", "s"]
assert wcs.pixel_axis_names == ["", "", ""]
assert wcs.world_axis_names == ["", "", ""]
assert_equal(
wcs.axis_correlation_matrix,
[[True, True, False], [True, True, False], [False, False, True]],
)
components = wcs.world_axis_object_components
assert components[0] == ("celestial", 1, "spherical.lat.degree")
assert components[1] == ("celestial", 0, "spherical.lon.degree")
assert components[2][:2] == ("time", 0)
assert callable(components[2][2])
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"] == (u.deg, u.deg)
assert wcs.world_axis_object_classes["time"][0] is Time
assert wcs.world_axis_object_classes["time"][1] == ()
assert wcs.world_axis_object_classes["time"][2] == {}
assert callable(wcs.world_axis_object_classes["time"][3])
assert_allclose(
wcs.pixel_to_world_values(-449.2, 2955.6, 0),
(14.8289418840003, 2.01824372640628, 2375.341),
)
assert_allclose(
wcs.array_index_to_world_values(0, 2955.6, -449.2),
(14.8289418840003, 2.01824372640628, 2375.341),
)
assert_allclose(
wcs.world_to_pixel_values(14.8289418840003, 2.01824372640628, 2375.341),
(-449.2, 2955.6, 0),
)
assert_equal(
wcs.world_to_array_index_values(14.8289418840003, 2.01824372640628, 2375.341),
(0, 2956, -449),
)
# High-level API
coord, time = wcs.pixel_to_world(29, 39, 44)
assert isinstance(coord, SkyCoord)
assert isinstance(coord.frame, ICRS)
assert_allclose(coord.ra.deg, 1.7323356692202325)
assert_allclose(coord.dec.deg, 14.783516054817797)
assert isinstance(time, Time)
assert_allclose(time.mjd, 54746.03429755324)
coord, time = wcs.array_index_to_world(44, 39, 29)
assert isinstance(coord, SkyCoord)
assert isinstance(coord.frame, ICRS)
assert_allclose(coord.ra.deg, 1.7323356692202325)
assert_allclose(coord.dec.deg, 14.783516054817797)
assert isinstance(time, Time)
assert_allclose(time.mjd, 54746.03429755324)
x, y, z = wcs.world_to_pixel(coord, time)
assert_allclose(x, 29.0)
assert_allclose(y, 39.0)
assert_allclose(z, 44.0)
# Order of world coordinates shouldn't matter
x, y, z = wcs.world_to_pixel(time, coord)
assert_allclose(x, 29.0)
assert_allclose(y, 39.0)
assert_allclose(z, 44.0)
i, j, k = wcs.world_to_array_index(coord, time)
assert_equal(i, 44)
assert_equal(j, 39)
assert_equal(k, 29)
# Order of world coordinates shouldn't matter
i, j, k = wcs.world_to_array_index(time, coord)
assert_equal(i, 44)
assert_equal(j, 39)
assert_equal(k, 29)
###############################################################################
# The following tests are to make sure that Time objects are constructed
# correctly for a variety of combinations of WCS keywords
###############################################################################
HEADER_TIME_1D = """
SIMPLE = T
BITPIX = -32
NAXIS = 1
NAXIS1 = 2048
TIMESYS = 'UTC'
TREFPOS = 'TOPOCENT'
MJDREF = 50002.6
CTYPE1 = 'UTC'
CRVAL1 = 5
CUNIT1 = 's'
CRPIX1 = 1.0
CDELT1 = 2
OBSGEO-L= -20
OBSGEO-B= -70
OBSGEO-H= 2530
"""
if Version(wcsver) >= Version("7.1"):
HEADER_TIME_1D += "DATEREF = '1995-10-12T14:24:00'\n"
@pytest.fixture
def header_time_1d():
return Header.fromstring(HEADER_TIME_1D, sep="\n")
def assert_time_at(header, position, jd1, jd2, scale, format):
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header)
time = wcs.pixel_to_world(position)
assert_allclose(time.jd1, jd1, rtol=1e-10)
assert_allclose(time.jd2, jd2, rtol=1e-10)
assert time.format == format
assert time.scale == scale
@pytest.mark.parametrize(
"scale", ("tai", "tcb", "tcg", "tdb", "tt", "ut1", "utc", "local")
)
def test_time_1d_values(header_time_1d, scale):
# Check that Time objects are instantiated with the correct values,
# scales, and formats.
header_time_1d["CTYPE1"] = scale.upper()
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, scale, "mjd")
def test_time_1d_values_gps(header_time_1d):
# Special treatment for GPS scale
header_time_1d["CTYPE1"] = "GPS"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + (7 + 19) / 3600 / 24, "tai", "mjd")
def test_time_1d_values_deprecated(header_time_1d):
# Deprecated (in FITS) scales
header_time_1d["CTYPE1"] = "TDT"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, "tt", "mjd")
header_time_1d["CTYPE1"] = "IAT"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, "tai", "mjd")
header_time_1d["CTYPE1"] = "GMT"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, "utc", "mjd")
header_time_1d["CTYPE1"] = "ET"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, "tt", "mjd")
def test_time_1d_values_time(header_time_1d):
header_time_1d["CTYPE1"] = "TIME"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, "utc", "mjd")
header_time_1d["TIMESYS"] = "TAI"
assert_time_at(header_time_1d, 1, 2450003, 0.1 + 7 / 3600 / 24, "tai", "mjd")
@pytest.mark.remote_data
@pytest.mark.parametrize("scale", ("tai", "tcb", "tcg", "tdb", "tt", "ut1", "utc"))
def test_time_1d_roundtrip(header_time_1d, scale):
# Check that coordinates round-trip
pixel_in = np.arange(3, 10)
header_time_1d["CTYPE1"] = scale.upper()
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header_time_1d)
# Simple test
time = wcs.pixel_to_world(pixel_in)
pixel_out = wcs.world_to_pixel(time)
assert_allclose(pixel_in, pixel_out)
# Test with an intermediate change to a different scale/format
time = wcs.pixel_to_world(pixel_in).tdb
time.format = "isot"
pixel_out = wcs.world_to_pixel(time)
assert_allclose(pixel_in, pixel_out)
def test_time_1d_high_precision(header_time_1d):
# Case where the MJDREF is split into two for high precision
del header_time_1d["MJDREF"]
header_time_1d["MJDREFI"] = 52000.0
header_time_1d["MJDREFF"] = 1e-11
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header_time_1d)
time = wcs.pixel_to_world(10)
# Here we have to use a very small rtol to really test that MJDREFF is
# taken into account
assert_allclose(time.jd1, 2452001.0, rtol=1e-12)
assert_allclose(time.jd2, -0.5 + 25 / 3600 / 24 + 1e-11, rtol=1e-13)
def test_time_1d_location_geodetic(header_time_1d):
# Make sure that the location is correctly returned (geodetic case)
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header_time_1d)
time = wcs.pixel_to_world(10)
lon, lat, alt = time.location.to_geodetic()
# FIXME: alt won't work for now because ERFA doesn't implement the IAU 1976
# ellipsoid (https://github.com/astropy/astropy/issues/9420)
assert_allclose(lon.degree, -20)
assert_allclose(lat.degree, -70)
# assert_allclose(alt.to_value(u.m), 2530.)
@pytest.fixture
def header_time_1d_no_obs():
header = Header.fromstring(HEADER_TIME_1D, sep="\n")
del header["OBSGEO-L"]
del header["OBSGEO-B"]
del header["OBSGEO-H"]
return header
def test_time_1d_location_geocentric(header_time_1d_no_obs):
# Make sure that the location is correctly returned (geocentric case)
header = header_time_1d_no_obs
header["OBSGEO-X"] = 10
header["OBSGEO-Y"] = -20
header["OBSGEO-Z"] = 30
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header)
time = wcs.pixel_to_world(10)
x, y, z = time.location.to_geocentric()
assert_allclose(x.to_value(u.m), 10)
assert_allclose(y.to_value(u.m), -20)
assert_allclose(z.to_value(u.m), 30)
def test_time_1d_location_geocenter(header_time_1d_no_obs):
header_time_1d_no_obs["TREFPOS"] = "GEOCENTER"
wcs = WCS(header_time_1d_no_obs)
time = wcs.pixel_to_world(10)
x, y, z = time.location.to_geocentric()
assert_allclose(x.to_value(u.m), 0)
assert_allclose(y.to_value(u.m), 0)
assert_allclose(z.to_value(u.m), 0)
def test_time_1d_location_missing(header_time_1d_no_obs):
# Check what happens when no location is present
wcs = WCS(header_time_1d_no_obs)
with pytest.warns(
UserWarning,
match=(
"Missing or incomplete observer location "
"information, setting location in Time to None"
),
):
time = wcs.pixel_to_world(10)
assert time.location is None
def test_time_1d_location_incomplete(header_time_1d_no_obs):
# Check what happens when location information is incomplete
header_time_1d_no_obs["OBSGEO-L"] = 10.0
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header_time_1d_no_obs)
with pytest.warns(
UserWarning,
match=(
"Missing or incomplete observer location "
"information, setting location in Time to None"
),
):
time = wcs.pixel_to_world(10)
assert time.location is None
def test_time_1d_location_unsupported(header_time_1d_no_obs):
# Check what happens when TREFPOS is unsupported
header_time_1d_no_obs["TREFPOS"] = "BARYCENTER"
wcs = WCS(header_time_1d_no_obs)
with pytest.warns(
UserWarning,
match=(
"Observation location 'barycenter' is not "
"supported, setting location in Time to None"
),
):
time = wcs.pixel_to_world(10)
assert time.location is None
def test_time_1d_unsupported_ctype(header_time_1d_no_obs):
# For cases that we don't support yet, e.g. UT(...), use Time and drop sub-scale
# Case where the MJDREF is split into two for high precision
header_time_1d_no_obs["CTYPE1"] = "UT(WWV)"
if PYTEST_LT_8_0:
ctx = nullcontext()
else:
ctx = pytest.warns(
UserWarning, match="Missing or incomplete observer location information"
)
wcs = WCS(header_time_1d_no_obs)
with (
pytest.warns(
UserWarning, match="Dropping unsupported sub-scale WWV from scale UT"
),
ctx,
):
time = wcs.pixel_to_world(10)
assert isinstance(time, Time)
###############################################################################
# Extra corner cases
###############################################################################
def test_unrecognized_unit():
# TODO: Determine whether the following behavior is desirable
wcs = WCS(naxis=1)
with pytest.warns(UnitsWarning):
wcs.wcs.cunit = ["bananas // sekonds"]
assert wcs.world_axis_units == ["bananas // sekonds"]
def test_distortion_correlations():
filename = get_pkg_data_filename("../../tests/data/sip.fits")
with pytest.warns(FITSFixedWarning):
w = WCS(filename)
assert_equal(w.axis_correlation_matrix, True)
# Changing PC to an identity matrix doesn't change anything since
# distortions are still present.
w.wcs.pc = [[1, 0], [0, 1]]
assert_equal(w.axis_correlation_matrix, True)
# Nor does changing the name of the axes to make them non-celestial
w.wcs.ctype = ["X", "Y"]
assert_equal(w.axis_correlation_matrix, True)
# However once we turn off the distortions the matrix changes
w.sip = None
assert_equal(w.axis_correlation_matrix, [[True, False], [False, True]])
# If we go back to celestial coordinates then the matrix is all True again
w.wcs.ctype = ["RA---TAN", "DEC--TAN"]
assert_equal(w.axis_correlation_matrix, True)
# Or if we change to X/Y but have a non-identity PC
w.wcs.pc = [[0.9, -0.1], [0.1, 0.9]]
w.wcs.ctype = ["X", "Y"]
assert_equal(w.axis_correlation_matrix, True)
def test_custom_ctype_to_ucd_mappings():
wcs = WCS(naxis=1)
wcs.wcs.ctype = ["SPAM"]
assert wcs.world_axis_physical_types == [None]
# Check simple behavior
with custom_ctype_to_ucd_mapping({"APPLE": "food.fruit"}):
assert wcs.world_axis_physical_types == [None]
with custom_ctype_to_ucd_mapping({"APPLE": "food.fruit", "SPAM": "food.spam"}):
assert wcs.world_axis_physical_types == ["food.spam"]
# Check nesting
with custom_ctype_to_ucd_mapping({"SPAM": "food.spam"}):
with custom_ctype_to_ucd_mapping({"APPLE": "food.fruit"}):
assert wcs.world_axis_physical_types == ["food.spam"]
with custom_ctype_to_ucd_mapping({"APPLE": "food.fruit"}):
with custom_ctype_to_ucd_mapping({"SPAM": "food.spam"}):
assert wcs.world_axis_physical_types == ["food.spam"]
# Check priority in nesting
with custom_ctype_to_ucd_mapping({"SPAM": "notfood"}):
with custom_ctype_to_ucd_mapping({"SPAM": "food.spam"}):
assert wcs.world_axis_physical_types == ["food.spam"]
with custom_ctype_to_ucd_mapping({"SPAM": "food.spam"}):
with custom_ctype_to_ucd_mapping({"SPAM": "notfood"}):
assert wcs.world_axis_physical_types == ["notfood"]
def test_caching_components_and_classes():
# Make sure that when we change the WCS object, the classes and components
# are updated (we use a cache internally, so we need to make sure the cache
# is invalidated if needed)
wcs = WCS_SIMPLE_CELESTIAL.deepcopy()
assert wcs.world_axis_object_components == [
("celestial", 0, "spherical.lon.degree"),
("celestial", 1, "spherical.lat.degree"),
]
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"] == (u.deg, u.deg)
wcs.wcs.radesys = "FK5"
frame = wcs.world_axis_object_classes["celestial"][2]["frame"]
assert isinstance(frame, FK5)
assert frame.equinox.jyear == 2000.0
wcs.wcs.equinox = 2010
frame = wcs.world_axis_object_classes["celestial"][2]["frame"]
assert isinstance(frame, FK5)
assert frame.equinox.jyear == 2010.0
def test_sub_wcsapi_attributes():
# Regression test for a bug that caused some of the WCS attributes to be
# incorrect when using WCS.sub or WCS.celestial (which is an alias for sub
# with lon/lat types).
wcs = WCS_SPECTRAL_CUBE.deepcopy()
wcs.pixel_shape = (30, 40, 50)
wcs.pixel_bounds = [(-1, 11), (-2, 18), (5, 15)]
# Use celestial shortcut
wcs_sub1 = wcs.celestial
assert wcs_sub1.pixel_n_dim == 2
assert wcs_sub1.world_n_dim == 2
assert wcs_sub1.array_shape == (50, 30)
assert wcs_sub1.pixel_shape == (30, 50)
assert wcs_sub1.pixel_bounds == [(-1, 11), (5, 15)]
assert wcs_sub1.world_axis_physical_types == [
"pos.galactic.lat",
"pos.galactic.lon",
]
assert wcs_sub1.world_axis_units == ["deg", "deg"]
assert wcs_sub1.world_axis_names == ["Latitude", "Longitude"]
# Try adding axes
wcs_sub2 = wcs.sub([0, 2, 0])
assert wcs_sub2.pixel_n_dim == 3
assert wcs_sub2.world_n_dim == 3
assert wcs_sub2.array_shape == (None, 40, None)
assert wcs_sub2.pixel_shape == (None, 40, None)
assert wcs_sub2.pixel_bounds == [None, (-2, 18), None]
assert wcs_sub2.world_axis_physical_types == [None, "em.freq", None]
assert wcs_sub2.world_axis_units == ["", "Hz", ""]
assert wcs_sub2.world_axis_names == ["", "Frequency", ""]
# Use strings
wcs_sub3 = wcs.sub(["longitude", "latitude"])
assert wcs_sub3.pixel_n_dim == 2
assert wcs_sub3.world_n_dim == 2
assert wcs_sub3.array_shape == (30, 50)
assert wcs_sub3.pixel_shape == (50, 30)
assert wcs_sub3.pixel_bounds == [(5, 15), (-1, 11)]
assert wcs_sub3.world_axis_physical_types == [
"pos.galactic.lon",
"pos.galactic.lat",
]
assert wcs_sub3.world_axis_units == ["deg", "deg"]
assert wcs_sub3.world_axis_names == ["Longitude", "Latitude"]
# Now try without CNAME set
wcs.wcs.cname = [""] * wcs.wcs.naxis
wcs_sub4 = wcs.sub(["longitude", "latitude"])
assert wcs_sub4.pixel_n_dim == 2
assert wcs_sub4.world_n_dim == 2
assert wcs_sub4.array_shape == (30, 50)
assert wcs_sub4.pixel_shape == (50, 30)
assert wcs_sub4.pixel_bounds == [(5, 15), (-1, 11)]
assert wcs_sub4.world_axis_physical_types == [
"pos.galactic.lon",
"pos.galactic.lat",
]
assert wcs_sub4.world_axis_units == ["deg", "deg"]
assert wcs_sub4.world_axis_names == ["", ""]
###############################################################################
# Spectral transformations
###############################################################################
HEADER_SPECTRAL_FRAMES = """
BUNIT = 'Jy/beam'
EQUINOX = 2.000000000E+03
CTYPE1 = 'RA---SIN'
CRVAL1 = 2.60108333333E+02
CDELT1 = -2.777777845E-04
CRPIX1 = 1.0
CUNIT1 = 'deg'
CTYPE2 = 'DEC--SIN'
CRVAL2 = -9.75000000000E-01
CDELT2 = 2.777777845E-04
CRPIX2 = 1.0
CUNIT2 = 'deg'
CTYPE3 = 'FREQ'
CRVAL3 = 1.37835117405E+09
CDELT3 = 9.765625000E+04
CRPIX3 = 32.0
CUNIT3 = 'Hz'
SPECSYS = 'TOPOCENT'
RESTFRQ = 1.420405752E+09 / [Hz]
RADESYS = 'FK5'
"""
@pytest.fixture
def header_spectral_frames():
return Header.fromstring(HEADER_SPECTRAL_FRAMES, sep="\n")
def test_spectralcoord_frame(header_spectral_frames):
# This is a test to check the numerical results of transformations between
# different velocity frames. We simply make sure that the returned
# SpectralCoords are in the right frame but don't check the transformations
# since this is already done in test_spectralcoord_accuracy
# in astropy.coordinates.
with iers.conf.set_temp("auto_download", False):
obstime = Time("2009-05-04T04:44:23", scale="utc")
header = header_spectral_frames.copy()
header["MJD-OBS"] = obstime.mjd
header["CRVAL1"] = 16.33211
header["CRVAL2"] = -34.2221
header["OBSGEO-L"] = 144.2
header["OBSGEO-B"] = -20.2
header["OBSGEO-H"] = 0.0
# We start off with a WCS defined in topocentric frequency
with pytest.warns(FITSFixedWarning):
wcs_topo = WCS(header)
# We convert a single pixel coordinate to world coordinates and keep only
# the second high level object - a SpectralCoord:
sc_topo = wcs_topo.pixel_to_world(0, 0, 31)[1]
# We check that this is in topocentric frame with zero velocities
assert isinstance(sc_topo, SpectralCoord)
assert isinstance(sc_topo.observer, ITRS)
assert sc_topo.observer.obstime.isot == obstime.isot
assert_equal(sc_topo.observer.data.differentials["s"].d_xyz.value, 0)
observatory = (
EarthLocation.from_geodetic(144.2, -20.2)
.get_itrs(obstime=obstime)
.transform_to(ICRS())
)
assert (
observatory.separation_3d(sc_topo.observer.transform_to(ICRS())) < 1 * u.km
)
for specsys, expected_frame in VELOCITY_FRAMES.items():
header["SPECSYS"] = specsys
with pytest.warns(FITSFixedWarning):
wcs = WCS(header)
sc = wcs.pixel_to_world(0, 0, 31)[1]
# Now transform to the expected velocity frame, which should leave
# the spectral coordinate unchanged
sc_check = sc.with_observer_stationary_relative_to(expected_frame)
assert_quantity_allclose(sc.quantity, sc_check.quantity)
@pytest.mark.parametrize(
("ctype3", "observer"),
product(["ZOPT", "BETA", "VELO", "VRAD", "VOPT"], [False, True]),
)
def test_different_ctypes(header_spectral_frames, ctype3, observer):
header = header_spectral_frames.copy()
header["CTYPE3"] = ctype3
header["CRVAL3"] = 0.1
header["CDELT3"] = 0.001
if ctype3[0] == "V":
header["CUNIT3"] = "m s-1"
else:
header["CUNIT3"] = ""
header["RESTWAV"] = 1.420405752e09
header["MJD-OBS"] = 55197
if observer:
header["OBSGEO-L"] = 144.2
header["OBSGEO-B"] = -20.2
header["OBSGEO-H"] = 0.0
header["SPECSYS"] = "BARYCENT"
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header)
skycoord, spectralcoord = wcs.pixel_to_world(0, 0, 31)
assert isinstance(spectralcoord, SpectralCoord)
if observer:
pix = wcs.world_to_pixel(skycoord, spectralcoord)
else:
with pytest.warns(AstropyUserWarning, match="No observer defined on WCS"):
pix = wcs.world_to_pixel(skycoord, spectralcoord)
assert_allclose(pix, [0, 0, 31], rtol=1e-6, atol=1e-9)
def test_non_convergence_warning():
"""Test case for issue #11446
Since we can't define a target accuracy when plotting a WCS `all_world2pix`
should not error but only warn when the default accuracy can't be reached.
"""
# define a minimal WCS where convergence fails for certain image positions
wcs = WCS(naxis=2)
crpix = [0, 0]
a = b = ap = bp = np.zeros((4, 4))
a[3, 0] = -1.20116753e-07
test_pos_x = [1000, 1]
test_pos_y = [0, 2]
wcs.sip = Sip(a, b, ap, bp, crpix)
# first make sure the WCS works when using a low accuracy
expected = wcs.all_world2pix(test_pos_x, test_pos_y, 0, tolerance=1e-3)
# then check that it fails when using the default accuracy
with pytest.raises(NoConvergence):
wcs.all_world2pix(test_pos_x, test_pos_y, 0)
# at last check that world_to_pixel_values raises a warning but returns
# the same 'low accuray' result
with pytest.warns(UserWarning):
assert_allclose(wcs.world_to_pixel_values(test_pos_x, test_pos_y), expected)
HEADER_SPECTRAL_1D = """
CTYPE1 = 'FREQ'
CRVAL1 = 1.37835117405E+09
CDELT1 = 9.765625000E+04
CRPIX1 = 32.0
CUNIT1 = 'Hz'
SPECSYS = 'TOPOCENT'
RESTFRQ = 1.420405752E+09 / [Hz]
RADESYS = 'FK5'
"""
@pytest.fixture
def header_spectral_1d():
return Header.fromstring(HEADER_SPECTRAL_1D, sep="\n")
@pytest.mark.parametrize(
("ctype1", "observer"),
product(["ZOPT", "BETA", "VELO", "VRAD", "VOPT"], [False, True]),
)
def test_spectral_1d(header_spectral_1d, ctype1, observer):
# This is a regression test for issues that happened with 1-d WCS
# where the target is not defined but observer is.
header = header_spectral_1d.copy()
header["CTYPE1"] = ctype1
header["CRVAL1"] = 0.1
header["CDELT1"] = 0.001
if ctype1[0] == "V":
header["CUNIT1"] = "m s-1"
else:
header["CUNIT1"] = ""
header["RESTWAV"] = 1.420405752e09
header["MJD-OBS"] = 55197
if observer:
header["OBSGEO-L"] = 144.2
header["OBSGEO-B"] = -20.2
header["OBSGEO-H"] = 0.0
header["SPECSYS"] = "BARYCENT"
with warnings.catch_warnings():
warnings.simplefilter("ignore", FITSFixedWarning)
wcs = WCS(header)
# First ensure that transformations round-trip
spectralcoord = wcs.pixel_to_world(31)
assert isinstance(spectralcoord, SpectralCoord)
assert spectralcoord.target is None
assert (spectralcoord.observer is not None) is observer
if observer:
expected_message = "No target defined on SpectralCoord"
else:
expected_message = "No observer defined on WCS"
with pytest.warns(AstropyUserWarning, match=expected_message):
pix = wcs.world_to_pixel(spectralcoord)
assert_allclose(pix, [31], rtol=1e-6)
# Also make sure that we can convert a SpectralCoord on which the observer
# is not defined but the target is.
with pytest.warns(AstropyUserWarning, match="No velocity defined on frame"):
spectralcoord_no_obs = SpectralCoord(
spectralcoord.quantity,
doppler_rest=spectralcoord.doppler_rest,
doppler_convention=spectralcoord.doppler_convention,
target=ICRS(10 * u.deg, 20 * u.deg, distance=1 * u.kpc),
)
if observer:
expected_message = "No observer defined on SpectralCoord"
else:
expected_message = "No observer defined on WCS"
with pytest.warns(AstropyUserWarning, match=expected_message):
pix2 = wcs.world_to_pixel(spectralcoord_no_obs)
assert_allclose(pix2, [31], rtol=1e-6)
# And finally check case when both observer and target are defined on the
# SpectralCoord
with pytest.warns(AstropyUserWarning, match="No velocity defined on frame"):
spectralcoord_no_obs = SpectralCoord(
spectralcoord.quantity,
doppler_rest=spectralcoord.doppler_rest,
doppler_convention=spectralcoord.doppler_convention,
observer=ICRS(10 * u.deg, 20 * u.deg, distance=0 * u.kpc),
target=ICRS(10 * u.deg, 20 * u.deg, distance=1 * u.kpc),
)
if observer:
pix3 = wcs.world_to_pixel(spectralcoord_no_obs)
else:
with pytest.warns(AstropyUserWarning, match="No observer defined on WCS"):
pix3 = wcs.world_to_pixel(spectralcoord_no_obs)
assert_allclose(pix3, [31], rtol=1e-6)
HEADER_SPECTRAL_WITH_TIME = """
WCSAXES = 3
CTYPE1 = 'RA---TAN'
CTYPE2 = 'DEC--TAN'
CTYPE3 = 'WAVE'
CRVAL1 = 98.83153
CRVAL2 = -66.818
CRVAL3 = 6.4205
CRPIX1 = 21.
CRPIX2 = 22.
CRPIX3 = 1.
CDELT1 = 3.6111E-05
CDELT2 = 3.6111E-05
CDELT3 = 0.001
CUNIT1 = 'deg'
CUNIT2 = 'deg'
CUNIT3 = 'um'
MJD-AVG = 59045.41466
RADESYS = 'ICRS'
SPECSYS = 'BARYCENT'
TIMESYS = 'UTC'
"""
@pytest.fixture
def header_spectral_with_time():
return Header.fromstring(HEADER_SPECTRAL_WITH_TIME, sep="\n")
def test_spectral_with_time_kw(header_spectral_with_time):
def check_wcs(header):
assert_allclose(w.all_pix2world(*w.wcs.crpix, 1), w.wcs.crval)
sky, spec = w.pixel_to_world(*w.wcs.crpix)
assert_allclose(
(sky.spherical.lon.degree, sky.spherical.lat.degree, spec.value),
w.wcs.crval,
rtol=1e-3,
)
# Check with MJD-AVG and TIMESYS
hdr = header_spectral_with_time.copy()
with warnings.catch_warnings():
warnings.simplefilter("ignore", (VerifyWarning, FITSFixedWarning))
w = WCS(hdr)
# Make sure the correct keyword is used in a test
assert ~np.isnan(w.wcs.mjdavg)
assert np.isnan(w.wcs.mjdobs)
check_wcs(w)
# Check fall back to MJD-OBS
hdr["MJD-OBS"] = hdr["MJD-AVG"]
del hdr["MJD-AVG"]
with warnings.catch_warnings():
warnings.simplefilter("ignore", (VerifyWarning, FITSFixedWarning))
w = WCS(hdr)
# Make sure the correct keyword is used in a test
assert ~np.isnan(w.wcs.mjdobs)
assert np.isnan(w.wcs.mjdavg)
check_wcs(w)
# Check fall back to DATE--OBS
hdr["DATE-OBS"] = "2020-07-15"
del hdr["MJD-OBS"]
with warnings.catch_warnings():
warnings.simplefilter("ignore", (VerifyWarning, FITSFixedWarning))
w = WCS(hdr)
w.wcs.mjdobs = np.nan
# Make sure the correct keyword is used in a test
assert np.isnan(w.wcs.mjdobs)
assert np.isnan(w.wcs.mjdavg)
assert w.wcs.dateobs != ""
check_wcs(hdr)
# Check fall back to scale='utc'
del hdr["TIMESYS"]
check_wcs(hdr)
def test_fits_tab_time_and_units():
"""
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(get_pkg_data_filename("data/example_4d_tab.fits")) as hdul,
pytest.warns(FITSFixedWarning),
):
w = WCS(header=hdul[0].header, fobj=hdul)
with warnings.catch_warnings():
warnings.filterwarnings("ignore", message=r".*dubious year \(Note \d\)")
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 u.allclose(world[0].l, 0.06475506 * u.deg)
assert u.allclose(world[0].b, -0.02430561 * u.deg)
assert isinstance(world[1], SpectralCoord)
assert u.allclose(world[1], 24.96 * u.Hz)
assert isinstance(world[2], Time)
assert world[2].scale == "utc"
assert u.allclose(world[2].mjd, 0.00032986111111110716)
################################################################################
# Tests with Stokes
################################################################################
HEADER_POLARIZED = """
CTYPE1 = 'HPLT-TAN'
CTYPE2 = 'HPLN-TAN'
CTYPE3 = 'STOKES'
"""
@pytest.fixture
def header_polarized():
return Header.fromstring(HEADER_POLARIZED, sep="\n")
@pytest.fixture()
def wcs_polarized(header_polarized):
return WCS(header_polarized)
def test_phys_type_polarization(wcs_polarized):
w = wcs_polarized
assert w.world_axis_physical_types[2] == "phys.polarization.stokes"
def test_pixel_to_world_stokes(wcs_polarized):
w = wcs_polarized
world = w.pixel_to_world(0, 0, 0)
assert world[2] == 1
assert isinstance(world[2], StokesCoord)
assert_equal(world[2].symbol, "I")
world = w.pixel_to_world(0, 0, [0, 1, 2, 3])
assert isinstance(world[2], StokesCoord)
assert_array_equal(world[2], [1, 2, 3, 4])
assert_array_equal(world[2].symbol, ["I", "Q", "U", "V"])