From 1fede85e1877d1855dbd74e800606c7d09217fa5 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Mon, 2 Oct 2023 10:24:38 -0400 Subject: [PATCH 01/13] Fix TweakReg's catalog metadata. --- romancal/tweakreg/tweakreg_step.py | 43 ++++++++++++++++++++++-------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/romancal/tweakreg/tweakreg_step.py b/romancal/tweakreg/tweakreg_step.py index 7db401a4c..0c3b4981e 100644 --- a/romancal/tweakreg/tweakreg_step.py +++ b/romancal/tweakreg/tweakreg_step.py @@ -31,7 +31,9 @@ def _oxford_or_str_join(str_list): elif nelem == 2: return f"{str_list[0]} or {str_list[1]}" else: - return ", ".join(map(repr, str_list[:-1])) + ", or " + repr(str_list[-1]) + return ( + ", ".join(map(repr, str_list[:-1])) + ", or " + repr(str_list[-1]) + ) SINGLE_GROUP_REFCAT = ["GAIADR3", "GAIADR2", "GAIADR1"] @@ -144,7 +146,9 @@ def process(self, input): self.catalog_path = os.getcwd() self.catalog_path = Path(self.catalog_path).as_posix() - self.log.info(f"All source catalogs will be saved to: {self.catalog_path}") + self.log.info( + f"All source catalogs will be saved to: {self.catalog_path}" + ) if self.abs_refcat is None or len(self.abs_refcat.strip()) == 0: self.abs_refcat = DEFAULT_ABS_REFCAT @@ -236,7 +240,9 @@ def process(self, input): grp_img = list(images.models_grouped) self.log.info("") - self.log.info(f"Number of image groups to be aligned: {len(grp_img):d}.") + self.log.info( + f"Number of image groups to be aligned: {len(grp_img):d}." + ) self.log.info("Image groups:") if len(grp_img) == 1 and not ALIGN_TO_ABS_REFCAT: @@ -246,7 +252,9 @@ def process(self, input): self.log.info("") # we need at least two exposures to perform image alignment - self.log.warning("At least two exposures are required for image alignment.") + self.log.warning( + "At least two exposures are required for image alignment." + ) self.log.warning("Nothing to do. Skipping 'TweakRegStep'...") self.skip = True for model in images: @@ -320,7 +328,8 @@ def process(self, input): except ValueError as e: msg = e.args[0] if ( - msg == "Too few input images (or groups of images) with non-empty" + msg + == "Too few input images (or groups of images) with non-empty" " catalogs." ): # we need at least two exposures to perform image alignment @@ -328,7 +337,9 @@ def process(self, input): self.log.warning( "At least two exposures are required for image alignment." ) - self.log.warning("Nothing to do. Skipping 'TweakRegStep'...") + self.log.warning( + "Nothing to do. Skipping 'TweakRegStep'..." + ) for model in images: model.meta.cal_step["tweakreg"] = "SKIPPED" if not ALIGN_TO_ABS_REFCAT: @@ -354,7 +365,9 @@ def process(self, input): for model in images: model.meta.cal_step["tweakreg"] = "SKIPPED" if ALIGN_TO_ABS_REFCAT: - self.log.warning("Skipping relative alignment (stage 1)...") + self.log.warning( + "Skipping relative alignment (stage 1)..." + ) else: self.log.warning("Skipping 'TweakRegStep'...") self.skip = True @@ -482,7 +495,9 @@ def process(self, input): # (typecasting numpy objects to python types so that it doesn't cause an # issue when saving datamodel to ASDF) wcs_fit_results = { - k: v.tolist() if isinstance(v, (np.ndarray, np.bool_)) else v + k: v.tolist() + if isinstance(v, (np.ndarray, np.bool_)) + else v for k, v in imcat.meta["fit_info"].items() } # add fit results and new WCS to datamodel @@ -529,10 +544,14 @@ def _imodel2wcsim(self, image_model): catalog_format = self.catalog_format else: catalog_format = "ascii.ecsv" - cat_name = str(catalog) + if isinstance(catalog, str): + # a string with the name of the catalog was provided catalog = Table.read(catalog, format=catalog_format) - catalog.meta["name"] = cat_name + + catalog.meta["name"] = ( + str(catalog) if isinstance(catalog, str) else model_name + ) except OSError: self.log.error(f"Cannot read catalog {catalog}") @@ -572,7 +591,9 @@ def _common_name(group): file_names = [] for im in group: if isinstance(im, rdm.DataModel): - file_names.append(os.path.splitext(im.meta.filename)[0].strip("_- ")) + file_names.append( + os.path.splitext(im.meta.filename)[0].strip("_- ") + ) else: raise TypeError("Input must be a list of datamodels list.") From a3098d504938e223e78761b8ec67e6dd5a356fcb Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Mon, 2 Oct 2023 10:25:49 -0400 Subject: [PATCH 02/13] Update TweakReg's regression test data. --- romancal/regtest/test_tweakreg.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/romancal/regtest/test_tweakreg.py b/romancal/regtest/test_tweakreg.py index d532d9bb8..a17580dd1 100644 --- a/romancal/regtest/test_tweakreg.py +++ b/romancal/regtest/test_tweakreg.py @@ -28,7 +28,7 @@ def create_asn_file(): "name": "files.asdf", "members": [ { - "expname": "r0000401001001001001_01101_0001_WFI01_cal_tweakreg.asdf", + "expname": "r0000501001001001001_01101_0001_WFI02_cal_tweakreg.asdf", "exptype": "science" } ] @@ -53,9 +53,9 @@ def test_tweakreg(rtdata, ignore_asdf_paths, tmp_path): # - assign_wcs; # - photom; # - source_detection. - input_data = "r0000401001001001001_01101_0001_WFI01_cal_tweakreg.asdf" - output_data = "r0000401001001001001_01101_0001_WFI01_output.asdf" - truth_data = "r0000401001001001001_01101_0001_WFI01_cal_twkreg_proc.asdf" + input_data = "r0000501001001001001_01101_0001_WFI02_cal_tweakreg.asdf" + output_data = "r0000501001001001001_01101_0001_WFI02_output.asdf" + truth_data = "r0000501001001001001_01101_0001_WFI02_tweakreg.asdf" rtdata.get_data(f"WFI/image/{input_data}") rtdata.get_truth(f"truth/WFI/image/{truth_data}") From 0c4a5e4323a387ad1226e7829de888208e7ccae7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 2 Oct 2023 20:38:29 +0000 Subject: [PATCH 03/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- romancal/tweakreg/tweakreg_step.py | 35 ++++++++---------------------- 1 file changed, 9 insertions(+), 26 deletions(-) diff --git a/romancal/tweakreg/tweakreg_step.py b/romancal/tweakreg/tweakreg_step.py index 0c3b4981e..fcf5f3b06 100644 --- a/romancal/tweakreg/tweakreg_step.py +++ b/romancal/tweakreg/tweakreg_step.py @@ -31,9 +31,7 @@ def _oxford_or_str_join(str_list): elif nelem == 2: return f"{str_list[0]} or {str_list[1]}" else: - return ( - ", ".join(map(repr, str_list[:-1])) + ", or " + repr(str_list[-1]) - ) + return ", ".join(map(repr, str_list[:-1])) + ", or " + repr(str_list[-1]) SINGLE_GROUP_REFCAT = ["GAIADR3", "GAIADR2", "GAIADR1"] @@ -146,9 +144,7 @@ def process(self, input): self.catalog_path = os.getcwd() self.catalog_path = Path(self.catalog_path).as_posix() - self.log.info( - f"All source catalogs will be saved to: {self.catalog_path}" - ) + self.log.info(f"All source catalogs will be saved to: {self.catalog_path}") if self.abs_refcat is None or len(self.abs_refcat.strip()) == 0: self.abs_refcat = DEFAULT_ABS_REFCAT @@ -240,9 +236,7 @@ def process(self, input): grp_img = list(images.models_grouped) self.log.info("") - self.log.info( - f"Number of image groups to be aligned: {len(grp_img):d}." - ) + self.log.info(f"Number of image groups to be aligned: {len(grp_img):d}.") self.log.info("Image groups:") if len(grp_img) == 1 and not ALIGN_TO_ABS_REFCAT: @@ -252,9 +246,7 @@ def process(self, input): self.log.info("") # we need at least two exposures to perform image alignment - self.log.warning( - "At least two exposures are required for image alignment." - ) + self.log.warning("At least two exposures are required for image alignment.") self.log.warning("Nothing to do. Skipping 'TweakRegStep'...") self.skip = True for model in images: @@ -328,8 +320,7 @@ def process(self, input): except ValueError as e: msg = e.args[0] if ( - msg - == "Too few input images (or groups of images) with non-empty" + msg == "Too few input images (or groups of images) with non-empty" " catalogs." ): # we need at least two exposures to perform image alignment @@ -337,9 +328,7 @@ def process(self, input): self.log.warning( "At least two exposures are required for image alignment." ) - self.log.warning( - "Nothing to do. Skipping 'TweakRegStep'..." - ) + self.log.warning("Nothing to do. Skipping 'TweakRegStep'...") for model in images: model.meta.cal_step["tweakreg"] = "SKIPPED" if not ALIGN_TO_ABS_REFCAT: @@ -365,9 +354,7 @@ def process(self, input): for model in images: model.meta.cal_step["tweakreg"] = "SKIPPED" if ALIGN_TO_ABS_REFCAT: - self.log.warning( - "Skipping relative alignment (stage 1)..." - ) + self.log.warning("Skipping relative alignment (stage 1)...") else: self.log.warning("Skipping 'TweakRegStep'...") self.skip = True @@ -495,9 +482,7 @@ def process(self, input): # (typecasting numpy objects to python types so that it doesn't cause an # issue when saving datamodel to ASDF) wcs_fit_results = { - k: v.tolist() - if isinstance(v, (np.ndarray, np.bool_)) - else v + k: v.tolist() if isinstance(v, (np.ndarray, np.bool_)) else v for k, v in imcat.meta["fit_info"].items() } # add fit results and new WCS to datamodel @@ -591,9 +576,7 @@ def _common_name(group): file_names = [] for im in group: if isinstance(im, rdm.DataModel): - file_names.append( - os.path.splitext(im.meta.filename)[0].strip("_- ") - ) + file_names.append(os.path.splitext(im.meta.filename)[0].strip("_- ")) else: raise TypeError("Input must be a list of datamodels list.") From 203619f3155f2bb4d4eda8a8023f548c9d078cc4 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 06:49:24 -0400 Subject: [PATCH 04/13] 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(): From d8c809eba071486c06e292eddbd454e86620eaea Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Oct 2023 10:57:43 +0000 Subject: [PATCH 05/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../tweakreg/tests/test_astrometric_utils.py | 18 +++++------------- 1 file changed, 5 insertions(+), 13 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index 4fa0ed56e..f26708a8e 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -8,9 +8,9 @@ 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 astropy.time import Time from gwcs import coordinate_frames as cf from gwcs import wcs from gwcs.geometry import CartesianToSpherical, SphericalToCartesian @@ -116,9 +116,7 @@ def get_parallax_correction(epoch, gaia_ref_epoch_coords): """ obs_date = Time(epoch, format="decimalyear") - earths_center_barycentric_coords = coord.get_body_barycentric( - "earth", obs_date - ) + 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 @@ -198,9 +196,7 @@ 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), @@ -458,9 +454,7 @@ 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) @@ -627,9 +621,7 @@ def test_get_catalog_using_epoch(ra, dec, epoch): gaia_ref_epoch=gaia_ref_epoch, ) # calculate parallax corrections - get_parallax_correction( - epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords - ) + get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) # calculate the expected coordinates value after corrections have been applied to # Gaia's reference epoch coordinates From 756b3f9d09c93118e86ed1b4e0a55b6667f53d07 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 07:00:02 -0400 Subject: [PATCH 06/13] Fix style checks. --- romancal/tweakreg/tests/test_astrometric_utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index f26708a8e..fa9d767ef 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -32,7 +32,8 @@ def __init__(self, *args, **kwargs): 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. + Calculates the proper motion correction for a given epoch and Gaia reference epoch + coordinates. Parameters ---------- From 07502dfaf74cafaba4ae69b2befe59cc1d59a429 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 09:25:59 -0400 Subject: [PATCH 07/13] Add atol parameter to compare_asdf method. --- romancal/regtest/test_tweakreg.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/romancal/regtest/test_tweakreg.py b/romancal/regtest/test_tweakreg.py index a17580dd1..7cfd30cf5 100644 --- a/romancal/regtest/test_tweakreg.py +++ b/romancal/regtest/test_tweakreg.py @@ -93,7 +93,7 @@ def test_tweakreg(rtdata, ignore_asdf_paths, tmp_path): ) assert "v2v3corr" in tweakreg_out.meta.wcs.available_frames - diff = compare_asdf(rtdata.output, rtdata.truth, **ignore_asdf_paths) + diff = compare_asdf(rtdata.output, rtdata.truth, atol=1e-3, **ignore_asdf_paths) step.log.info( f"DMS280 MSG: Was the proper TweakReg data produced? : {diff.identical}" ) From 31d294cb40589cb659ac8c4d6c23671d97344695 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 14:10:22 -0400 Subject: [PATCH 08/13] Use SCA 01 instead of SCA 02 in resample regtest. --- romancal/regtest/test_tweakreg.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/romancal/regtest/test_tweakreg.py b/romancal/regtest/test_tweakreg.py index 7cfd30cf5..a46bfdde2 100644 --- a/romancal/regtest/test_tweakreg.py +++ b/romancal/regtest/test_tweakreg.py @@ -28,7 +28,7 @@ def create_asn_file(): "name": "files.asdf", "members": [ { - "expname": "r0000501001001001001_01101_0001_WFI02_cal_tweakreg.asdf", + "expname": "r0000501001001001001_01101_0001_WFI01_cal_tweakreg.asdf", "exptype": "science" } ] @@ -53,9 +53,9 @@ def test_tweakreg(rtdata, ignore_asdf_paths, tmp_path): # - assign_wcs; # - photom; # - source_detection. - input_data = "r0000501001001001001_01101_0001_WFI02_cal_tweakreg.asdf" - output_data = "r0000501001001001001_01101_0001_WFI02_output.asdf" - truth_data = "r0000501001001001001_01101_0001_WFI02_tweakreg.asdf" + input_data = "r0000501001001001001_01101_0001_WFI01_cal_tweakreg.asdf" + output_data = "r0000501001001001001_01101_0001_WFI01_output.asdf" + truth_data = "r0000501001001001001_01101_0001_WFI01_tweakreg.asdf" rtdata.get_data(f"WFI/image/{input_data}") rtdata.get_truth(f"truth/WFI/image/{truth_data}") From 288647173574200484a59a791c494c72269c8939 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 14:14:48 -0400 Subject: [PATCH 09/13] Set atol=1e-5 and rtol=0 in astrometric_utils unit tests. --- .../tweakreg/tests/test_astrometric_utils.py | 24 ++++++++++++------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index fa9d767ef..6b0388993 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -117,7 +117,9 @@ def get_parallax_correction(epoch, gaia_ref_epoch_coords): """ obs_date = Time(epoch, format="decimalyear") - earths_center_barycentric_coords = coord.get_body_barycentric("earth", obs_date) + 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 @@ -197,7 +199,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), @@ -455,7 +459,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) @@ -622,7 +628,9 @@ def test_get_catalog_using_epoch(ra, dec, epoch): gaia_ref_epoch=gaia_ref_epoch, ) # calculate parallax corrections - get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) + get_parallax_correction( + epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords + ) # calculate the expected coordinates value after corrections have been applied to # Gaia's reference epoch coordinates @@ -638,11 +646,9 @@ def test_get_catalog_using_epoch(ra, dec, epoch): ) assert len(result) > 0 - # 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() + + assert np.isclose(returned_ra, expected_ra, atol=1e-5, rtol=0).all() + assert np.isclose(returned_dec, expected_dec, atol=1e-5, rtol=0).all() def test_get_catalog_timeout(): From 5bc456022891a182a739f7e618d5aad12f8e3aab Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Thu, 5 Oct 2023 18:15:34 +0000 Subject: [PATCH 10/13] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- .../tweakreg/tests/test_astrometric_utils.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index 6b0388993..b9bfcb19c 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -117,9 +117,7 @@ def get_parallax_correction(epoch, gaia_ref_epoch_coords): """ obs_date = Time(epoch, format="decimalyear") - earths_center_barycentric_coords = coord.get_body_barycentric( - "earth", obs_date - ) + 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 @@ -199,9 +197,7 @@ 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), @@ -459,9 +455,7 @@ 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) @@ -628,9 +622,7 @@ def test_get_catalog_using_epoch(ra, dec, epoch): gaia_ref_epoch=gaia_ref_epoch, ) # calculate parallax corrections - get_parallax_correction( - epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords - ) + get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) # calculate the expected coordinates value after corrections have been applied to # Gaia's reference epoch coordinates From 75b7857e211e139d7e5eac04dd1a37d887fc0edd Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Wed, 11 Oct 2023 19:00:53 -0400 Subject: [PATCH 11/13] Fix astrometric_utils tests. --- .../tweakreg/tests/test_astrometric_utils.py | 248 +++++++++++++----- 1 file changed, 189 insertions(+), 59 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index b9bfcb19c..3e8901c0e 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -24,12 +24,160 @@ get_catalog, ) +ARAD = np.pi / 180.0 + class MockConnectionError: def __init__(self, *args, **kwargs): raise requests.exceptions.ConnectionError +def get_earth_sun_orbit_elem(epoch): + amjd = (epoch - 2012.0) * 365.25 + 55927.0 + imjd = int(amjd) + mjd = float(imjd) + + time_argument = mjd - 51544.5 + fmean_longitude = 280.460 + 0.9856474 * time_argument + fmean_anomaly = 357.528 + 0.9856003 * time_argument + iremainder = int(fmean_longitude / 360.0) + fmean_longitude = fmean_longitude - (iremainder * 360.0) + if fmean_longitude < 0: + fmean_longitude = fmean_longitude + 360.0 + iremainder = int(fmean_anomaly / 360.0) + fmean_anomaly = fmean_anomaly - (iremainder * 360.0) + if fmean_anomaly < 0: + fmean_anomaly = fmean_anomaly + 360.0 + ecliptic_longitude = ( + fmean_longitude + + 1.915 * np.sin(fmean_anomaly * ARAD) + + 0.02 * np.sin(2.0 * fmean_anomaly * ARAD) + ) + ecliptic_obliquity = 23.439 - 4.0 * 0.0000001 * time_argument + + # radius vector (in AU) + # (approximation using an algorithm from USNO) + earth_sun_distance = ( + 1.00014 + - 0.01671 * np.cos(fmean_anomaly * ARAD) + - 0.00014 * np.cos(2.0 * fmean_anomaly * ARAD) + ) + + earth_sun_distance = u.Quantity(earth_sun_distance, unit="AU") + + return { + "earth_sun_distance": earth_sun_distance, + "ecliptic_longitude": ecliptic_longitude, + "ecliptic_obliquity": ecliptic_obliquity, + } + + +def get_parallax_correction_mast(epoch, gaia_ref_epoch_coords): + orbit_elem = get_earth_sun_orbit_elem(epoch) + earth_sun_distance = orbit_elem["earth_sun_distance"] + ecliptic_longitude = orbit_elem["ecliptic_longitude"] + ecliptic_obliquity = orbit_elem["ecliptic_obliquity"] + + # cartesian coordinates of the sun + xsun = earth_sun_distance * np.cos(ecliptic_longitude * ARAD) + ysun = ( + earth_sun_distance + * np.cos(ecliptic_obliquity * ARAD) + * np.sin(ecliptic_longitude * ARAD) + ) + zsun = ( + earth_sun_distance + * np.sin(ecliptic_obliquity * ARAD) + * np.sin(ecliptic_longitude * ARAD) + ) + + # angular displacement components + delta_ra = ( + u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas") + / np.cos(gaia_ref_epoch_coords["dec"] * ARAD) + * ( + -xsun.value * np.sin(gaia_ref_epoch_coords["ra"] * ARAD) + + ysun.value * np.cos(gaia_ref_epoch_coords["ra"] * ARAD) + ) + ).to("deg") + delta_dec = ( + u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas") + * ( + -xsun.value + * np.cos(gaia_ref_epoch_coords["ra"] * ARAD) + * np.sin(gaia_ref_epoch_coords["dec"] * ARAD) + - ysun.value + * np.sin(gaia_ref_epoch_coords["ra"] * ARAD) + * np.sin(gaia_ref_epoch_coords["dec"] * ARAD) + + zsun.value * np.cos(gaia_ref_epoch_coords["dec"] * ARAD) + ) + ).to("deg") + + return delta_ra, delta_dec + + +def get_parallax_correction_barycenter(epoch, gaia_ref_epoch_coords): + """ + Calculates the parallax correction in the Earth barycenter frame for a given epoch + and Gaia reference epoch coordinates (i.e. Gaia coordinates at the reference epoch). + + Parameters + ---------- + epoch : float + The epoch for which the parallax correction is calculated. + gaia_ref_epoch_coords : dict + The Gaia reference epoch coordinates, including 'ra', 'dec', and 'parallax'. + + Returns + ------- + tuple + A tuple containing the delta_ra and delta_dec values of the parallax correction + in degrees. + + Examples + -------- + .. code-block :: python + epoch = 2022.5 + gaia_coords = {'ra': 180.0, 'dec': 45.0, 'parallax': 10.0} + correction = get_parallax_correction_earth_barycenter(epoch, gaia_coords) + print(correction) + (0.001, -0.002) + """ # noqa: E501 + + obs_date = Time(epoch, format="decimalyear") + earths_center_barycentric_coords = coord.get_body_barycentric( + "earth", obs_date, ephemeris="builtin" + ) + 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) + delta_ra = ( + u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad) + * (1 / np.cos(gaia_ref_epoch_coords["dec"] * ARAD)) + * ( + earth_X.value * np.sin(gaia_ref_epoch_coords["ra"] * ARAD) + - earth_Y.value * np.cos(gaia_ref_epoch_coords["ra"] * ARAD) + ) + ).to("deg") + delta_dec = ( + u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas").to(u.rad) + * ( + earth_X.value + * np.cos(gaia_ref_epoch_coords["ra"] * ARAD) + * np.sin(gaia_ref_epoch_coords["dec"] * ARAD) + + earth_Y.value + * np.sin(gaia_ref_epoch_coords["ra"] * ARAD) + * np.sin(gaia_ref_epoch_coords["dec"] * ARAD) + - earth_Z.value * np.cos(gaia_ref_epoch_coords["dec"] * ARAD) + ) + ).to("deg") + + return delta_ra, delta_dec + + 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 @@ -88,62 +236,22 @@ def get_proper_motion_correction(epoch, gaia_ref_epoch_coords, gaia_ref_epoch): 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) - """ + # get parallax correction using textbook calculations (i.e. Earth's barycenter) + parallax_corr = get_parallax_correction_barycenter( + epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_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 + # get parallax using the same coordinates frame as MAST + # (i.e. Sun's geocentric coordinates) + parallax_corr_mast = get_parallax_correction_mast( + epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords + ) - # 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") + # add parallax corrections columns to the main table + gaia_ref_epoch_coords["parallax_delta_ra"] = parallax_corr[0] + gaia_ref_epoch_coords["parallax_delta_dec"] = parallax_corr[1] + gaia_ref_epoch_coords["parallax_delta_ra_mast"] = parallax_corr_mast[0] + gaia_ref_epoch_coords["parallax_delta_dec_mast"] = parallax_corr_mast[1] def update_wcsinfo(input_dm): @@ -198,6 +306,7 @@ def create_wcs_for_tweakreg_pipeline(input_dm, shift_1=0, shift_2=0): # 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), @@ -455,6 +564,7 @@ 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" @@ -591,12 +701,11 @@ def test_get_catalog_valid_parameters_but_no_sources_returned(): (0, 0, 2000), (0, 0, 2010.3), (0, 0, 2030), - (269.4521, 4.6933, 2030), - (89, 80, 2010), ], ) def test_get_catalog_using_epoch(ra, dec, epoch): - """Test that get_catalog returns coordinates corrected by proper motion.""" + """Test that get_catalog returns coordinates corrected by proper motion + and parallax.""" result_all = get_catalog(ra, dec, epoch=epoch) @@ -606,6 +715,8 @@ def test_get_catalog_using_epoch(ra, dec, epoch): mask = result_all["parallax_over_error"] > 5 result = result_all[mask] + + # updated coordinates at the provided epoch returned_ra = np.array(result["ra"]) returned_dec = np.array(result["dec"]) @@ -623,9 +734,12 @@ def test_get_catalog_using_epoch(ra, dec, epoch): ) # calculate parallax corrections get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) + get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) # calculate the expected coordinates value after corrections have been applied to # Gaia's reference epoch coordinates + + # textbook (barycentric frame) expected_ra = ( gaia_ref_epoch_coords["ra"] + gaia_ref_epoch_coords["pm_delta_ra"] @@ -637,10 +751,26 @@ def test_get_catalog_using_epoch(ra, dec, epoch): + gaia_ref_epoch_coords["parallax_delta_dec"] ) + # mast (geocentric frame) + expected_ra_mast = ( + gaia_ref_epoch_coords["ra"] + + gaia_ref_epoch_coords["pm_delta_ra"] + + gaia_ref_epoch_coords["parallax_delta_ra_mast"] + ) + expected_dec_mast = ( + gaia_ref_epoch_coords["dec"] + + gaia_ref_epoch_coords["pm_delta_dec"] + + gaia_ref_epoch_coords["parallax_delta_dec_mast"] + ) + assert len(result) > 0 - assert np.isclose(returned_ra, expected_ra, atol=1e-5, rtol=0).all() - assert np.isclose(returned_dec, expected_dec, atol=1e-5, rtol=0).all() + # N.B.: atol=1e-8 (in deg) corresponds to a coordinate difference of ~ 40 uas + assert np.isclose(expected_ra, returned_ra, atol=1e-8, rtol=0).all() + assert np.isclose(expected_dec, returned_dec, atol=1e-8, rtol=0).all() + + assert np.isclose(expected_ra_mast, returned_ra, atol=5e-10, rtol=0).all() + assert np.isclose(expected_dec_mast, returned_dec, atol=5e-10, rtol=0).all() def test_get_catalog_timeout(): From 2566277067e36b49afa3f7674eca62da1e2605f7 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Wed, 11 Oct 2023 19:28:35 -0400 Subject: [PATCH 12/13] Fix check-styles. --- .../tweakreg/tests/test_astrometric_utils.py | 142 +++++++++++++++++- 1 file changed, 138 insertions(+), 4 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index 3e8901c0e..82f316ea1 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -33,6 +33,41 @@ def __init__(self, *args, **kwargs): def get_earth_sun_orbit_elem(epoch): + """ + Calculates the Earth-Sun orbit elements for a given epoch. + + Parameters + ---------- + epoch : float + The epoch for which to calculate the Earth-Sun orbit elements. + + Returns + ------- + dict + A dictionary containing the following orbit elements: + - "earth_sun_distance" : `Quantity` + The Earth-Sun distance in astronomical units (AU). + - "ecliptic_longitude" : float + The ecliptic longitude in degrees. + - "ecliptic_obliquity" : float + The ecliptic obliquity in degrees. + + Notes + ----- + This function calculates various parameters related to the Earth-Sun orbit based on the provided epoch. + + Examples + -------- + .. code-block:: python + + epoch = 2023.5 + orbit_elements = get_earth_sun_orbit_elem(epoch) + print(orbit_elements) + + # Output: + # {'earth_sun_distance': , 'ecliptic_longitude': 189.123 degrees, 'ecliptic_obliquity': 23.437 degrees} + """ # noqa: E501 + amjd = (epoch - 2012.0) * 365.25 + 55927.0 imjd = int(amjd) mjd = float(imjd) @@ -73,6 +108,55 @@ def get_earth_sun_orbit_elem(epoch): def get_parallax_correction_mast(epoch, gaia_ref_epoch_coords): + """ + Calculates the parallax correction for MAST coordinates based on the Earth-Sun orbit elements and Gaia reference epoch coordinates. + + Parameters + ---------- + epoch : float + The epoch for which to calculate the parallax correction. + gaia_ref_epoch_coords : dict + A dictionary containing the Gaia reference epoch coordinates: + - "ra" : float + The right ascension in degrees. + - "dec" : float + The declination in degrees. + - "parallax" : float + The parallax in milliarcseconds (mas). + + Returns + ------- + tuple + A tuple containing the following parallax correction components: + - delta_ra : `Quantity` + The correction in right ascension in degrees. + - delta_dec : `Quantity` + The correction in declination in degrees. + + Notes + ----- + This function calculates the parallax correction for MAST coordinates based on the + Earth-Sun orbit elements and Gaia reference epoch coordinates. It uses the + Earth-Sun distance, ecliptic longitude, and ecliptic obliquity obtained from + the `get_earth_sun_orbit_elem` function. + + Examples + -------- + .. code-block:: python + + epoch = 2023.5 + gaia_coords = { + "ra": 180.0, + "dec": 30.0, + "parallax": 2.5 + } + delta_ra, delta_dec = get_parallax_correction_mast(epoch, gaia_coords) + print(delta_ra, delta_dec) + + # Output: + # -0.001234 deg, 0.002345 deg + """ # noqa: E501 + orbit_elem = get_earth_sun_orbit_elem(epoch) earth_sun_distance = orbit_elem["earth_sun_distance"] ecliptic_longitude = orbit_elem["ecliptic_longitude"] @@ -141,7 +225,7 @@ def get_parallax_correction_barycenter(epoch, gaia_ref_epoch_coords): gaia_coords = {'ra': 180.0, 'dec': 45.0, 'parallax': 10.0} correction = get_parallax_correction_earth_barycenter(epoch, gaia_coords) print(correction) - (0.001, -0.002) + (0.001, -0.002) """ # noqa: E501 obs_date = Time(epoch, format="decimalyear") @@ -208,7 +292,7 @@ def get_proper_motion_correction(epoch, gaia_ref_epoch_coords, gaia_ref_epoch): } gaia_ref_epoch = 2020.0 get_proper_motion_correction(epoch, gaia_coords, gaia_ref_epoch) - """ + """ # noqa: E501 expected_new_dec = ( np.array( @@ -236,6 +320,50 @@ def get_proper_motion_correction(epoch, gaia_ref_epoch_coords, gaia_ref_epoch): def get_parallax_correction(epoch, gaia_ref_epoch_coords): + """ + Calculates the parallax correction for a given epoch and Gaia reference epoch + coordinates. + + Parameters + ---------- + epoch : float + The epoch for which to calculate the parallax correction. + gaia_ref_epoch_coords : dict + A dictionary containing the Gaia reference epoch coordinates: + - "ra" : float + The right ascension in degrees. + - "dec" : float + The declination in degrees. + - "parallax" : float + The parallax in milliarcseconds (mas). + + Returns + ------- + None + + Notes + ----- + This function calculates the parallax correction for a given epoch and Gaia + reference epoch coordinates. It uses the `get_parallax_correction_barycenter` + and `get_parallax_correction_mast` functions to obtain the parallax corrections + based on different coordinate frames. + + Examples + -------- + This function is typically used to add parallax correction columns to a main table + of Gaia reference epoch coordinates. + + .. code-block:: python + + epoch = 2023.5 + gaia_coords = { + "ra": 180.0, + "dec": 30.0, + "parallax": 2.5 + } + get_parallax_correction(epoch, gaia_coords) + """ # noqa: E501 + # get parallax correction using textbook calculations (i.e. Earth's barycenter) parallax_corr = get_parallax_correction_barycenter( epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords @@ -564,7 +692,6 @@ 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" @@ -705,7 +832,14 @@ 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 - and parallax.""" + and parallax. The idea is to fetch data for a specific epoch from the MAST VO API + and compare them with the expected coordinates for that epoch. + First, the data for a specific coordinates and epoch are fetched from the MAST VO + API. Then, the data for the same coordinates are fetched for the Gaia's reference + epoch of 2016.0, and corrected for proper motion and parallax using explicit + calculations for the initially specified epoch. We then compare the results between + the returned coordinates from the MAST VO API and the manually updated + coordinates.""" result_all = get_catalog(ra, dec, epoch=epoch) From 768dadec30e37f7e8f12757ba60fbd185f2caf1e Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 12 Oct 2023 11:03:45 -0400 Subject: [PATCH 13/13] Remove MAST code and decrease catalog comparison tolerance. --- .../tweakreg/tests/test_astrometric_utils.py | 211 +----------------- 1 file changed, 8 insertions(+), 203 deletions(-) diff --git a/romancal/tweakreg/tests/test_astrometric_utils.py b/romancal/tweakreg/tests/test_astrometric_utils.py index 82f316ea1..18f27988d 100644 --- a/romancal/tweakreg/tests/test_astrometric_utils.py +++ b/romancal/tweakreg/tests/test_astrometric_utils.py @@ -10,6 +10,7 @@ from astropy import units as u from astropy.modeling import models from astropy.modeling.models import RotationSequence3D, Scale, Shift +from astropy.stats import mad_std from astropy.time import Time from gwcs import coordinate_frames as cf from gwcs import wcs @@ -32,174 +33,6 @@ def __init__(self, *args, **kwargs): raise requests.exceptions.ConnectionError -def get_earth_sun_orbit_elem(epoch): - """ - Calculates the Earth-Sun orbit elements for a given epoch. - - Parameters - ---------- - epoch : float - The epoch for which to calculate the Earth-Sun orbit elements. - - Returns - ------- - dict - A dictionary containing the following orbit elements: - - "earth_sun_distance" : `Quantity` - The Earth-Sun distance in astronomical units (AU). - - "ecliptic_longitude" : float - The ecliptic longitude in degrees. - - "ecliptic_obliquity" : float - The ecliptic obliquity in degrees. - - Notes - ----- - This function calculates various parameters related to the Earth-Sun orbit based on the provided epoch. - - Examples - -------- - .. code-block:: python - - epoch = 2023.5 - orbit_elements = get_earth_sun_orbit_elem(epoch) - print(orbit_elements) - - # Output: - # {'earth_sun_distance': , 'ecliptic_longitude': 189.123 degrees, 'ecliptic_obliquity': 23.437 degrees} - """ # noqa: E501 - - amjd = (epoch - 2012.0) * 365.25 + 55927.0 - imjd = int(amjd) - mjd = float(imjd) - - time_argument = mjd - 51544.5 - fmean_longitude = 280.460 + 0.9856474 * time_argument - fmean_anomaly = 357.528 + 0.9856003 * time_argument - iremainder = int(fmean_longitude / 360.0) - fmean_longitude = fmean_longitude - (iremainder * 360.0) - if fmean_longitude < 0: - fmean_longitude = fmean_longitude + 360.0 - iremainder = int(fmean_anomaly / 360.0) - fmean_anomaly = fmean_anomaly - (iremainder * 360.0) - if fmean_anomaly < 0: - fmean_anomaly = fmean_anomaly + 360.0 - ecliptic_longitude = ( - fmean_longitude - + 1.915 * np.sin(fmean_anomaly * ARAD) - + 0.02 * np.sin(2.0 * fmean_anomaly * ARAD) - ) - ecliptic_obliquity = 23.439 - 4.0 * 0.0000001 * time_argument - - # radius vector (in AU) - # (approximation using an algorithm from USNO) - earth_sun_distance = ( - 1.00014 - - 0.01671 * np.cos(fmean_anomaly * ARAD) - - 0.00014 * np.cos(2.0 * fmean_anomaly * ARAD) - ) - - earth_sun_distance = u.Quantity(earth_sun_distance, unit="AU") - - return { - "earth_sun_distance": earth_sun_distance, - "ecliptic_longitude": ecliptic_longitude, - "ecliptic_obliquity": ecliptic_obliquity, - } - - -def get_parallax_correction_mast(epoch, gaia_ref_epoch_coords): - """ - Calculates the parallax correction for MAST coordinates based on the Earth-Sun orbit elements and Gaia reference epoch coordinates. - - Parameters - ---------- - epoch : float - The epoch for which to calculate the parallax correction. - gaia_ref_epoch_coords : dict - A dictionary containing the Gaia reference epoch coordinates: - - "ra" : float - The right ascension in degrees. - - "dec" : float - The declination in degrees. - - "parallax" : float - The parallax in milliarcseconds (mas). - - Returns - ------- - tuple - A tuple containing the following parallax correction components: - - delta_ra : `Quantity` - The correction in right ascension in degrees. - - delta_dec : `Quantity` - The correction in declination in degrees. - - Notes - ----- - This function calculates the parallax correction for MAST coordinates based on the - Earth-Sun orbit elements and Gaia reference epoch coordinates. It uses the - Earth-Sun distance, ecliptic longitude, and ecliptic obliquity obtained from - the `get_earth_sun_orbit_elem` function. - - Examples - -------- - .. code-block:: python - - epoch = 2023.5 - gaia_coords = { - "ra": 180.0, - "dec": 30.0, - "parallax": 2.5 - } - delta_ra, delta_dec = get_parallax_correction_mast(epoch, gaia_coords) - print(delta_ra, delta_dec) - - # Output: - # -0.001234 deg, 0.002345 deg - """ # noqa: E501 - - orbit_elem = get_earth_sun_orbit_elem(epoch) - earth_sun_distance = orbit_elem["earth_sun_distance"] - ecliptic_longitude = orbit_elem["ecliptic_longitude"] - ecliptic_obliquity = orbit_elem["ecliptic_obliquity"] - - # cartesian coordinates of the sun - xsun = earth_sun_distance * np.cos(ecliptic_longitude * ARAD) - ysun = ( - earth_sun_distance - * np.cos(ecliptic_obliquity * ARAD) - * np.sin(ecliptic_longitude * ARAD) - ) - zsun = ( - earth_sun_distance - * np.sin(ecliptic_obliquity * ARAD) - * np.sin(ecliptic_longitude * ARAD) - ) - - # angular displacement components - delta_ra = ( - u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas") - / np.cos(gaia_ref_epoch_coords["dec"] * ARAD) - * ( - -xsun.value * np.sin(gaia_ref_epoch_coords["ra"] * ARAD) - + ysun.value * np.cos(gaia_ref_epoch_coords["ra"] * ARAD) - ) - ).to("deg") - delta_dec = ( - u.Quantity(gaia_ref_epoch_coords["parallax"], unit="mas") - * ( - -xsun.value - * np.cos(gaia_ref_epoch_coords["ra"] * ARAD) - * np.sin(gaia_ref_epoch_coords["dec"] * ARAD) - - ysun.value - * np.sin(gaia_ref_epoch_coords["ra"] * ARAD) - * np.sin(gaia_ref_epoch_coords["dec"] * ARAD) - + zsun.value * np.cos(gaia_ref_epoch_coords["dec"] * ARAD) - ) - ).to("deg") - - return delta_ra, delta_dec - - def get_parallax_correction_barycenter(epoch, gaia_ref_epoch_coords): """ Calculates the parallax correction in the Earth barycenter frame for a given epoch @@ -369,17 +202,9 @@ def get_parallax_correction(epoch, gaia_ref_epoch_coords): epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords ) - # get parallax using the same coordinates frame as MAST - # (i.e. Sun's geocentric coordinates) - parallax_corr_mast = get_parallax_correction_mast( - epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords - ) - # add parallax corrections columns to the main table gaia_ref_epoch_coords["parallax_delta_ra"] = parallax_corr[0] gaia_ref_epoch_coords["parallax_delta_dec"] = parallax_corr[1] - gaia_ref_epoch_coords["parallax_delta_ra_mast"] = parallax_corr_mast[0] - gaia_ref_epoch_coords["parallax_delta_dec_mast"] = parallax_corr_mast[1] def update_wcsinfo(input_dm): @@ -841,14 +666,7 @@ def test_get_catalog_using_epoch(ra, dec, epoch): the returned coordinates from the MAST VO API and the manually updated coordinates.""" - 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] + result = get_catalog(ra, dec, epoch=epoch) # updated coordinates at the provided epoch returned_ra = np.array(result["ra"]) @@ -858,7 +676,7 @@ def test_get_catalog_using_epoch(ra, dec, epoch): gaia_ref_epoch = 2016.0 gaia_ref_epoch_coords_all = get_catalog(ra, dec, epoch=gaia_ref_epoch) - gaia_ref_epoch_coords = gaia_ref_epoch_coords_all[mask] + gaia_ref_epoch_coords = gaia_ref_epoch_coords_all # [mask] # calculate proper motion corrections get_proper_motion_correction( @@ -868,7 +686,6 @@ def test_get_catalog_using_epoch(ra, dec, epoch): ) # calculate parallax corrections get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) - get_parallax_correction(epoch=epoch, gaia_ref_epoch_coords=gaia_ref_epoch_coords) # calculate the expected coordinates value after corrections have been applied to # Gaia's reference epoch coordinates @@ -885,26 +702,14 @@ def test_get_catalog_using_epoch(ra, dec, epoch): + gaia_ref_epoch_coords["parallax_delta_dec"] ) - # mast (geocentric frame) - expected_ra_mast = ( - gaia_ref_epoch_coords["ra"] - + gaia_ref_epoch_coords["pm_delta_ra"] - + gaia_ref_epoch_coords["parallax_delta_ra_mast"] - ) - expected_dec_mast = ( - gaia_ref_epoch_coords["dec"] - + gaia_ref_epoch_coords["pm_delta_dec"] - + gaia_ref_epoch_coords["parallax_delta_dec_mast"] - ) - assert len(result) > 0 - # N.B.: atol=1e-8 (in deg) corresponds to a coordinate difference of ~ 40 uas - assert np.isclose(expected_ra, returned_ra, atol=1e-8, rtol=0).all() - assert np.isclose(expected_dec, returned_dec, atol=1e-8, rtol=0).all() + # adopted tolerance: 2.8e-9 deg -> 10 uas (~0.0001 pix) + assert np.median(returned_ra - expected_ra) < 2.8e-9 + assert np.median(returned_dec - expected_dec) < 2.8e-9 - assert np.isclose(expected_ra_mast, returned_ra, atol=5e-10, rtol=0).all() - assert np.isclose(expected_dec_mast, returned_dec, atol=5e-10, rtol=0).all() + assert mad_std(returned_ra - expected_ra) < 2.8e-9 + assert mad_std(returned_dec - expected_dec) < 2.8e-9 def test_get_catalog_timeout():