# -*- coding: utf-8 -*- # Licensed under a 3-clause BSD style license - see LICENSE.rst import warnings import pytest import numpy as np from numpy import testing as npt import erfa from astropy import units as u from astropy.time import Time from astropy.coordinates.builtin_frames import ICRS, AltAz from astropy.coordinates.builtin_frames.utils import get_jd12 from astropy.coordinates import EarthLocation from astropy.coordinates import SkyCoord from astropy.utils import iers from astropy.coordinates.angle_utilities import golden_spiral_grid # These fixtures are used in test_iau_fullstack @pytest.fixture(scope="function") def fullstack_icrs(): rep = golden_spiral_grid(size=1000) return ICRS(rep) @pytest.fixture(scope="function") def fullstack_fiducial_altaz(fullstack_icrs): altazframe = AltAz(location=EarthLocation(lat=0*u.deg, lon=0*u.deg, height=0*u.m), obstime=Time('J2000')) with warnings.catch_warnings(): # Ignore remote_data warning warnings.simplefilter('ignore') result = fullstack_icrs.transform_to(altazframe) return result @pytest.fixture(scope="function", params=['J2000.1', 'J2010']) def fullstack_times(request): return Time(request.param) @pytest.fixture(scope="function", params=[(0, 0, 0), (23, 0, 0), (-70, 0, 0), (0, 100, 0), (23, 0, 3000)]) def fullstack_locations(request): return EarthLocation(lat=request.param[0]*u.deg, lon=request.param[0]*u.deg, height=request.param[0]*u.m) @pytest.fixture(scope="function", params=[(0*u.bar, 0*u.deg_C, 0, 1*u.micron), (1*u.bar, 0*u.deg_C, 0*u.one, 1*u.micron), (1*u.bar, 10*u.deg_C, 0, 1*u.micron), (1*u.bar, 0*u.deg_C, 50*u.percent, 1*u.micron), (1*u.bar, 0*u.deg_C, 0, 21*u.cm)]) def fullstack_obsconditions(request): return request.param def _erfa_check(ira, idec, astrom): """ This function does the same thing the astropy layer is supposed to do, but all in erfa """ cra, cdec = erfa.atciq(ira, idec, 0, 0, 0, 0, astrom) az, zen, ha, odec, ora = erfa.atioq(cra, cdec, astrom) alt = np.pi/2-zen cra2, cdec2 = erfa.atoiq('A', az, zen, astrom) ira2, idec2 = erfa.aticq(cra2, cdec2, astrom) dct = locals() del dct['astrom'] return dct def test_iau_fullstack(fullstack_icrs, fullstack_fiducial_altaz, fullstack_times, fullstack_locations, fullstack_obsconditions): """ Test the full transform from ICRS <-> AltAz """ # create the altaz frame altazframe = AltAz(obstime=fullstack_times, location=fullstack_locations, pressure=fullstack_obsconditions[0], temperature=fullstack_obsconditions[1], relative_humidity=fullstack_obsconditions[2], obswl=fullstack_obsconditions[3]) aacoo = fullstack_icrs.transform_to(altazframe) # compare aacoo to the fiducial AltAz - should always be different assert np.all(np.abs(aacoo.alt - fullstack_fiducial_altaz.alt) > 50*u.milliarcsecond) assert np.all(np.abs(aacoo.az - fullstack_fiducial_altaz.az) > 50*u.milliarcsecond) # if the refraction correction is included, we *only* do the comparisons # where altitude >5 degrees. The SOFA guides imply that below 5 is where # where accuracy gets more problematic, and testing reveals that alt<~0 # gives garbage round-tripping, and <10 can give ~1 arcsec uncertainty if fullstack_obsconditions[0].value == 0: # but if there is no refraction correction, check everything msk = slice(None) tol = 5*u.microarcsecond else: msk = aacoo.alt > 5*u.deg # most of them aren't this bad, but some of those at low alt are offset # this much. For alt > 10, this is always better than 100 masec tol = 750*u.milliarcsecond # now make sure the full stack round-tripping works icrs2 = aacoo.transform_to(ICRS()) adras = np.abs(fullstack_icrs.ra - icrs2.ra)[msk] addecs = np.abs(fullstack_icrs.dec - icrs2.dec)[msk] assert np.all(adras < tol), f'largest RA change is {np.max(adras.arcsec * 1000)} mas, > {tol}' assert np.all(addecs < tol), f'largest Dec change is {np.max(addecs.arcsec * 1000)} mas, > {tol}' # check that we're consistent with the ERFA alt/az result iers_tab = iers.earth_orientation_table.get() xp, yp = u.Quantity(iers_tab.pm_xy(fullstack_times)).to_value(u.radian) lon = fullstack_locations.geodetic[0].to_value(u.radian) lat = fullstack_locations.geodetic[1].to_value(u.radian) height = fullstack_locations.geodetic[2].to_value(u.m) jd1, jd2 = get_jd12(fullstack_times, 'utc') pressure = fullstack_obsconditions[0].to_value(u.hPa) temperature = fullstack_obsconditions[1].to_value(u.deg_C) # Relative humidity can be a quantity or a number. relative_humidity = u.Quantity(fullstack_obsconditions[2], u.one).value obswl = fullstack_obsconditions[3].to_value(u.micron) astrom, eo = erfa.apco13(jd1, jd2, fullstack_times.delta_ut1_utc, lon, lat, height, xp, yp, pressure, temperature, relative_humidity, obswl) erfadct = _erfa_check(fullstack_icrs.ra.rad, fullstack_icrs.dec.rad, astrom) npt.assert_allclose(erfadct['alt'], aacoo.alt.radian, atol=1e-7) npt.assert_allclose(erfadct['az'], aacoo.az.radian, atol=1e-7) def test_fiducial_roudtrip(fullstack_icrs, fullstack_fiducial_altaz): """ Test the full transform from ICRS <-> AltAz """ aacoo = fullstack_icrs.transform_to(fullstack_fiducial_altaz) # make sure the round-tripping works icrs2 = aacoo.transform_to(ICRS()) npt.assert_allclose(fullstack_icrs.ra.deg, icrs2.ra.deg) npt.assert_allclose(fullstack_icrs.dec.deg, icrs2.dec.deg) def test_future_altaz(): """ While this does test the full stack, it is mostly meant to check that a warning is raised when attempting to get to AltAz in the future (beyond IERS tables) """ from astropy.utils.exceptions import AstropyWarning # this is an ugly hack to get the warning to show up even if it has already # appeared from astropy.coordinates.builtin_frames import utils if hasattr(utils, '__warningregistry__'): utils.__warningregistry__.clear() location = EarthLocation(lat=0*u.deg, lon=0*u.deg) t = Time('J2161') # check that these message(s) appear among any other warnings. If tests are run with # --remote-data then the IERS table will be an instance of IERS_Auto which is # assured of being "fresh". In this case getting times outside the range of the # table does not raise an exception. Only if using IERS_B (which happens without # --remote-data, i.e. for all CI testing) do we expect another warning. with pytest.warns(AstropyWarning, match=r"Tried to get polar motions for " "times after IERS data is valid.*") as found_warnings: SkyCoord(1*u.deg, 2*u.deg).transform_to(AltAz(location=location, obstime=t)) if isinstance(iers.earth_orientation_table.get(), iers.IERS_B): messages_found = ["(some) times are outside of range covered by IERS " "table." in str(w.message) for w in found_warnings] assert any(messages_found)