From 50209266621abe1e437cf9e5f72a9deb831740f3 Mon Sep 17 00:00:00 2001 From: Joachim Moeyens Date: Mon, 6 May 2024 12:19:39 -0700 Subject: [PATCH] Add ADESObservations --- adam_core/observations/__init__.py | 8 + adam_core/observations/ades.py | 291 ++++++++++++++++++ adam_core/observations/tests/test_ades.py | 235 ++++++++++++++ .../observations/tests/testdata/__init__.py | 0 .../tests/testdata/sample_ades.psv | 104 +++++++ 5 files changed, 638 insertions(+) create mode 100644 adam_core/observations/ades.py create mode 100644 adam_core/observations/tests/test_ades.py create mode 100644 adam_core/observations/tests/testdata/__init__.py create mode 100644 adam_core/observations/tests/testdata/sample_ades.psv diff --git a/adam_core/observations/__init__.py b/adam_core/observations/__init__.py index 3e787817..b3b4715d 100644 --- a/adam_core/observations/__init__.py +++ b/adam_core/observations/__init__.py @@ -1,6 +1,14 @@ # flake8: noqa: F401 # TODO: associations - should it be in the API? its useful for test helpers, but kind of niche +from .ades import ( + ADESObservations, + ObsContext, + ObservatoryObsContext, + SoftwareObsContext, + SubmitterObsContext, + TelescopeObsContext, +) from .associations import Associations from .detections import PointSourceDetections from .exposures import Exposures diff --git a/adam_core/observations/ades.py b/adam_core/observations/ades.py new file mode 100644 index 00000000..e16ecbee --- /dev/null +++ b/adam_core/observations/ades.py @@ -0,0 +1,291 @@ +import os +from dataclasses import asdict, dataclass +from typing import Optional + +import numpy as np +import pyarrow.compute as pc +import quivr as qv +from astropy.time import Time + +from ..time import Timestamp + +STRING100 = 100 +STRING25 = 25 + + +@dataclass +class ObservatoryObsContext: + mpcCode: str + name: Optional[str] = None + + def __post_init__(self): + assert len(self.mpcCode) in [3, 4] + if self.name is not None: + assert len(self.name) <= STRING100 + + +@dataclass +class SubmitterObsContext: + name: str + institution: Optional[str] = None + + def __post_init__(self): + assert len(self.name) <= STRING100 + if self.institution is not None: + assert len(self.institution) <= STRING100 + + +@dataclass +class TelescopeObsContext: + name: str + design: str + aperture: Optional[float] = None + detector: Optional[str] = None + fRatio: Optional[float] = None + filter: Optional[str] = None + arraySize: Optional[str] = None + pixelSize: Optional[float] = None + + def __post_init__(self): + assert len(self.name) <= STRING100 + assert len(self.design) <= STRING25 + if self.aperture is not None: + assert self.aperture > 0 + if self.detector is not None: + assert len(self.detector) <= STRING25 + if self.fRatio is not None: + assert self.fRatio > 0 + if self.filter is not None: + assert len(self.filter) <= STRING25 + if self.arraySize is not None: + assert len(self.arraySize) <= STRING25 + if self.pixelSize is not None: + assert self.pixelSize > 0 + + +@dataclass +class SoftwareObsContext: + astrometry: Optional[str] = None + fitOrder: Optional[str] = None + photometry: Optional[str] = None + objectDetection: Optional[str] = None + + def __post_init__(self): + if self.astrometry is not None: + assert len(self.astrometry) <= STRING100 + if self.fitOrder is not None: + assert len(self.fitOrder) <= STRING25 + if self.photometry is not None: + assert len(self.photometry) <= STRING100 + if self.objectDetection is not None: + assert len(self.objectDetection) <= STRING100 + + +@dataclass +class ObsContext: + observatory: ObservatoryObsContext + submitter: SubmitterObsContext + observers: list[str] + measurers: list[str] + telescope: TelescopeObsContext + software: Optional[SoftwareObsContext] = None + coinvestigators: Optional[list[str]] = None + collaborators: Optional[list[str]] = None + fundingSource: Optional[str] = None + comments: Optional[list[str]] = None + + def __post_init__(self): + assert len(self.observers) > 0 + for observer in self.observers: + assert len(observer) <= STRING100 + assert len(self.measurers) > 0 + for measurer in self.measurers: + assert len(measurer) <= STRING100 + if self.coinvestigators is not None: + assert len(self.coinvestigators) > 0 + for coinvestigator in self.coinvestigators: + assert len(coinvestigator) <= STRING100 + if self.collaborators is not None: + assert len(self.collaborators) > 0 + for collaborator in self.collaborators: + assert len(collaborator) <= STRING100 + if self.fundingSource is not None: + assert len(self.fundingSource) <= STRING100 + if self.comments is not None: + assert len(self.comments) > 0 + for comment in self.comments: + assert len(comment) <= STRING100 + + def to_string(self) -> str: + lines = [] + for k, v in asdict(self).items(): + if isinstance(v, dict): + lines.append(f"# {k}") + for k2, v2 in v.items(): + if v2 is not None: + lines.append(f"! {k2} {v2}") + else: + if v is not None: + if k in [ + "observers", + "measurers", + "coinvestigators", + "collaborators", + ]: + lines.append(f"# {k}") + for name in v: + lines.append(f"! name {name}") + elif k == "fundingSource": + lines.append(f"# fundingSource {v}") + elif k == "comments": + lines.append("# comment") + for comment in v: + lines.append(f"! line {comment}") + return "\n".join(lines) + "\n" + + +class ADESObservations(qv.Table): + + permID = qv.LargeStringColumn(nullable=True) + provID = qv.LargeStringColumn(nullable=True) + trkSub = qv.LargeStringColumn(nullable=True) + obsSubID = qv.LargeStringColumn(nullable=True) + obsTime = Timestamp.as_column() + ra = qv.Float64Column() + dec = qv.Float64Column() + rmsRA = qv.Float64Column(nullable=True) + rmsDec = qv.Float64Column(nullable=True) + mag = qv.Float64Column(nullable=True) + rmsMag = qv.Float64Column(nullable=True) + band = qv.LargeStringColumn(nullable=True) + stn = qv.LargeStringColumn() + mode = qv.LargeStringColumn() + astCat = qv.LargeStringColumn() + remarks = qv.LargeStringColumn(nullable=True) + + def to_psv( + self, + file_out: str, + obs_contexts: dict[str, ObsContext], + seconds_precision: int = 3, + columns_precision: dict[str, int] = { + "ra": 8, + "dec": 8, + "rmsRA": 4, + "rmsDec": 4, + "mag": 2, + "rmsMag": 2, + }, + ): + """ + Save observations to a MPC-submittable ADES psv file. + + Parameters + ---------- + file_out : str + The file to save the observations to. + obs_contexts : dict[str, ObsContext] + A dictionary of observatory codes and their corresponding ObsContexts to use + as the context headers for the different observatory codes in the observations. + seconds_precision : int, optional + The precision to use for the seconds in the obsTime field, by default 3. + columns_precision : dict[str, int], optional + A dictionary of column names and their corresponding precision to use when writing + the observations to the file, by default { + "ra": 8, + "dec": 8, + "rmsRA" : 4, + "rmsDec" : 4, + "mag": 2, + "rmsMag": 2, + } + The MPC enforces strict limits on these and submitters may need permission to send + high-precision data. + """ + if os.path.exists(file_out): + raise FileExistsError(f"{file_out} already exists.") + + with open(file_out, "a") as f: + + unique_observatories = self.stn.unique().to_numpy(zero_copy_only=False) + unique_observatories.sort() + + f.writelines(["# version=2022\n"]) + + for obs in unique_observatories: + if obs not in obs_contexts: + raise ValueError(f"Observatory {obs} not found in obs_contexts") + + observations_obscode = self.select("stn", obs) + observations_obscode = observations_obscode.sort_by( + [ + ("provID", "ascending"), + ("permID", "ascending"), + ("trkSub", "ascending"), + ("obsTime.days", "ascending"), + ("obsTime.nanos", "ascending"), + ] + ) + + id_present = False + if not pc.all(pc.is_null(observations_obscode.permID)).as_py(): + id_present = True + if not pc.all(pc.is_null(observations_obscode.provID)).as_py(): + id_present = True + if not pc.all(pc.is_null(observations_obscode.trkSub)).as_py(): + id_present = True + + if not id_present: + err = ( + "At least one of permID, provID, or trkSub should\n" + "be present in observations." + ) + raise ValueError(err) + + # Write the observatory context block + obs_context = obs_contexts[obs] + f.writelines([obs_context.to_string()]) + + # Write the observations block + ades = observations_obscode.to_dataframe() + + # Convert the timestamp to ISOT with the desired precision + observation_times = Time( + observations_obscode.obsTime.rescale("utc") + .mjd() + .to_numpy(zero_copy_only=False), + format="mjd", + precision=seconds_precision, + ) + ades.insert( + 4, + "obsTime", + np.array([i + "Z" for i in observation_times.utc.isot]), + ) + ades.drop(columns=["obsTime.days", "obsTime.nanos"], inplace=True) + + # Multiply rmsRA by cos(dec) since ADES wants the random component in rmsRAcosDec + ades.loc[:, "rmsRA"] *= np.cos(np.radians(ades["dec"])) + + # Convert rmsRA and rmsDec to arcseconds + ades.loc[:, "rmsRA"] *= 3600 + ades.loc[:, "rmsDec"] *= 3600 + + ades.dropna(how="all", axis=1, inplace=True) + + # Change the precision of some of the columns to conform + # to MPC standards + for col, prec_col in columns_precision.items(): + print(col) + if col in ades.columns: + ades[col] = [ + f"{i:.{prec_col}f}" + if i is not None or not np.isnan(i) + else "" + for i in ades[col] + ] + + # This will append to the file since the file is opened in append mode + ades.to_csv(f, sep="|", header=True, index=False, float_format="%.16f") + + return diff --git a/adam_core/observations/tests/test_ades.py b/adam_core/observations/tests/test_ades.py new file mode 100644 index 00000000..b977bb01 --- /dev/null +++ b/adam_core/observations/tests/test_ades.py @@ -0,0 +1,235 @@ +from importlib.resources import files + +import pytest + +from ...time import Timestamp +from ..ades import ( + ADESObservations, + ObsContext, + ObservatoryObsContext, + SoftwareObsContext, + SubmitterObsContext, + TelescopeObsContext, +) + + +@pytest.fixture +def ades_obscontext(): + # Nearly real metadata used for ADAM::THOR observations + measurers = [ + "J. Moeyens", + "M. Juric", + "S. Nelson", + "A. Koumjian", + "K. Kiker", + "N. Tellis", + "D. Veronese-Milin", + "A. Posner", + "E. Lu", + "C. Fiaschetti", + "D. Remy", + ] + + software_context = SoftwareObsContext( + objectDetection="Asteroid Discovery Analysis and Mapping + Tracklet-less Heliocentric Orbit Recovery (ADAM::THOR)", + ) + + submitter_context = SubmitterObsContext( + name="J. Moeyens", institution="B612 Asteroid Institute" + ) + + fundingSource = "WilliamBowes, McGregorGirand, Tito's Vodka, PRawls, SKraus, Yishan/KWong, SGalitsky, Google" + comments = [ + "THIS IS A TEST FILE CONTAINING FAKE OBSERVATIONS", + "Discovery candidates found by members of the Asteroid Institute, a program of B612 Foundation,", + "using the institute's ADAM::THOR discovery service running on Google Cloud.", + ] + + obs_contexts = { + "W84": ObsContext( + observatory=ObservatoryObsContext( + mpcCode="W84", name="Cerro Tololo - Blanco + DECam" + ), + submitter=submitter_context, + observers=["D. E. Survey"], + measurers=measurers, + telescope=TelescopeObsContext( + name="Blanco 4m", design="Reflector", detector="CCD", aperture=4.0 + ), + software=software_context, + fundingSource=fundingSource, + comments=comments, + ), + "695": ObsContext( + observatory=ObservatoryObsContext( + mpcCode="695", name="Kitt Peak National Observatory - Mayall + Mosaic3" + ), + submitter=submitter_context, + observers=["M. L. Survey"], + measurers=measurers, + telescope=TelescopeObsContext( + name="Mayall 4m", + design="Reflector", + detector="CCD", + aperture=4.0, + ), + software=software_context, + fundingSource=fundingSource, + comments=comments, + ), + "V00": ObsContext( + observatory=ObservatoryObsContext( + mpcCode="V00", name="Kitt Peak National Observatory - Bok + 90Prime" + ), + submitter=submitter_context, + observers=["B. A. S. Survey"], + measurers=measurers, + telescope=TelescopeObsContext( + name="Bok 2.3m", + design="Reflector", + detector="CCD", + aperture=2.3, + ), + software=software_context, + fundingSource=fundingSource, + comments=comments, + ), + } + + return obs_contexts + + +@pytest.fixture +def ades_observations(): + + observations = ADESObservations.from_kwargs( + permID=["3000", "3000", "3001", "3001"], + trkSub=["a1234b", "a1234b", "a2345b", "a2345b"], + obsSubID=["obs01", "obs02", "obs03", "obs04"], + obsTime=Timestamp.from_mjd( + [60434.0, 60434.1, 60435.0, 60435.2], + scale="utc", + ), + ra=[240.00, 240.05, 15.00, 15.05], + dec=[-15.00, -15.05, 10.00, 10.05], + rmsRA=[1 / 3600, 1 / 3600, None, None], + rmsDec=[1 / 3600, 1 / 3600, None, None], + mag=[20.0, 20.3, None, 21.4], + band=["r", "g", None, "r"], + stn=["W84", "W84", "V00", "695"], + mode=["CCD", "CCD", "CCD", "CCD"], + astCat=["Gaia2", "Gaia2", "Gaia2", "Gaia2"], + remarks=[ + "This is a dummy observation", + "This is another dummy observation", + None, + "This is the fourth dummy observation", + ], + ) + return observations + + +def test_ObsContext_to_string(ades_obscontext): + # Test that we can convert an ObsContext to a string representation + W84 = ades_obscontext["W84"] + string = W84.to_string() + + assert ( + string + == """# observatory +! mpcCode W84 +! name Cerro Tololo - Blanco + DECam +# submitter +! name J. Moeyens +! institution B612 Asteroid Institute +# observers +! name D. E. Survey +# measurers +! name J. Moeyens +! name M. Juric +! name S. Nelson +! name A. Koumjian +! name K. Kiker +! name N. Tellis +! name D. Veronese-Milin +! name A. Posner +! name E. Lu +! name C. Fiaschetti +! name D. Remy +# telescope +! name Blanco 4m +! design Reflector +! aperture 4.0 +! detector CCD +# software +! objectDetection Asteroid Discovery Analysis and Mapping + Tracklet-less Heliocentric Orbit Recovery (ADAM::THOR) +# fundingSource WilliamBowes, McGregorGirand, Tito's Vodka, PRawls, SKraus, Yishan/KWong, SGalitsky, Google +# comment +! line THIS IS A TEST FILE CONTAINING FAKE OBSERVATIONS +! line Discovery candidates found by members of the Asteroid Institute, a program of B612 Foundation, +! line using the institute's ADAM::THOR discovery service running on Google Cloud. +""" + ) + + V00 = ades_obscontext["V00"] + string = V00.to_string() + + assert ( + string + == """# observatory +! mpcCode V00 +! name Kitt Peak National Observatory - Bok + 90Prime +# submitter +! name J. Moeyens +! institution B612 Asteroid Institute +# observers +! name B. A. S. Survey +# measurers +! name J. Moeyens +! name M. Juric +! name S. Nelson +! name A. Koumjian +! name K. Kiker +! name N. Tellis +! name D. Veronese-Milin +! name A. Posner +! name E. Lu +! name C. Fiaschetti +! name D. Remy +# telescope +! name Bok 2.3m +! design Reflector +! aperture 2.3 +! detector CCD +# software +! objectDetection Asteroid Discovery Analysis and Mapping + Tracklet-less Heliocentric Orbit Recovery (ADAM::THOR) +# fundingSource WilliamBowes, McGregorGirand, Tito's Vodka, PRawls, SKraus, Yishan/KWong, SGalitsky, Google +# comment +! line THIS IS A TEST FILE CONTAINING FAKE OBSERVATIONS +! line Discovery candidates found by members of the Asteroid Institute, a program of B612 Foundation, +! line using the institute's ADAM::THOR discovery service running on Google Cloud. +""" + ) + + +@pytest.fixture +def ades_psv(ades_obscontext, ades_observations): + actual = files("adam_core.observations.tests.testdata").joinpath( + "sample_ades_actual.psv" + ) + ades_observations.to_psv(actual, ades_obscontext) + yield actual + actual.unlink() + + +def test_ADESObservations_to_psv(ades_psv): + # Test that we can convert ADESObservations to a PSV file + desired = files("adam_core.observations.tests.testdata").joinpath("sample_ades.psv") + + with open(desired, "r") as f: + desired_lines = f.readlines() + + with open(ades_psv, "r") as f: + actual_lines = f.readlines() + + assert desired_lines == actual_lines diff --git a/adam_core/observations/tests/testdata/__init__.py b/adam_core/observations/tests/testdata/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/adam_core/observations/tests/testdata/sample_ades.psv b/adam_core/observations/tests/testdata/sample_ades.psv new file mode 100644 index 00000000..573e3bf0 --- /dev/null +++ b/adam_core/observations/tests/testdata/sample_ades.psv @@ -0,0 +1,104 @@ +# version=2022 +# observatory +! mpcCode 695 +! name Kitt Peak National Observatory - Mayall + Mosaic3 +# submitter +! name J. Moeyens +! institution B612 Asteroid Institute +# observers +! name M. L. Survey +# measurers +! name J. Moeyens +! name M. Juric +! name S. Nelson +! name A. Koumjian +! name K. Kiker +! name N. Tellis +! name D. Veronese-Milin +! name A. Posner +! name E. Lu +! name C. Fiaschetti +! name D. Remy +# telescope +! name Mayall 4m +! design Reflector +! aperture 4.0 +! detector CCD +# software +! objectDetection Asteroid Discovery Analysis and Mapping + Tracklet-less Heliocentric Orbit Recovery (ADAM::THOR) +# fundingSource WilliamBowes, McGregorGirand, Tito's Vodka, PRawls, SKraus, Yishan/KWong, SGalitsky, Google +# comment +! line THIS IS A TEST FILE CONTAINING FAKE OBSERVATIONS +! line Discovery candidates found by members of the Asteroid Institute, a program of B612 Foundation, +! line using the institute's ADAM::THOR discovery service running on Google Cloud. +permID|trkSub|obsSubID|obsTime|ra|dec|mag|band|stn|mode|astCat|remarks +3001|a2345b|obs04|2024-05-05T04:48:00.000Z|15.05000000|10.05000000|21.40|r|695|CCD|Gaia2|This is the fourth dummy observation +# observatory +! mpcCode V00 +! name Kitt Peak National Observatory - Bok + 90Prime +# submitter +! name J. Moeyens +! institution B612 Asteroid Institute +# observers +! name B. A. S. Survey +# measurers +! name J. Moeyens +! name M. Juric +! name S. Nelson +! name A. Koumjian +! name K. Kiker +! name N. Tellis +! name D. Veronese-Milin +! name A. Posner +! name E. Lu +! name C. Fiaschetti +! name D. Remy +# telescope +! name Bok 2.3m +! design Reflector +! aperture 2.3 +! detector CCD +# software +! objectDetection Asteroid Discovery Analysis and Mapping + Tracklet-less Heliocentric Orbit Recovery (ADAM::THOR) +# fundingSource WilliamBowes, McGregorGirand, Tito's Vodka, PRawls, SKraus, Yishan/KWong, SGalitsky, Google +# comment +! line THIS IS A TEST FILE CONTAINING FAKE OBSERVATIONS +! line Discovery candidates found by members of the Asteroid Institute, a program of B612 Foundation, +! line using the institute's ADAM::THOR discovery service running on Google Cloud. +permID|trkSub|obsSubID|obsTime|ra|dec|stn|mode|astCat +3001|a2345b|obs03|2024-05-05T00:00:00.000Z|15.00000000|10.00000000|V00|CCD|Gaia2 +# observatory +! mpcCode W84 +! name Cerro Tololo - Blanco + DECam +# submitter +! name J. Moeyens +! institution B612 Asteroid Institute +# observers +! name D. E. Survey +# measurers +! name J. Moeyens +! name M. Juric +! name S. Nelson +! name A. Koumjian +! name K. Kiker +! name N. Tellis +! name D. Veronese-Milin +! name A. Posner +! name E. Lu +! name C. Fiaschetti +! name D. Remy +# telescope +! name Blanco 4m +! design Reflector +! aperture 4.0 +! detector CCD +# software +! objectDetection Asteroid Discovery Analysis and Mapping + Tracklet-less Heliocentric Orbit Recovery (ADAM::THOR) +# fundingSource WilliamBowes, McGregorGirand, Tito's Vodka, PRawls, SKraus, Yishan/KWong, SGalitsky, Google +# comment +! line THIS IS A TEST FILE CONTAINING FAKE OBSERVATIONS +! line Discovery candidates found by members of the Asteroid Institute, a program of B612 Foundation, +! line using the institute's ADAM::THOR discovery service running on Google Cloud. +permID|trkSub|obsSubID|obsTime|ra|dec|rmsRA|rmsDec|mag|band|stn|mode|astCat|remarks +3000|a1234b|obs01|2024-05-04T00:00:00.000Z|240.00000000|-15.00000000|0.9659|1.0000|20.00|r|W84|CCD|Gaia2|This is a dummy observation +3000|a1234b|obs02|2024-05-04T02:24:00.000Z|240.05000000|-15.05000000|0.9657|1.0000|20.30|g|W84|CCD|Gaia2|This is another dummy observation