1# Licensed under a 3-clause BSD style license - see LICENSE.rst 2"""Accuracy tests for GCRS coordinate transformations, primarily to/from AltAz. 3 4""" 5import os 6from importlib import metadata 7 8 9import pytest 10import numpy as np 11import erfa 12 13from astropy import units as u 14from astropy.tests.helper import assert_quantity_allclose as assert_allclose 15from astropy.time import Time 16from astropy.coordinates import ( 17 EarthLocation, get_sun, ICRS, GCRS, CIRS, ITRS, AltAz, HADec, 18 PrecessedGeocentric, CartesianRepresentation, SkyCoord, 19 CartesianDifferential, SphericalRepresentation, UnitSphericalRepresentation, 20 HCRS, HeliocentricMeanEcliptic, TEME, TETE, Angle) 21from astropy.coordinates.solar_system import _apparent_position_in_true_coordinates, get_body 22from astropy.utils import iers 23from astropy.utils.exceptions import AstropyWarning, AstropyDeprecationWarning 24from astropy.utils.compat.optional_deps import HAS_JPLEPHEM 25 26from astropy.coordinates.angle_utilities import golden_spiral_grid 27from astropy.coordinates.builtin_frames.intermediate_rotation_transforms import ( 28 get_location_gcrs, tete_to_itrs_mat, gcrs_to_cirs_mat, cirs_to_itrs_mat) 29from astropy.coordinates.builtin_frames.utils import get_jd12 30from astropy.coordinates import solar_system_ephemeris 31from astropy.units import allclose 32 33CI = os.environ.get('CI', False) == "true" 34 35 36def test_icrs_cirs(): 37 """ 38 Check a few cases of ICRS<->CIRS for consistency. 39 40 Also includes the CIRS<->CIRS transforms at different times, as those go 41 through ICRS 42 """ 43 usph = golden_spiral_grid(200) 44 dist = np.linspace(0., 1, len(usph)) * u.pc 45 inod = ICRS(usph) 46 iwd = ICRS(ra=usph.lon, dec=usph.lat, distance=dist) 47 48 cframe1 = CIRS() 49 cirsnod = inod.transform_to(cframe1) # uses the default time 50 # first do a round-tripping test 51 inod2 = cirsnod.transform_to(ICRS()) 52 assert_allclose(inod.ra, inod2.ra) 53 assert_allclose(inod.dec, inod2.dec) 54 55 # now check that a different time yields different answers 56 cframe2 = CIRS(obstime=Time('J2005')) 57 cirsnod2 = inod.transform_to(cframe2) 58 assert not allclose(cirsnod.ra, cirsnod2.ra, rtol=1e-8) 59 assert not allclose(cirsnod.dec, cirsnod2.dec, rtol=1e-8) 60 61 # parallax effects should be included, so with and w/o distance should be different 62 cirswd = iwd.transform_to(cframe1) 63 assert not allclose(cirswd.ra, cirsnod.ra, rtol=1e-8) 64 assert not allclose(cirswd.dec, cirsnod.dec, rtol=1e-8) 65 # and the distance should transform at least somehow 66 assert not allclose(cirswd.distance, iwd.distance, rtol=1e-8) 67 68 # now check that the cirs self-transform works as expected 69 cirsnod3 = cirsnod.transform_to(cframe1) # should be a no-op 70 assert_allclose(cirsnod.ra, cirsnod3.ra) 71 assert_allclose(cirsnod.dec, cirsnod3.dec) 72 73 cirsnod4 = cirsnod.transform_to(cframe2) # should be different 74 assert not allclose(cirsnod4.ra, cirsnod.ra, rtol=1e-8) 75 assert not allclose(cirsnod4.dec, cirsnod.dec, rtol=1e-8) 76 77 cirsnod5 = cirsnod4.transform_to(cframe1) # should be back to the same 78 assert_allclose(cirsnod.ra, cirsnod5.ra) 79 assert_allclose(cirsnod.dec, cirsnod5.dec) 80 81 82usph = golden_spiral_grid(200) 83dist = np.linspace(0.5, 1, len(usph)) * u.pc 84icrs_coords = [ICRS(usph), ICRS(usph.lon, usph.lat, distance=dist)] 85gcrs_frames = [GCRS(), GCRS(obstime=Time('J2005'))] 86 87 88@pytest.mark.parametrize('icoo', icrs_coords) 89def test_icrs_gcrs(icoo): 90 """ 91 Check ICRS<->GCRS for consistency 92 """ 93 gcrscoo = icoo.transform_to(gcrs_frames[0]) # uses the default time 94 # first do a round-tripping test 95 icoo2 = gcrscoo.transform_to(ICRS()) 96 assert_allclose(icoo.distance, icoo2.distance) 97 assert_allclose(icoo.ra, icoo2.ra) 98 assert_allclose(icoo.dec, icoo2.dec) 99 assert isinstance(icoo2.data, icoo.data.__class__) 100 101 # now check that a different time yields different answers 102 gcrscoo2 = icoo.transform_to(gcrs_frames[1]) 103 assert not allclose(gcrscoo.ra, gcrscoo2.ra, rtol=1e-8, atol=1e-10*u.deg) 104 assert not allclose(gcrscoo.dec, gcrscoo2.dec, rtol=1e-8, atol=1e-10*u.deg) 105 106 # now check that the cirs self-transform works as expected 107 gcrscoo3 = gcrscoo.transform_to(gcrs_frames[0]) # should be a no-op 108 assert_allclose(gcrscoo.ra, gcrscoo3.ra) 109 assert_allclose(gcrscoo.dec, gcrscoo3.dec) 110 111 gcrscoo4 = gcrscoo.transform_to(gcrs_frames[1]) # should be different 112 assert not allclose(gcrscoo4.ra, gcrscoo.ra, rtol=1e-8, atol=1e-10*u.deg) 113 assert not allclose(gcrscoo4.dec, gcrscoo.dec, rtol=1e-8, atol=1e-10*u.deg) 114 115 gcrscoo5 = gcrscoo4.transform_to(gcrs_frames[0]) # should be back to the same 116 assert_allclose(gcrscoo.ra, gcrscoo5.ra, rtol=1e-8, atol=1e-10*u.deg) 117 assert_allclose(gcrscoo.dec, gcrscoo5.dec, rtol=1e-8, atol=1e-10*u.deg) 118 119 # also make sure that a GCRS with a different geoloc/geovel gets a different answer 120 # roughly a moon-like frame 121 gframe3 = GCRS(obsgeoloc=[385000., 0, 0]*u.km, obsgeovel=[1, 0, 0]*u.km/u.s) 122 gcrscoo6 = icoo.transform_to(gframe3) # should be different 123 assert not allclose(gcrscoo.ra, gcrscoo6.ra, rtol=1e-8, atol=1e-10*u.deg) 124 assert not allclose(gcrscoo.dec, gcrscoo6.dec, rtol=1e-8, atol=1e-10*u.deg) 125 icooviag3 = gcrscoo6.transform_to(ICRS()) # and now back to the original 126 assert_allclose(icoo.ra, icooviag3.ra) 127 assert_allclose(icoo.dec, icooviag3.dec) 128 129 130@pytest.mark.parametrize('gframe', gcrs_frames) 131def test_icrs_gcrs_dist_diff(gframe): 132 """ 133 Check that with and without distance give different ICRS<->GCRS answers 134 """ 135 gcrsnod = icrs_coords[0].transform_to(gframe) 136 gcrswd = icrs_coords[1].transform_to(gframe) 137 138 # parallax effects should be included, so with and w/o distance should be different 139 assert not allclose(gcrswd.ra, gcrsnod.ra, rtol=1e-8, atol=1e-10*u.deg) 140 assert not allclose(gcrswd.dec, gcrsnod.dec, rtol=1e-8, atol=1e-10*u.deg) 141 # and the distance should transform at least somehow 142 assert not allclose(gcrswd.distance, icrs_coords[1].distance, rtol=1e-8, 143 atol=1e-10*u.pc) 144 145 146def test_cirs_to_altaz(): 147 """ 148 Check the basic CIRS<->AltAz transforms. More thorough checks implicitly 149 happen in `test_iau_fullstack` 150 """ 151 from astropy.coordinates import EarthLocation 152 153 usph = golden_spiral_grid(200) 154 dist = np.linspace(0.5, 1, len(usph)) * u.pc 155 cirs = CIRS(usph, obstime='J2000') 156 crepr = SphericalRepresentation(lon=usph.lon, lat=usph.lat, distance=dist) 157 cirscart = CIRS(crepr, obstime=cirs.obstime, representation_type=CartesianRepresentation) 158 159 loc = EarthLocation(lat=0*u.deg, lon=0*u.deg, height=0*u.m) 160 altazframe = AltAz(location=loc, obstime=Time('J2005')) 161 162 cirs2 = cirs.transform_to(altazframe).transform_to(cirs) 163 cirs3 = cirscart.transform_to(altazframe).transform_to(cirs) 164 165 # check round-tripping 166 assert_allclose(cirs.ra, cirs2.ra) 167 assert_allclose(cirs.dec, cirs2.dec) 168 assert_allclose(cirs.ra, cirs3.ra) 169 assert_allclose(cirs.dec, cirs3.dec) 170 171 172def test_cirs_to_hadec(): 173 """ 174 Check the basic CIRS<->HADec transforms. 175 """ 176 from astropy.coordinates import EarthLocation 177 178 usph = golden_spiral_grid(200) 179 dist = np.linspace(0.5, 1, len(usph)) * u.pc 180 cirs = CIRS(usph, obstime='J2000') 181 crepr = SphericalRepresentation(lon=usph.lon, lat=usph.lat, distance=dist) 182 cirscart = CIRS(crepr, obstime=cirs.obstime, representation_type=CartesianRepresentation) 183 184 loc = EarthLocation(lat=0*u.deg, lon=0*u.deg, height=0*u.m) 185 hadecframe = HADec(location=loc, obstime=Time('J2005')) 186 187 cirs2 = cirs.transform_to(hadecframe).transform_to(cirs) 188 cirs3 = cirscart.transform_to(hadecframe).transform_to(cirs) 189 190 # check round-tripping 191 assert_allclose(cirs.ra, cirs2.ra) 192 assert_allclose(cirs.dec, cirs2.dec) 193 assert_allclose(cirs.ra, cirs3.ra) 194 assert_allclose(cirs.dec, cirs3.dec) 195 196 197def test_gcrs_itrs(): 198 """ 199 Check basic GCRS<->ITRS transforms for round-tripping. 200 """ 201 usph = golden_spiral_grid(200) 202 gcrs = GCRS(usph, obstime='J2000') 203 gcrs6 = GCRS(usph, obstime='J2006') 204 205 gcrs2 = gcrs.transform_to(ITRS()).transform_to(gcrs) 206 gcrs6_2 = gcrs6.transform_to(ITRS()).transform_to(gcrs) 207 208 assert_allclose(gcrs.ra, gcrs2.ra) 209 assert_allclose(gcrs.dec, gcrs2.dec) 210 # these should be different: 211 assert not allclose(gcrs.ra, gcrs6_2.ra, rtol=1e-8) 212 assert not allclose(gcrs.dec, gcrs6_2.dec, rtol=1e-8) 213 214 # also try with the cartesian representation 215 gcrsc = gcrs.realize_frame(gcrs.data) 216 gcrsc.representation_type = CartesianRepresentation 217 gcrsc2 = gcrsc.transform_to(ITRS()).transform_to(gcrsc) 218 assert_allclose(gcrsc.spherical.lon, gcrsc2.ra) 219 assert_allclose(gcrsc.spherical.lat, gcrsc2.dec) 220 221 222def test_cirs_itrs(): 223 """ 224 Check basic CIRS<->ITRS transforms for round-tripping. 225 """ 226 usph = golden_spiral_grid(200) 227 cirs = CIRS(usph, obstime='J2000') 228 cirs6 = CIRS(usph, obstime='J2006') 229 230 cirs2 = cirs.transform_to(ITRS()).transform_to(cirs) 231 cirs6_2 = cirs6.transform_to(ITRS()).transform_to(cirs) # different obstime 232 233 # just check round-tripping 234 assert_allclose(cirs.ra, cirs2.ra) 235 assert_allclose(cirs.dec, cirs2.dec) 236 assert not allclose(cirs.ra, cirs6_2.ra) 237 assert not allclose(cirs.dec, cirs6_2.dec) 238 239 240def test_gcrs_cirs(): 241 """ 242 Check GCRS<->CIRS transforms for round-tripping. More complicated than the 243 above two because it's multi-hop 244 """ 245 usph = golden_spiral_grid(200) 246 gcrs = GCRS(usph, obstime='J2000') 247 gcrs6 = GCRS(usph, obstime='J2006') 248 249 gcrs2 = gcrs.transform_to(CIRS()).transform_to(gcrs) 250 gcrs6_2 = gcrs6.transform_to(CIRS()).transform_to(gcrs) 251 252 assert_allclose(gcrs.ra, gcrs2.ra) 253 assert_allclose(gcrs.dec, gcrs2.dec) 254 # these should be different: 255 assert not allclose(gcrs.ra, gcrs6_2.ra, rtol=1e-8) 256 assert not allclose(gcrs.dec, gcrs6_2.dec, rtol=1e-8) 257 258 # now try explicit intermediate pathways and ensure they're all consistent 259 gcrs3 = gcrs.transform_to(ITRS()).transform_to(CIRS()).transform_to(ITRS()).transform_to(gcrs) 260 assert_allclose(gcrs.ra, gcrs3.ra) 261 assert_allclose(gcrs.dec, gcrs3.dec) 262 263 gcrs4 = gcrs.transform_to(ICRS()).transform_to(CIRS()).transform_to(ICRS()).transform_to(gcrs) 264 assert_allclose(gcrs.ra, gcrs4.ra) 265 assert_allclose(gcrs.dec, gcrs4.dec) 266 267 268def test_gcrs_altaz(): 269 """ 270 Check GCRS<->AltAz transforms for round-tripping. Has multiple paths 271 """ 272 from astropy.coordinates import EarthLocation 273 274 usph = golden_spiral_grid(128) 275 gcrs = GCRS(usph, obstime='J2000')[None] # broadcast with times below 276 277 # check array times sure N-d arrays work 278 times = Time(np.linspace(2456293.25, 2456657.25, 51) * u.day, 279 format='jd')[:, None] 280 281 loc = EarthLocation(lon=10 * u.deg, lat=80. * u.deg) 282 aaframe = AltAz(obstime=times, location=loc) 283 284 aa1 = gcrs.transform_to(aaframe) 285 aa2 = gcrs.transform_to(ICRS()).transform_to(CIRS()).transform_to(aaframe) 286 aa3 = gcrs.transform_to(ITRS()).transform_to(CIRS()).transform_to(aaframe) 287 288 # make sure they're all consistent 289 assert_allclose(aa1.alt, aa2.alt) 290 assert_allclose(aa1.az, aa2.az) 291 assert_allclose(aa1.alt, aa3.alt) 292 assert_allclose(aa1.az, aa3.az) 293 294 295def test_gcrs_hadec(): 296 """ 297 Check GCRS<->HADec transforms for round-tripping. Has multiple paths 298 """ 299 from astropy.coordinates import EarthLocation 300 301 usph = golden_spiral_grid(128) 302 gcrs = GCRS(usph, obstime='J2000') # broadcast with times below 303 304 # check array times sure N-d arrays work 305 times = Time(np.linspace(2456293.25, 2456657.25, 51) * u.day, 306 format='jd')[:, np.newaxis] 307 308 loc = EarthLocation(lon=10 * u.deg, lat=80. * u.deg) 309 hdframe = HADec(obstime=times, location=loc) 310 311 hd1 = gcrs.transform_to(hdframe) 312 hd2 = gcrs.transform_to(ICRS()).transform_to(CIRS()).transform_to(hdframe) 313 hd3 = gcrs.transform_to(ITRS()).transform_to(CIRS()).transform_to(hdframe) 314 315 # make sure they're all consistent 316 assert_allclose(hd1.dec, hd2.dec) 317 assert_allclose(hd1.ha, hd2.ha) 318 assert_allclose(hd1.dec, hd3.dec) 319 assert_allclose(hd1.ha, hd3.ha) 320 321 322def test_precessed_geocentric(): 323 assert PrecessedGeocentric().equinox.jd == Time('J2000').jd 324 325 gcrs_coo = GCRS(180*u.deg, 2*u.deg, distance=10000*u.km) 326 pgeo_coo = gcrs_coo.transform_to(PrecessedGeocentric()) 327 assert np.abs(gcrs_coo.ra - pgeo_coo.ra) > 10*u.marcsec 328 assert np.abs(gcrs_coo.dec - pgeo_coo.dec) > 10*u.marcsec 329 assert_allclose(gcrs_coo.distance, pgeo_coo.distance) 330 331 gcrs_roundtrip = pgeo_coo.transform_to(GCRS()) 332 assert_allclose(gcrs_coo.ra, gcrs_roundtrip.ra) 333 assert_allclose(gcrs_coo.dec, gcrs_roundtrip.dec) 334 assert_allclose(gcrs_coo.distance, gcrs_roundtrip.distance) 335 336 pgeo_coo2 = gcrs_coo.transform_to(PrecessedGeocentric(equinox='B1850')) 337 assert np.abs(gcrs_coo.ra - pgeo_coo2.ra) > 1.5*u.deg 338 assert np.abs(gcrs_coo.dec - pgeo_coo2.dec) > 0.5*u.deg 339 assert_allclose(gcrs_coo.distance, pgeo_coo2.distance) 340 341 gcrs2_roundtrip = pgeo_coo2.transform_to(GCRS()) 342 assert_allclose(gcrs_coo.ra, gcrs2_roundtrip.ra) 343 assert_allclose(gcrs_coo.dec, gcrs2_roundtrip.dec) 344 assert_allclose(gcrs_coo.distance, gcrs2_roundtrip.distance) 345 346 347def test_precessed_geocentric_different_obstime(): 348 # Create two PrecessedGeocentric frames with different obstime 349 precessedgeo1 = PrecessedGeocentric(obstime='2021-09-07') 350 precessedgeo2 = PrecessedGeocentric(obstime='2021-06-07') 351 352 # GCRS->PrecessedGeocentric should give different results for the two frames 353 gcrs_coord = GCRS(10*u.deg, 20*u.deg, 3*u.AU, obstime=precessedgeo1.obstime) 354 pg_coord1 = gcrs_coord.transform_to(precessedgeo1) 355 pg_coord2 = gcrs_coord.transform_to(precessedgeo2) 356 assert not pg_coord1.is_equivalent_frame(pg_coord2) 357 assert not allclose(pg_coord1.cartesian.xyz, pg_coord2.cartesian.xyz) 358 359 # Looping back to GCRS should return the original coordinate 360 loopback1 = pg_coord1.transform_to(gcrs_coord) 361 loopback2 = pg_coord2.transform_to(gcrs_coord) 362 assert loopback1.is_equivalent_frame(gcrs_coord) 363 assert loopback2.is_equivalent_frame(gcrs_coord) 364 assert_allclose(loopback1.cartesian.xyz, gcrs_coord.cartesian.xyz) 365 assert_allclose(loopback2.cartesian.xyz, gcrs_coord.cartesian.xyz) 366 367 368# shared by parametrized tests below. Some use the whole AltAz, others use just obstime 369totest_frames = [AltAz(location=EarthLocation(-90*u.deg, 65*u.deg), 370 obstime=Time('J2000')), # J2000 is often a default so this might work when others don't 371 AltAz(location=EarthLocation(120*u.deg, -35*u.deg), 372 obstime=Time('J2000')), 373 AltAz(location=EarthLocation(-90*u.deg, 65*u.deg), 374 obstime=Time('2014-01-01 00:00:00')), 375 AltAz(location=EarthLocation(-90*u.deg, 65*u.deg), 376 obstime=Time('2014-08-01 08:00:00')), 377 AltAz(location=EarthLocation(120*u.deg, -35*u.deg), 378 obstime=Time('2014-01-01 00:00:00')) 379 ] 380MOONDIST = 385000*u.km # approximate moon semi-major orbit axis of moon 381MOONDIST_CART = CartesianRepresentation(3**-0.5*MOONDIST, 3**-0.5*MOONDIST, 3**-0.5*MOONDIST) 382EARTHECC = 0.017 + 0.005 # roughly earth orbital eccentricity, but with an added tolerance 383 384 385@pytest.mark.parametrize('testframe', totest_frames) 386def test_gcrs_altaz_sunish(testframe): 387 """ 388 Sanity-check that the sun is at a reasonable distance from any altaz 389 """ 390 sun = get_sun(testframe.obstime) 391 392 assert sun.frame.name == 'gcrs' 393 394 # the .to(u.au) is not necessary, it just makes the asserts on failure more readable 395 assert (EARTHECC - 1)*u.au < sun.distance.to(u.au) < (EARTHECC + 1)*u.au 396 397 sunaa = sun.transform_to(testframe) 398 assert (EARTHECC - 1)*u.au < sunaa.distance.to(u.au) < (EARTHECC + 1)*u.au 399 400 401@pytest.mark.parametrize('testframe', totest_frames) 402def test_gcrs_altaz_moonish(testframe): 403 """ 404 Sanity-check that an object resembling the moon goes to the right place with 405 a GCRS->AltAz transformation 406 """ 407 moon = GCRS(MOONDIST_CART, obstime=testframe.obstime) 408 409 moonaa = moon.transform_to(testframe) 410 411 # now check that the distance change is similar to earth radius 412 assert 1000*u.km < np.abs(moonaa.distance - moon.distance).to(u.au) < 7000*u.km 413 414 # now check that it round-trips 415 moon2 = moonaa.transform_to(moon) 416 assert_allclose(moon.cartesian.xyz, moon2.cartesian.xyz) 417 418 # also should add checks that the alt/az are different for different earth locations 419 420 421@pytest.mark.parametrize('testframe', totest_frames) 422def test_gcrs_altaz_bothroutes(testframe): 423 """ 424 Repeat of both the moonish and sunish tests above to make sure the two 425 routes through the coordinate graph are consistent with each other 426 """ 427 sun = get_sun(testframe.obstime) 428 sunaa_viaicrs = sun.transform_to(ICRS()).transform_to(testframe) 429 sunaa_viaitrs = sun.transform_to(ITRS(obstime=testframe.obstime)).transform_to(testframe) 430 431 moon = GCRS(MOONDIST_CART, obstime=testframe.obstime) 432 moonaa_viaicrs = moon.transform_to(ICRS()).transform_to(testframe) 433 moonaa_viaitrs = moon.transform_to(ITRS(obstime=testframe.obstime)).transform_to(testframe) 434 435 assert_allclose(sunaa_viaicrs.cartesian.xyz, sunaa_viaitrs.cartesian.xyz) 436 assert_allclose(moonaa_viaicrs.cartesian.xyz, moonaa_viaitrs.cartesian.xyz) 437 438 439@pytest.mark.parametrize('testframe', totest_frames) 440def test_cirs_altaz_moonish(testframe): 441 """ 442 Sanity-check that an object resembling the moon goes to the right place with 443 a CIRS<->AltAz transformation 444 """ 445 moon = CIRS(MOONDIST_CART, obstime=testframe.obstime) 446 447 moonaa = moon.transform_to(testframe) 448 assert 1000*u.km < np.abs(moonaa.distance - moon.distance).to(u.km) < 7000*u.km 449 450 # now check that it round-trips 451 moon2 = moonaa.transform_to(moon) 452 assert_allclose(moon.cartesian.xyz, moon2.cartesian.xyz) 453 454 455@pytest.mark.parametrize('testframe', totest_frames) 456def test_cirs_altaz_nodist(testframe): 457 """ 458 Check that a UnitSphericalRepresentation coordinate round-trips for the 459 CIRS<->AltAz transformation. 460 """ 461 coo0 = CIRS(UnitSphericalRepresentation(10*u.deg, 20*u.deg), obstime=testframe.obstime) 462 463 # check that it round-trips 464 coo1 = coo0.transform_to(testframe).transform_to(coo0) 465 assert_allclose(coo0.cartesian.xyz, coo1.cartesian.xyz) 466 467 468@pytest.mark.parametrize('testframe', totest_frames) 469def test_cirs_icrs_moonish(testframe): 470 """ 471 check that something like the moon goes to about the right distance from the 472 ICRS origin when starting from CIRS 473 """ 474 moonish = CIRS(MOONDIST_CART, obstime=testframe.obstime) 475 moonicrs = moonish.transform_to(ICRS()) 476 477 assert 0.97*u.au < moonicrs.distance < 1.03*u.au 478 479 480@pytest.mark.parametrize('testframe', totest_frames) 481def test_gcrs_icrs_moonish(testframe): 482 """ 483 check that something like the moon goes to about the right distance from the 484 ICRS origin when starting from GCRS 485 """ 486 moonish = GCRS(MOONDIST_CART, obstime=testframe.obstime) 487 moonicrs = moonish.transform_to(ICRS()) 488 489 assert 0.97*u.au < moonicrs.distance < 1.03*u.au 490 491 492@pytest.mark.parametrize('testframe', totest_frames) 493def test_icrs_gcrscirs_sunish(testframe): 494 """ 495 check that the ICRS barycenter goes to about the right distance from various 496 ~geocentric frames (other than testframe) 497 """ 498 # slight offset to avoid divide-by-zero errors 499 icrs = ICRS(0*u.deg, 0*u.deg, distance=10*u.km) 500 501 gcrs = icrs.transform_to(GCRS(obstime=testframe.obstime)) 502 assert (EARTHECC - 1)*u.au < gcrs.distance.to(u.au) < (EARTHECC + 1)*u.au 503 504 cirs = icrs.transform_to(CIRS(obstime=testframe.obstime)) 505 assert (EARTHECC - 1)*u.au < cirs.distance.to(u.au) < (EARTHECC + 1)*u.au 506 507 itrs = icrs.transform_to(ITRS(obstime=testframe.obstime)) 508 assert (EARTHECC - 1)*u.au < itrs.spherical.distance.to(u.au) < (EARTHECC + 1)*u.au 509 510 511@pytest.mark.parametrize('testframe', totest_frames) 512def test_icrs_altaz_moonish(testframe): 513 """ 514 Check that something expressed in *ICRS* as being moon-like goes to the 515 right AltAz distance 516 """ 517 # we use epv00 instead of get_sun because get_sun includes aberration 518 earth_pv_helio, earth_pv_bary = erfa.epv00(*get_jd12(testframe.obstime, 'tdb')) 519 earth_icrs_xyz = earth_pv_bary[0]*u.au 520 moonoffset = [0, 0, MOONDIST.value]*MOONDIST.unit 521 moonish_icrs = ICRS(CartesianRepresentation(earth_icrs_xyz + moonoffset)) 522 moonaa = moonish_icrs.transform_to(testframe) 523 524 # now check that the distance change is similar to earth radius 525 assert 1000*u.km < np.abs(moonaa.distance - MOONDIST).to(u.au) < 7000*u.km 526 527 528def test_gcrs_self_transform_closeby(): 529 """ 530 Tests GCRS self transform for objects which are nearby and thus 531 have reasonable parallax. 532 533 Moon positions were originally created using JPL DE432s ephemeris. 534 535 The two lunar positions (one geocentric, one at a defined location) 536 are created via a transformation from ICRS to two different GCRS frames. 537 538 We test that the GCRS-GCRS self transform can correctly map one GCRS 539 frame onto the other. 540 """ 541 t = Time("2014-12-25T07:00") 542 moon_geocentric = SkyCoord(GCRS(318.10579159*u.deg, 543 -11.65281165*u.deg, 544 365042.64880308*u.km, obstime=t)) 545 546 # this is the location of the Moon as seen from La Palma 547 obsgeoloc = [-5592982.59658935, -63054.1948592, 3059763.90102216]*u.m 548 obsgeovel = [4.59798494, -407.84677071, 0.]*u.m/u.s 549 moon_lapalma = SkyCoord(GCRS(318.7048445*u.deg, 550 -11.98761996*u.deg, 551 369722.8231031*u.km, 552 obstime=t, 553 obsgeoloc=obsgeoloc, 554 obsgeovel=obsgeovel)) 555 556 transformed = moon_geocentric.transform_to(moon_lapalma.frame) 557 delta = transformed.separation_3d(moon_lapalma) 558 assert_allclose(delta, 0.0*u.m, atol=1*u.m) 559 560 561def test_teme_itrf(): 562 """ 563 Test case transform from TEME to ITRF. 564 565 Test case derives from example on appendix C of Vallado, Crawford, Hujsak & Kelso (2006). 566 See https://celestrak.com/publications/AIAA/2006-6753/AIAA-2006-6753-Rev2.pdf 567 """ 568 v_itrf = CartesianDifferential(-3.225636520, -2.872451450, 5.531924446, 569 unit=u.km/u.s) 570 p_itrf = CartesianRepresentation(-1033.479383, 7901.2952740, 6380.35659580, 571 unit=u.km, differentials={'s': v_itrf}) 572 t = Time("2004-04-06T07:51:28.386") 573 574 teme = ITRS(p_itrf, obstime=t).transform_to(TEME(obstime=t)) 575 v_teme = CartesianDifferential(-4.746131487, 0.785818041, 5.531931288, 576 unit=u.km/u.s) 577 p_teme = CartesianRepresentation(5094.18016210, 6127.64465050, 6380.34453270, 578 unit=u.km, differentials={'s': v_teme}) 579 580 assert_allclose(teme.cartesian.without_differentials().xyz, 581 p_teme.without_differentials().xyz, atol=30*u.cm) 582 583 assert_allclose(teme.cartesian.differentials['s'].d_xyz, 584 p_teme.differentials['s'].d_xyz, atol=1.0*u.cm/u.s) 585 586 # test round trip 587 itrf = teme.transform_to(ITRS(obstime=t)) 588 assert_allclose( 589 itrf.cartesian.without_differentials().xyz, 590 p_itrf.without_differentials().xyz, 591 atol=100*u.cm 592 ) 593 assert_allclose( 594 itrf.cartesian.differentials['s'].d_xyz, 595 p_itrf.differentials['s'].d_xyz, 596 atol=1*u.cm/u.s 597 ) 598 599 600def test_precessedgeocentric_loopback(): 601 from_coo = PrecessedGeocentric(1*u.deg, 2*u.deg, 3*u.AU, 602 obstime='2001-01-01', equinox='2001-01-01') 603 604 # Change just the obstime 605 to_frame = PrecessedGeocentric(obstime='2001-06-30', equinox='2001-01-01') 606 607 explicit_coo = from_coo.transform_to(ICRS()).transform_to(to_frame) 608 implicit_coo = from_coo.transform_to(to_frame) 609 610 # Confirm that the explicit transformation changes the coordinate 611 assert not allclose(explicit_coo.ra, from_coo.ra, rtol=1e-10) 612 assert not allclose(explicit_coo.dec, from_coo.dec, rtol=1e-10) 613 assert not allclose(explicit_coo.distance, from_coo.distance, rtol=1e-10) 614 615 # Confirm that the loopback matches the explicit transformation 616 assert_allclose(explicit_coo.ra, implicit_coo.ra, rtol=1e-10) 617 assert_allclose(explicit_coo.dec, implicit_coo.dec, rtol=1e-10) 618 assert_allclose(explicit_coo.distance, implicit_coo.distance, rtol=1e-10) 619 620 # Change just the equinox 621 to_frame = PrecessedGeocentric(obstime='2001-01-01', equinox='2001-06-30') 622 623 explicit_coo = from_coo.transform_to(ICRS()).transform_to(to_frame) 624 implicit_coo = from_coo.transform_to(to_frame) 625 626 # Confirm that the explicit transformation changes the direction but not the distance 627 assert not allclose(explicit_coo.ra, from_coo.ra, rtol=1e-10) 628 assert not allclose(explicit_coo.dec, from_coo.dec, rtol=1e-10) 629 assert allclose(explicit_coo.distance, from_coo.distance, rtol=1e-10) 630 631 # Confirm that the loopback matches the explicit transformation 632 assert_allclose(explicit_coo.ra, implicit_coo.ra, rtol=1e-10) 633 assert_allclose(explicit_coo.dec, implicit_coo.dec, rtol=1e-10) 634 assert_allclose(explicit_coo.distance, implicit_coo.distance, rtol=1e-10) 635 636 637def test_teme_loopback(): 638 from_coo = TEME(1*u.AU, 2*u.AU, 3*u.AU, obstime='2001-01-01') 639 to_frame = TEME(obstime='2001-06-30') 640 641 explicit_coo = from_coo.transform_to(ICRS()).transform_to(to_frame) 642 implicit_coo = from_coo.transform_to(to_frame) 643 644 # Confirm that the explicit transformation changes the coordinate 645 assert not allclose(explicit_coo.cartesian.xyz, from_coo.cartesian.xyz, rtol=1e-10) 646 647 # Confirm that the loopback matches the explicit transformation 648 assert_allclose(explicit_coo.cartesian.xyz, implicit_coo.cartesian.xyz, rtol=1e-10) 649 650 651@pytest.mark.remote_data 652def test_earth_orientation_table(monkeypatch): 653 """Check that we can set the IERS table used as Earth Reference. 654 655 Use the here and now to be sure we get a difference. 656 """ 657 monkeypatch.setattr('astropy.utils.iers.conf.auto_download', True) 658 t = Time.now() 659 location = EarthLocation(lat=0*u.deg, lon=0*u.deg) 660 altaz = AltAz(location=location, obstime=t) 661 sc = SkyCoord(1*u.deg, 2*u.deg) 662 # Default: uses IERS_Auto, which will give a prediction. 663 # Note: tests run with warnings turned into errors, so it is 664 # meaningful if this passes. 665 with pytest.warns(None) as warning_lines: 666 altaz_auto = sc.transform_to(altaz) 667 668 # Server occasionally blocks IERS download in CI. 669 n_warnings = len(warning_lines) 670 if CI: 671 assert n_warnings <= 1, f'Expected at most one warning but got {n_warnings}' 672 if n_warnings == 1: 673 w_msg = str(warning_lines[0].message) 674 # This also captures unclosed socket warning that is ignored in setup.cfg 675 assert 'using local IERS-B' in w_msg or 'unclosed' in w_msg, f'Got unexpected warning: {w_msg}' 676 else: 677 assert n_warnings == 0, f'Expected no warning but got {n_warnings}' 678 679 with iers.earth_orientation_table.set(iers.IERS_B.open()): 680 with pytest.warns(AstropyWarning, match='after IERS data'): 681 altaz_b = sc.transform_to(altaz) 682 683 sep_b_auto = altaz_b.separation(altaz_auto) 684 assert_allclose(sep_b_auto, 0.0*u.deg, atol=1*u.arcsec) 685 assert sep_b_auto > 10*u.microarcsecond 686 687 # Check we returned to regular IERS system. 688 altaz_auto2 = sc.transform_to(altaz) 689 assert altaz_auto2.separation(altaz_auto) == 0. 690 691 692@pytest.mark.remote_data 693@pytest.mark.skipif(not HAS_JPLEPHEM, reason='requires jplephem') 694def test_ephemerides(): 695 """ 696 We test that using different ephemerides gives very similar results 697 for transformations 698 """ 699 t = Time("2014-12-25T07:00") 700 moon = SkyCoord(GCRS(318.10579159*u.deg, 701 -11.65281165*u.deg, 702 365042.64880308*u.km, obstime=t)) 703 704 icrs_frame = ICRS() 705 hcrs_frame = HCRS(obstime=t) 706 ecl_frame = HeliocentricMeanEcliptic(equinox=t) 707 cirs_frame = CIRS(obstime=t) 708 709 moon_icrs_builtin = moon.transform_to(icrs_frame) 710 moon_hcrs_builtin = moon.transform_to(hcrs_frame) 711 moon_helioecl_builtin = moon.transform_to(ecl_frame) 712 moon_cirs_builtin = moon.transform_to(cirs_frame) 713 714 with solar_system_ephemeris.set('jpl'): 715 moon_icrs_jpl = moon.transform_to(icrs_frame) 716 moon_hcrs_jpl = moon.transform_to(hcrs_frame) 717 moon_helioecl_jpl = moon.transform_to(ecl_frame) 718 moon_cirs_jpl = moon.transform_to(cirs_frame) 719 720 # most transformations should differ by an amount which is 721 # non-zero but of order milliarcsecs 722 sep_icrs = moon_icrs_builtin.separation(moon_icrs_jpl) 723 sep_hcrs = moon_hcrs_builtin.separation(moon_hcrs_jpl) 724 sep_helioecl = moon_helioecl_builtin.separation(moon_helioecl_jpl) 725 sep_cirs = moon_cirs_builtin.separation(moon_cirs_jpl) 726 727 assert_allclose([sep_icrs, sep_hcrs, sep_helioecl], 0.0*u.deg, atol=10*u.mas) 728 assert all(sep > 10*u.microarcsecond for sep in (sep_icrs, sep_hcrs, sep_helioecl)) 729 730 # CIRS should be the same 731 assert_allclose(sep_cirs, 0.0*u.deg, atol=1*u.microarcsecond) 732 733 734def test_tete_transforms(): 735 """ 736 We test the TETE transforms for proper behaviour here. 737 738 The TETE transforms are tested for accuracy against JPL Horizons in 739 test_solar_system.py. Here we are looking to check for consistency and 740 errors in the self transform. 741 """ 742 loc = EarthLocation.from_geodetic("-22°57'35.1", "-67°47'14.1", 5186*u.m) 743 time = Time('2020-04-06T00:00') 744 p, v = loc.get_gcrs_posvel(time) 745 746 gcrs_frame = GCRS(obstime=time, obsgeoloc=p, obsgeovel=v) 747 moon = SkyCoord(169.24113968*u.deg, 10.86086666*u.deg, 358549.25381755*u.km, frame=gcrs_frame) 748 749 tete_frame = TETE(obstime=time, location=loc) 750 # need to set obsgeoloc/vel explicitly or skycoord behaviour over-writes 751 tete_geo = TETE(obstime=time, location=EarthLocation(*([0, 0, 0]*u.km))) 752 753 # test self-transform by comparing to GCRS-TETE-ITRS-TETE route 754 tete_coo1 = moon.transform_to(tete_frame) 755 tete_coo2 = moon.transform_to(tete_geo) 756 assert_allclose(tete_coo1.separation_3d(tete_coo2), 0*u.mm, atol=1*u.mm) 757 758 # test TETE-ITRS transform by comparing GCRS-CIRS-ITRS to GCRS-TETE-ITRS 759 itrs1 = moon.transform_to(CIRS()).transform_to(ITRS()) 760 itrs2 = moon.transform_to(TETE()).transform_to(ITRS()) 761 assert_allclose(itrs1.separation_3d(itrs2), 0*u.mm, atol=1*u.mm) 762 763 # test round trip GCRS->TETE->GCRS 764 new_moon = moon.transform_to(TETE()).transform_to(moon) 765 assert_allclose(new_moon.separation_3d(moon), 0*u.mm, atol=1*u.mm) 766 767 # test round trip via ITRS 768 tete_rt = tete_coo1.transform_to(ITRS(obstime=time)).transform_to(tete_coo1) 769 assert_allclose(tete_rt.separation_3d(tete_coo1), 0*u.mm, atol=1*u.mm) 770 771 # ensure deprecated routine remains consistent 772 # make sure test raises warning! 773 with pytest.warns(AstropyDeprecationWarning, match='The use of'): 774 tete_alt = _apparent_position_in_true_coordinates(moon) 775 assert_allclose(tete_coo1.separation_3d(tete_alt), 0*u.mm, atol=100*u.mm) 776 777 778def test_straight_overhead(): 779 """ 780 With a precise CIRS<->AltAz transformation this should give Alt=90 exactly 781 782 If the CIRS self-transform breaks it won't, due to improper treatment of aberration 783 """ 784 t = Time('J2010') 785 obj = EarthLocation(-1*u.deg, 52*u.deg, height=10.*u.km) 786 home = EarthLocation(-1*u.deg, 52*u.deg, height=0.*u.km) 787 788 # An object that appears straight overhead - FOR A GEOCENTRIC OBSERVER. 789 # Note, this won't be overhead for a topocentric observer because of 790 # aberration. 791 cirs_geo = obj.get_itrs(t).transform_to(CIRS(obstime=t)) 792 793 # now get the Geocentric CIRS position of observatory 794 obsrepr = home.get_itrs(t).transform_to(CIRS(obstime=t)).cartesian 795 796 # topocentric CIRS position of a straight overhead object 797 cirs_repr = cirs_geo.cartesian - obsrepr 798 799 # create a CIRS object that appears straight overhead for a TOPOCENTRIC OBSERVER 800 topocentric_cirs_frame = CIRS(obstime=t, location=home) 801 cirs_topo = topocentric_cirs_frame.realize_frame(cirs_repr) 802 803 # Check AltAz (though Azimuth can be anything so is not tested). 804 aa = cirs_topo.transform_to(AltAz(obstime=t, location=home)) 805 assert_allclose(aa.alt, 90*u.deg, atol=1*u.uas, rtol=0) 806 807 # Check HADec. 808 hd = cirs_topo.transform_to(HADec(obstime=t, location=home)) 809 assert_allclose(hd.ha, 0*u.hourangle, atol=1*u.uas, rtol=0) 810 assert_allclose(hd.dec, 52*u.deg, atol=1*u.uas, rtol=0) 811 812 813def jplephem_ge(minversion): 814 """Check if jplephem is installed and has version >= minversion.""" 815 # This is a separate routine since somehow with pyinstaller the stanza 816 # not HAS_JPLEPHEM or metadata.version('jplephem') < '2.15' 817 # leads to a module not found error. 818 try: 819 return HAS_JPLEPHEM and metadata.version('jplephem') >= minversion 820 except Exception: 821 return False 822 823 824@pytest.mark.remote_data 825@pytest.mark.skipif(not jplephem_ge('2.15'), reason='requires jplephem >= 2.15') 826def test_aa_hd_high_precision(): 827 """These tests are provided by @mkbrewer - see issue #10356. 828 829 The code that produces them agrees very well (<0.5 mas) with SkyField once Polar motion 830 is turned off, but SkyField does not include polar motion, so a comparison to Skyfield 831 or JPL Horizons will be ~1" off. 832 833 The absence of polar motion within Skyfield and the disagreement between Skyfield and Horizons 834 make high precision comparisons to those codes difficult. 835 836 Updated 2020-11-29, after the comparison between codes became even better, 837 down to 100 nas. 838 839 NOTE: the agreement reflects consistency in approach between two codes, 840 not necessarily absolute precision. If this test starts failing, the 841 tolerance can and should be weakened *if* it is clear that the change is 842 due to an improvement (e.g., a new IAU precession model). 843 844 """ 845 lat = -22.959748*u.deg 846 lon = -67.787260*u.deg 847 elev = 5186*u.m 848 loc = EarthLocation.from_geodetic(lon, lat, elev) 849 # Note: at this level of precision for the comparison, we have to include 850 # the location in the time, as it influences the transformation to TDB. 851 t = Time('2017-04-06T00:00:00.0', location=loc) 852 with solar_system_ephemeris.set('de430'): 853 moon = get_body('moon', t, loc) 854 moon_aa = moon.transform_to(AltAz(obstime=t, location=loc)) 855 moon_hd = moon.transform_to(HADec(obstime=t, location=loc)) 856 857 # Numbers from 858 # https://github.com/astropy/astropy/pull/11073#issuecomment-735486271 859 # updated in https://github.com/astropy/astropy/issues/11683 860 TARGET_AZ, TARGET_EL = 15.032673509956*u.deg, 50.303110133923*u.deg 861 TARGET_DISTANCE = 376252883.247239*u.m 862 assert_allclose(moon_aa.az, TARGET_AZ, atol=0.1*u.uas, rtol=0) 863 assert_allclose(moon_aa.alt, TARGET_EL, atol=0.1*u.uas, rtol=0) 864 assert_allclose(moon_aa.distance, TARGET_DISTANCE, atol=0.1*u.mm, rtol=0) 865 ha, dec = erfa.ae2hd(moon_aa.az.to_value(u.radian), moon_aa.alt.to_value(u.radian), 866 lat.to_value(u.radian)) 867 ha = u.Quantity(ha, u.radian, copy=False) 868 dec = u.Quantity(dec, u.radian, copy=False) 869 assert_allclose(moon_hd.ha, ha, atol=0.1*u.uas, rtol=0) 870 assert_allclose(moon_hd.dec, dec, atol=0.1*u.uas, rtol=0) 871 872 873def test_aa_high_precision_nodata(): 874 """ 875 These tests are designed to ensure high precision alt-az transforms. 876 877 They are a slight fudge since the target values come from astropy itself. They are generated 878 with a version of the code that passes the tests above, but for the internal solar system 879 ephemerides to avoid the use of remote data. 880 """ 881 # Last updated when switching to erfa 2.0.0 and its moon98 function. 882 TARGET_AZ, TARGET_EL = 15.03231495*u.deg, 50.3027193*u.deg 883 lat = -22.959748*u.deg 884 lon = -67.787260*u.deg 885 elev = 5186*u.m 886 loc = EarthLocation.from_geodetic(lon, lat, elev) 887 t = Time('2017-04-06T00:00:00.0') 888 889 moon = get_body('moon', t, loc) 890 moon_aa = moon.transform_to(AltAz(obstime=t, location=loc)) 891 assert_allclose(moon_aa.az - TARGET_AZ, 0*u.mas, atol=0.5*u.mas) 892 assert_allclose(moon_aa.alt - TARGET_EL, 0*u.mas, atol=0.5*u.mas) 893 894 895class TestGetLocationGCRS: 896 # TETE and CIRS use get_location_gcrs to get obsgeoloc and obsgeovel 897 # with knowledge of some of the matrices. Check that this is consistent 898 # with a direct transformation. 899 def setup_class(cls): 900 cls.loc = loc = EarthLocation.from_geodetic( 901 np.linspace(0, 360, 6)*u.deg, np.linspace(-90, 90, 6)*u.deg, 100*u.m) 902 cls.obstime = obstime = Time(np.linspace(2000, 2010, 6), format='jyear') 903 # Get comparison via a full transformation. We do not use any methods 904 # of EarthLocation, since those depend on the fast transform. 905 loc_itrs = ITRS(loc.x, loc.y, loc.z, obstime=obstime) 906 zeros = np.broadcast_to(0. * (u.km / u.s), (3,) + loc_itrs.shape, subok=True) 907 loc_itrs.data.differentials['s'] = CartesianDifferential(zeros) 908 loc_gcrs_cart = loc_itrs.transform_to(GCRS(obstime=obstime)).cartesian 909 cls.obsgeoloc = loc_gcrs_cart.without_differentials() 910 cls.obsgeovel = loc_gcrs_cart.differentials['s'].to_cartesian() 911 912 def check_obsgeo(self, obsgeoloc, obsgeovel): 913 assert_allclose(obsgeoloc.xyz, self.obsgeoloc.xyz, atol=.1*u.um, rtol=0.) 914 assert_allclose(obsgeovel.xyz, self.obsgeovel.xyz, atol=.1*u.mm/u.s, rtol=0.) 915 916 def test_get_gcrs_posvel(self): 917 # Really just a sanity check 918 self.check_obsgeo(*self.loc.get_gcrs_posvel(self.obstime)) 919 920 def test_tete_quick(self): 921 # Following copied from intermediate_rotation_transforms.gcrs_to_tete 922 rbpn = erfa.pnm06a(*get_jd12(self.obstime, 'tt')) 923 loc_gcrs_frame = get_location_gcrs(self.loc, self.obstime, 924 tete_to_itrs_mat(self.obstime, rbpn=rbpn), 925 rbpn) 926 self.check_obsgeo(loc_gcrs_frame.obsgeoloc, loc_gcrs_frame.obsgeovel) 927 928 def test_cirs_quick(self): 929 cirs_frame = CIRS(location=self.loc, obstime=self.obstime) 930 # Following copied from intermediate_rotation_transforms.gcrs_to_cirs 931 pmat = gcrs_to_cirs_mat(cirs_frame.obstime) 932 loc_gcrs_frame = get_location_gcrs(self.loc, self.obstime, 933 cirs_to_itrs_mat(cirs_frame.obstime), pmat) 934 self.check_obsgeo(loc_gcrs_frame.obsgeoloc, loc_gcrs_frame.obsgeovel) 935