From e835fdb6c3a7d146ab8f5f8ef600743fb4fb3cb4 Mon Sep 17 00:00:00 2001 From: D Davis <49163225+ddavis-stsci@users.noreply.github.com> Date: Wed, 11 Oct 2023 09:42:03 -0400 Subject: [PATCH 01/16] RCAL-641 Add FOV association generation (#931) Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> --- CHANGES.rst | 5 ++ romancal/associations/generate.py | 11 ++- romancal/associations/lib/dms_base.py | 87 ++++++------------- romancal/associations/lib/rules_elpp_base.py | 84 +++++++++++++++++- romancal/associations/lib/rules_level2.py | 46 ++++++++-- .../associations/tests/test_level2_basics.py | 5 +- .../tests/test_level2_candidates.py | 4 +- 7 files changed, 165 insertions(+), 77 deletions(-) diff --git a/CHANGES.rst b/CHANGES.rst index 2f88c6aa7..9de551250 100644 --- a/CHANGES.rst +++ b/CHANGES.rst @@ -1,6 +1,11 @@ 0.13.0 (unreleased) =================== +associations +------------ + +- Add FOV associations to the code [#931] + general ------- diff --git a/romancal/associations/generate.py b/romancal/associations/generate.py index b7184da9d..ed35dff30 100644 --- a/romancal/associations/generate.py +++ b/romancal/associations/generate.py @@ -98,10 +98,13 @@ def generate(pool, rules, version_id=None, finalize=True): # Finalize found associations logger.debug("# associations before finalization: %d", len(associations)) - try: - finalized_asns = rules.callback.reduce("finalize", associations) - except KeyError: - finalized_asns = associations + finalized_asns = associations + if finalize: + logger.debug("Performing association finalization.") + try: + finalized_asns = rules.callback.reduce("finalize", associations) + except KeyError as exception: + logger.debug("Finalization failed for reason: %s", exception) return finalized_asns diff --git a/romancal/associations/lib/dms_base.py b/romancal/associations/lib/dms_base.py index bf10623ec..88bde16e3 100644 --- a/romancal/associations/lib/dms_base.py +++ b/romancal/associations/lib/dms_base.py @@ -21,96 +21,59 @@ # Acquisition and Confirmation images ACQ_EXP_TYPES = ( - "mir_tacq", - "mir_taconfirm", - "nis_taconfirm", - "nis_tacq", "nrc_taconfirm", "nrc_tacq", - "nrs_confirm", - "nrs_msata", - "nrs_taconfirm", - "nrs_tacq", - "nrs_taslit", - "nrs_verify", - "nrs_wata", ) # Exposure EXP_TYPE to Association EXPTYPE mapping # flake8: noqa: E241 EXPTYPE_MAP = { - "mir_darkall": "dark", - "mir_darkimg": "dark", - "mir_darkmrs": "dark", - "mir_flatimage": "flat", - "mir_flatmrs": "flat", - "mir_flatimage-ext": "flat", - "mir_flatmrs-ext": "flat", - "mir_tacq": "target_acquisition", - "mir_taconfirm": "target_acquisition", - "nis_dark": "dark", - "nis_focus": "engineering", - "nis_lamp": "engineering", - "nis_tacq": "target_acquisition", - "nis_taconfirm": "target_acquisition", "nrc_dark": "dark", "nrc_flat": "flat", "nrc_focus": "engineering", "nrc_led": "engineering", "nrc_tacq": "target_acquisition", "nrc_taconfirm": "target_acquisition", - "nrs_autoflat": "autoflat", - "nrs_autowave": "autowave", - "nrs_confirm": "target_acquisition", - "nrs_dark": "dark", - "nrs_focus": "engineering", - "nrs_image": "engineering", - "nrs_lamp": "engineering", - "nrs_msata": "target_acquisition", - "nrs_tacq": "target_acquisition", - "nrs_taconfirm": "target_acquisition", - "nrs_taslit": "target_acquisition", - "nrs_wata": "target_acquisition", } # Coronographic exposures CORON_EXP_TYPES = ["mir_4qpm", "mir_lyot", "nrc_coron"] +# Roman WFI detectors +WFI_DETECTORS = [ + "wfi01", + "wfi02", + "wfi03", + "wfi04", + "wfi05", + "wfi06", + "wfi07", + "wfi08", + "wfi09", + "wfi10", + "wfi11", + "wfi12", + "wfi13", + "wfi14", + "wfi15", + "wfi16", + "wfi17", + "wfi18", +] + # Exposures that get Level2b processing IMAGE2_SCIENCE_EXP_TYPES = [ "wfi_image", - "mir_4qpm", - "mir_image", - "mir_lyot", - "nis_ami", - "nis_image", - "nrc_coron", - "nrc_image", - "nrs_mimf", - "nrc_tsimage", ] IMAGE2_NONSCIENCE_EXP_TYPES = [ - "mir_coroncal", - "nis_focus", - "nrc_focus", - "nrs_focus", - "nrs_image", + "wfi_focus", ] IMAGE2_NONSCIENCE_EXP_TYPES.extend(ACQ_EXP_TYPES) SPEC2_SCIENCE_EXP_TYPES = [ - "mir_lrs-fixedslit", - "mir_lrs-slitless", - "mir_mrs", - "nis_soss", - "nis_wfss", - "nrc_tsgrism", - "nrc_wfss", - "nrs_fixedslit", - "nrs_ifu", - "nrs_msaspec", - "nrs_brightobj", + "wfi_grism", + "wfi_prism", ] SPECIAL_EXPOSURE_MODIFIERS = { diff --git a/romancal/associations/lib/rules_elpp_base.py b/romancal/associations/lib/rules_elpp_base.py index 8cbf7d87d..cf5917ea7 100644 --- a/romancal/associations/lib/rules_elpp_base.py +++ b/romancal/associations/lib/rules_elpp_base.py @@ -1,9 +1,10 @@ """Base classes which define the ELPP Associations""" +import copy import logging import re from collections import defaultdict -from os.path import basename +from os.path import basename, split, splitext from stpipe.format_template import FormatTemplate @@ -18,6 +19,7 @@ IMAGE2_NONSCIENCE_EXP_TYPES, IMAGE2_SCIENCE_EXP_TYPES, SPEC2_SCIENCE_EXP_TYPES, + WFI_DETECTORS, DMSAttrConstraint, DMSBaseMixin, ) @@ -36,6 +38,7 @@ "AsnMixin_AuxData", "AsnMixin_Science", "AsnMixin_Spectrum", + "AsnMixin_Lv2FOV", "AsnMixin_Lv2Image", "AsnMixin_Lv2GBTDSfull", "AsnMixin_Lv2GBTDSpass", @@ -56,6 +59,7 @@ "Constraint_Single_Science", "Constraint_Spectral_Science", "Constraint_Target", + "Constraint_Filename", "DMS_ELPP_Base", "DMSAttrConstraint", "ProcessList", @@ -318,6 +322,43 @@ def make_member(self, item): ) return member + def make_fov_asn(self): + """Take the association with an single exposure with _WFI_ in the name + and expand that to include all 18 detectors. + + Returns + ------- + associations : [association[, ...]] + List of new members to be used in place of + the current one. + """ + results = [] + + # expand the products from _wfi_ to _wfi{det}_ + for product in self["products"]: + for member in product["members"]: + asn = copy.deepcopy(self) + asn.data["products"] = None + product_name = ( + splitext( + split(self.data["products"][0]["members"][0]["expname"])[1] + )[0].rsplit("_", 1)[0] + + "_drzl" + ) + asn.new_product(product_name) + new_members = asn.current_product["members"] + if "_wfi_" in member["expname"]: + # Make and add a member for each detector + for det in WFI_DETECTORS: + new_member = copy.deepcopy(member) + new_member["expname"] = member["expname"].replace("wfi", det) + new_members.append(new_member) + if asn.is_valid: + results.append(asn) + return results + else: + return None + def _init_hook(self, item): """Post-check and pre-add initialization""" super()._init_hook(item) @@ -646,6 +687,16 @@ def __init__(self): ) +class Constraint_Filename(DMSAttrConstraint): + """Select on visit number""" + + def __init__(self): + super().__init__( + name="Filename", + sources=["filename"], + ) + + class Constraint_Expos(DMSAttrConstraint): """Select on exposure number""" @@ -653,8 +704,8 @@ def __init__(self): super().__init__( name="exposure_number", sources=["nexpsur"], - # force_unique=True, - # required=True, + force_unique=True, + required=True, ) @@ -957,6 +1008,33 @@ def _init_hook(self, item): # --------------------------------------------- # Mixins to define the broad category of rules. # --------------------------------------------- +class AsnMixin_Lv2FOV: + """Level 2 Image association base""" + + def _init_hook(self, item): + """Post-check and pre-add initialization""" + + super()._init_hook(item) + self.data["asn_type"] = "FOV" + + def finalize(self): + """Finalize association + + + Returns + ------- + associations: [association[, ...]] or None + List of fully-qualified associations that this association + represents. + `None` if a complete association cannot be produced. + + """ + if self.is_valid: + return self.make_fov_asn() + else: + return None + + class AsnMixin_Lv2Image: """Level 2 Image association base""" diff --git a/romancal/associations/lib/rules_level2.py b/romancal/associations/lib/rules_level2.py index 03fb7dbee..80a324909 100644 --- a/romancal/associations/lib/rules_level2.py +++ b/romancal/associations/lib/rules_level2.py @@ -6,7 +6,14 @@ from romancal.associations.lib.rules_elpp_base import * from romancal.associations.registry import RegistryMarker -__all__ = ["Asn_Lv2Image", "Asn_Lv2GBTDSPass", "Asn_Lv2GBTDSFull", "AsnMixin_Lv2Image"] +__all__ = [ + "Asn_Lv2FOV", + "Asn_Lv2Image", + "Asn_Lv2GBTDSPass", + "Asn_Lv2GBTDSFull", + "AsnMixin_Lv2Image", + "AsnMinxin_Lv2FOV", +] # Configure logging logger = logging.getLogger(__name__) @@ -17,16 +24,40 @@ # -------------------------------- # Start of the User-level rules # -------------------------------- +@RegistryMarker.rule +class Asn_Lv2FOV(AsnMixin_Lv2FOV, DMS_ELPP_Base): + """Level2b Non-TSO Science Image Association + + Characteristics: + - Association type: ``FOV`` + - Pipeline: ``mosaic`` + - Image-based science exposures + - Science exposures for all 18 detectors + """ + + def __init__(self, *args, **kwargs): + # Setup constraints + self.constraints = Constraint( + [ + Constraint_Base(), + Constraint_Target(), + Constraint_Filename(), + ] + ) + + # Now check and continue initialization. + super().__init__(*args, **kwargs) + + @RegistryMarker.rule class Asn_Lv2Image(AsnMixin_Lv2Image, DMS_ELPP_Base): """Level2b Non-TSO Science Image Association Characteristics: - - Association type: ``image2`` - - Pipeline: ``calwebb_image2`` + - Association type: ``image`` + - Pipeline: ``ELPP`` - Image-based science exposures - Single science exposure - - Non-TSO """ def __init__(self, *args, **kwargs): @@ -35,7 +66,12 @@ def __init__(self, *args, **kwargs): [ Constraint_Base(), Constraint_Target(), - Constraint_Expos(), + Constraint( + [ + Constraint_Expos(), + ], + reduce=Constraint.any, + ), Constraint_Optical_Path(), Constraint_Sequence(), Constraint_Pass(), diff --git a/romancal/associations/tests/test_level2_basics.py b/romancal/associations/tests/test_level2_basics.py index 125c6f4d4..663452150 100644 --- a/romancal/associations/tests/test_level2_basics.py +++ b/romancal/associations/tests/test_level2_basics.py @@ -59,7 +59,10 @@ def test_level2_productname(): for member in product["members"] if member["exptype"] == "science" or member["exptype"] == "wfi_image" ] - assert len(science) == 2 + if asn["asn_rule"] == "Asn_Lv2Image": + assert len(science) == 2 + if asn["asn_rule"] == "Asn_Lv2FOV": + assert len(science) == 18 # match = re.match(REGEX_LEVEL2, science[0]['expname']) diff --git a/romancal/associations/tests/test_level2_candidates.py b/romancal/associations/tests/test_level2_candidates.py index 6110114ef..bb85000a5 100644 --- a/romancal/associations/tests/test_level2_candidates.py +++ b/romancal/associations/tests/test_level2_candidates.py @@ -17,11 +17,11 @@ # Basic observation ACIDs (["-i", "o001"], 0), # Whole program - ([], 2), + ([], 5), # Discovered only (["--discover"], 0), # Candidates only - (["--all-candidates"], 2), + (["--all-candidates"], 5), ], ) def test_candidate_observation(partial_args, n_asns): From 98db5549760e290ae63d43f58cee015292d63730 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Mon, 2 Oct 2023 10:24:38 -0400 Subject: [PATCH 02/16] 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 879ba9ef7a06ed734c3554ae91827320f4f7e1c5 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Mon, 2 Oct 2023 10:25:49 -0400 Subject: [PATCH 03/16] 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 cffe8bd44797154d9d281e93c3a3cb076b121ec8 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 04/16] [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 59951d43130ab3e66f7e84f4399c94da1fadabbe Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 06:49:24 -0400 Subject: [PATCH 05/16] 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 144e0daccf5f143b9e8e57ef552de67355098c95 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 06/16] [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 6c93ef5721aa375e6405501a4af315d1f0d66fed Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 07:00:02 -0400 Subject: [PATCH 07/16] 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 5cb0b0d1c01650f43131cfd4b18b9503f70f0555 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 09:25:59 -0400 Subject: [PATCH 08/16] 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 881f00cacc151396e9c59d9c082ea60429a42096 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 14:10:22 -0400 Subject: [PATCH 09/16] 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 b6254930cb58c3b47016cd166536af4f2eb963a6 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 5 Oct 2023 14:14:48 -0400 Subject: [PATCH 10/16] 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 ea17d2aeb47c8d34467e11d4c90f690d2737e7a8 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 11/16] [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 f93c56b994d153b3f27c49fed8f078b764a3ee61 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Wed, 11 Oct 2023 19:00:53 -0400 Subject: [PATCH 12/16] 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 9fd104e29cdf5bf015003ccbce3b57c03fc8accf Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Wed, 11 Oct 2023 19:28:35 -0400 Subject: [PATCH 13/16] 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 9d59ab77b755c839450e4dbae36364c91430f882 Mon Sep 17 00:00:00 2001 From: "M. Teodoro" Date: Thu, 12 Oct 2023 11:03:45 -0400 Subject: [PATCH 14/16] 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(): From 525019b6e129522eaa9814c7a46691a9e4908931 Mon Sep 17 00:00:00 2001 From: Zach Burnett Date: Fri, 13 Oct 2023 08:15:28 -0400 Subject: [PATCH 15/16] add instructions for downloading WebbPSF data (#937) --- README.md | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README.md b/README.md index 852037c76..a10982257 100644 --- a/README.md +++ b/README.md @@ -185,6 +185,8 @@ $ crds sync --contexts roman-edit The CRDS_READONLY_CACHE variable should not be set, since references will need to be downloaded to your local cache as they are requested. +Additionally, currently WebbPSF data is also required. Follow [these instructions to download the data files / point to existing files on the shared internal network](https://webbpsf.readthedocs.io/en/latest/installation.html#data-install). + ### Running tests Unit tests can be run via `pytest`. Within the top level of your local `roman` repo checkout: From 21f9e7c72964f791289fc310df62531acd9ff4fd Mon Sep 17 00:00:00 2001 From: Zach Burnett Date: Fri, 13 Oct 2023 09:20:38 -0400 Subject: [PATCH 16/16] add `.git` to dev dependency specifications (#932) --- requirements-dev-st.txt | 16 ++++++++-------- requirements-dev-thirdparty.txt | 12 ++++++------ 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/requirements-dev-st.txt b/requirements-dev-st.txt index f913d6e6e..c26e35579 100644 --- a/requirements-dev-st.txt +++ b/requirements-dev-st.txt @@ -1,13 +1,13 @@ # Roman upstream packages -git+https://github.com/spacetelescope/roman_datamodels -git+https://github.com/spacetelescope/rad +git+https://github.com/spacetelescope/roman_datamodels.git +git+https://github.com/spacetelescope/rad.git # shared upstream packages -git+https://github.com/spacetelescope/stcal -git+https://github.com/spacetelescope/stpipe +git+https://github.com/spacetelescope/stcal.git +git+https://github.com/spacetelescope/stpipe.git # Other important upstream packages -git+https://github.com/spacetelescope/crds -git+https://github.com/spacetelescope/gwcs -git+https://github.com/spacetelescope/metrics_logger -git+https://github.com/spacetelescope/tweakwcs +git+https://github.com/spacetelescope/crds.git +git+https://github.com/spacetelescope/gwcs.git +git+https://github.com/spacetelescope/metrics_logger.git +git+https://github.com/spacetelescope/tweakwcs.git diff --git a/requirements-dev-thirdparty.txt b/requirements-dev-thirdparty.txt index 0568dde4b..7e3c8d34f 100644 --- a/requirements-dev-thirdparty.txt +++ b/requirements-dev-thirdparty.txt @@ -1,12 +1,12 @@ # ASDF upstream packages -git+https://github.com/asdf-format/asdf-standard -git+https://github.com/asdf-format/asdf -git+https://github.com/asdf-format/asdf-transform-schemas -git+https://github.com/asdf-format/asdf-coordinates-schemas -git+https://github.com/asdf-format/asdf-wcs-schemas +git+https://github.com/asdf-format/asdf-standard.git +git+https://github.com/asdf-format/asdf.git +git+https://github.com/asdf-format/asdf-transform-schemas.git +git+https://github.com/asdf-format/asdf-coordinates-schemas.git +git+https://github.com/asdf-format/asdf-wcs-schemas.git # Use weekly astropy dev build -git+https://github.com/astropy/asdf-astropy +git+https://github.com/astropy/asdf-astropy.git --extra-index-url https://pypi.anaconda.org/astropy/simple astropy --pre git+https://github.com/astropy/photutils.git