From cb3ef7a3a8715d7614869c92156c90d82132ae52 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 06:49:24 -0400 Subject: [PATCH] Include parallax correction to unit tests. --- .../tweakreg/tests/test_astrometric_utils.py | 184 +++++++++++++++--- 1 file changed, 160 insertions(+), 24 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index ad887a28f..4fa0ed56e 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -8,6 +8,7 @@ from astropy import coordinates as coord from astropy import table from astropy import units as u +from astropy.time import Time from astropy.modeling import models from astropy.modeling.models import RotationSequence3D, Scale, Shift from gwcs import coordinate_frames as cf @@ -29,6 +30,123 @@ def __init__(self, *args, **kwargs): raise requests.exceptions.ConnectionError +def get_proper_motion_correction(epoch, gaia_ref_epoch_coords, gaia_ref_epoch): + """ + Calculates the proper motion correction for a given epoch and Gaia reference epoch coordinates. + + Parameters + ---------- + epoch : float + The epoch for which the proper motion correction is calculated. + gaia_ref_epoch_coords : dict + A dictionary containing Gaia reference epoch coordinates. + gaia_ref_epoch : float + The Gaia reference epoch. + + Returns + ------- + None + + Examples + -------- + .. code-block:: python + epoch = 2022.5 + gaia_coords = { + "ra": 180.0, + "dec": 45.0, + "pmra": 2.0, + "pmdec": 1.5 + } + gaia_ref_epoch = 2020.0 + get_proper_motion_correction(epoch, gaia_coords, gaia_ref_epoch) + """ + + expected_new_dec = ( + np.array( + gaia_ref_epoch_coords["dec"] * 3600 + + (epoch - gaia_ref_epoch) * gaia_ref_epoch_coords["pmdec"] / 1000 + ) + / 3600 + ) + average_dec = np.array( + [ + np.mean([new, old]) + for new, old in zip(expected_new_dec, gaia_ref_epoch_coords["dec"]) + ] + ) + pmra = gaia_ref_epoch_coords["pmra"] / np.cos(np.deg2rad(average_dec)) + + # angular displacement components + gaia_ref_epoch_coords["pm_delta_dec"] = u.Quantity( + (epoch - gaia_ref_epoch) * gaia_ref_epoch_coords["pmdec"] / 1000, + unit=u.arcsec, + ).to(u.deg) + gaia_ref_epoch_coords["pm_delta_ra"] = u.Quantity( + (epoch - gaia_ref_epoch) * (pmra / 1000), unit=u.arcsec + ).to(u.deg) + + +def get_parallax_correction(epoch, gaia_ref_epoch_coords): + """ + Calculates the parallax correction for a given epoch and Gaia reference epoch + coordinates. + Calculations based on Chapter 8 of "Spherical Astronomy" by Robin M. Green. + + Parameters + ---------- + epoch : float + The epoch for which the parallax correction is calculated. + gaia_ref_epoch_coords : dict + A dictionary containing Gaia reference epoch coordinates. + + Returns + ------- + None + + Examples + -------- + .. code-block :: python + epoch = 2022.5 + gaia_coords = { + "ra": 180.0, + "dec": 45.0, + "parallax": 10.0 + } + get_parallax_correction(epoch, gaia_coords) + """ + + obs_date = Time(epoch, format="decimalyear") + earths_center_barycentric_coords = coord.get_body_barycentric( + "earth", obs_date + ) + earth_X = earths_center_barycentric_coords.x + earth_Y = earths_center_barycentric_coords.y + earth_Z = earths_center_barycentric_coords.z + + # angular displacement components + # (see eq. 8.15 of "Spherical Astronomy" by Robert M. Green) + gaia_ref_epoch_coords["parallax_delta_ra"] = ( + u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad) + * (1 / np.cos(gaia_ref_epoch_coords["dec"] / 180 * np.pi)) + * ( + earth_X.value * np.sin(gaia_ref_epoch_coords["ra"] / 180 * np.pi) + - earth_Y.value * np.cos(gaia_ref_epoch_coords["ra"] / 180 * np.pi) + ) + ).to("deg") + gaia_ref_epoch_coords["parallax_delta_dec"] = ( + u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad) + * ( + earth_X.value + * np.cos(gaia_ref_epoch_coords["ra"] / 180 * np.pi) + * np.sin(gaia_ref_epoch_coords["dec"] / 180 * np.pi) + + earth_Y.value + * np.sin(gaia_ref_epoch_coords["ra"] / 180 * np.pi) + * np.sin(gaia_ref_epoch_coords["dec"] / 180 * np.pi) + - earth_Z.value * np.cos(gaia_ref_epoch_coords["dec"] / 180 * np.pi) + ) + ).to("deg") + + def update_wcsinfo(input_dm): """ Update WCSInfo with realistic data (i.e. obtained from romanisim simulations). @@ -80,7 +198,9 @@ def create_wcs_for_tweakreg_pipeline(input_dm, shift_1=0, shift_2=0): tel2sky = _create_tel2sky_model(input_dm) # create required frames - detector = cf.Frame2D(name="detector", axes_order=(0, 1), unit=(u.pix, u.pix)) + detector = cf.Frame2D( + name="detector", axes_order=(0, 1), unit=(u.pix, u.pix) + ) v2v3 = cf.Frame2D( name="v2v3", axes_order=(0, 1), @@ -338,7 +458,9 @@ def test_create_astrometric_catalog_write_results_to_disk(tmp_path, base_image): ("GAIADR3", None), ], ) -def test_create_astrometric_catalog_using_epoch(tmp_path, catalog, epoch, request): +def test_create_astrometric_catalog_using_epoch( + tmp_path, catalog, epoch, request +): """Test fetching data from supported catalogs for a specific epoch.""" output_filename = "ref_cat.ecsv" img = request.getfixturevalue("base_image")(shift_1=1000, shift_2=1000) @@ -481,39 +603,53 @@ def test_get_catalog_valid_parameters_but_no_sources_returned(): def test_get_catalog_using_epoch(ra, dec, epoch): """Test that get_catalog returns coordinates corrected by proper motion.""" - result = get_catalog(ra, dec, epoch=epoch) + result_all = get_catalog(ra, dec, epoch=epoch) + + # select sources with reliable astrometric solutions based on the + # parallax_over_error parameter as discussed in Fabricius et al. 2021 + # (https://www.aanda.org/articles/aa/full_html/2021/05/aa39834-20/aa39834-20.html) + mask = result_all["parallax_over_error"] > 5 + + result = result_all[mask] returned_ra = np.array(result["ra"]) returned_dec = np.array(result["dec"]) # get GAIA data and update coords to requested epoch using pm measurements gaia_ref_epoch = 2016.0 - gaia_ref_epoch_coords = get_catalog(ra, dec, epoch=gaia_ref_epoch) + gaia_ref_epoch_coords_all = get_catalog(ra, dec, epoch=gaia_ref_epoch) - expected_new_dec = ( - np.array( - gaia_ref_epoch_coords["dec"] * 3600 - + (epoch - gaia_ref_epoch) * gaia_ref_epoch_coords["pmdec"] / 1000 - ) - / 3600 + gaia_ref_epoch_coords = gaia_ref_epoch_coords_all[mask] + + # calculate proper motion corrections + get_proper_motion_correction( + epoch=epoch, + gaia_ref_epoch_coords=gaia_ref_epoch_coords, + gaia_ref_epoch=gaia_ref_epoch, ) - average_dec = np.array( - [ - np.mean([new, old]) - for new, old in zip(expected_new_dec, gaia_ref_epoch_coords["dec"]) - ] + # calculate parallax corrections + get_parallax_correction( + epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords ) - pmra = gaia_ref_epoch_coords["pmra"] / np.cos(np.deg2rad(average_dec)) - expected_new_ra = ( - np.array( - gaia_ref_epoch_coords["ra"] * 3600 - + (epoch - gaia_ref_epoch) * (pmra / 1000) - ) - / 3600 + + # calculate the expected coordinates value after corrections have been applied to + # Gaia's reference epoch coordinates + expected_ra = ( + gaia_ref_epoch_coords["ra"] + + gaia_ref_epoch_coords["pm_delta_ra"] + + gaia_ref_epoch_coords["parallax_delta_ra"] + ) + expected_dec = ( + gaia_ref_epoch_coords["dec"] + + gaia_ref_epoch_coords["pm_delta_dec"] + + gaia_ref_epoch_coords["parallax_delta_dec"] ) assert len(result) > 0 - assert np.isclose(returned_ra, expected_new_ra, atol=1e-10, rtol=1e-9).all() - assert np.isclose(returned_dec, expected_new_dec, atol=1e-10, rtol=1e-9).all() + # choosing atol=0 corresponds to abs(a - b) / abs(b) <= rtol, + # where a is the returned value from the VO API and b is the expected + # value from applying the calculated PM and parallax corrections. + assert np.isclose(returned_ra, expected_ra, atol=0, rtol=1e-5).all() + assert np.isclose(returned_dec, expected_dec, atol=0, rtol=1e-5).all() def test_get_catalog_timeout():