DM-42575: switched solar system association to MPSky

python/lsst/ap/association/__init__.py

@@ -28,3 +28,4 @@
 from .transformDiaSourceCatalog import *
 from .utils import *
 from .version import *
+from .mpSkyEphemerisQuery import *

python/lsst/ap/association/dailySolarSystemPrecompute.disabled

@@ -0,0 +1,146 @@
+# This file is part of ap_association.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+
+import os
+import sys
+from collections import defaultdict
+
+from astropy.time import Time
+import pickle
+import pandas as pd
+import healpy as hp
+import numpy as np
+
+from sorcha.ephemeris.simulation_parsing import (Observatory, mjd_tai_to_epoch)
+from sorcha.ephemeris.simulation_setup import create_assist_ephemeris, furnish_spiceypy, generate_simulations, precompute_pointing_information
+from sorcha.ephemeris.simulation_driver import EphemerisGeometryParameters, calculate_rates_and_geometry
+from sorcha.ephemeris.simulation_geometry import integrate_light_time
+from sorcha.utilities.sorchaArguments import sorchaArguments
+from sorcha.modules.PPConfigParser import PPConfigFileParser
+from sorcha.modules.PPReadPointingDatabase import PPReadPointingDatabase
+
+
+import lsst.pex.config
+import lsst.pipe.base 
+
+class DailySolarSystemPrecomputeConfig(lsst.pex.config.Config):
+ pass
+
+class DailySolarSystemPrecomputeTask(lsst.pipe.base.Task):
+ _DefaultName = "dailySolarSystemPrecompute"
+ ConfigClass = DailySolarSystemPrecomputeConfig
+
+ def __init__(self, **kwargs):
+ super().__init__(**kwargs)
+ #Make the database connection
+ self.mpcorb_db = None
+
+ def run(self, date):
+ orbits = self._load_mpcorb()
+ state_vectors = self._compute_positions(orbits, date)
+ self._write_apdb(state_vectors)
+
+ def _compute_positions(self, orbits, date):
+ """Summary
+
+ Parameters
+ ----------
+ orbits : `pandas.DataFrame`
+ MPCORB 6-parameter orbits, plus epoch and ObjID.
+ date : `astropy.time.Time`
+ Time at which to compute states.
+ """
+ chunksize = 100000
+
+ configs = {"ar_obs_code":"X05", "ar_ang_fov": 0, "ar_fov_buffer": 0, 
+ "ar_healpix_order": 0, "t": date.mjd, "outfile": ""}
+
+ cmd_args_dict = {
+ "paramsinput": None,
+ "orbinfile": "",
+ "oifoutput": None,
+ "configfile": None,
+ "outpath": ".",
+ "surveyname": "LSST",
+ "outfilestem": "precompute_practice",
+ "verbose": False,
+ "pointing_database": None,
+ "makeTemporaryEphemerisDatabase": False,
+ "readTemporaryEphemerisDatabase": False,
+ "deleteTemporaryEphemerisDatabase": False
+ }
+
+ args = sorchaArguments(cmd_args_dict)
+
+ pointings = pd.DataFrame()
+ pointings["observationStartMJD_TAI"] = [date.mjd]
+ pointings["visitTime"] = 30
+ pointings["fieldRA"] = 0
+ pointings["fieldDec"] = 0
+ pointings["FieldID"] = 0
+ pointings["observationMidpointMJD_TAI"] = pointings["observationStartMJD_TAI"] + (
+ (pointings["visitTime"] / 2.0) / 86400.0) 
+ pointings = precompute_pointing_information(pointings, args, configs)
+
+
+ ephem, gm_sun, gm_total = create_assist_ephemeris(args)
+ furnish_spiceypy(args)
+ sim_dict = generate_simulations(ephem, gm_sun, gm_total, orbits, args)
+ observatories = Observatory(args)
+
+ nside = 128
+ states = []
+ names = []
+ radecs = []
+
+ for _, pointing in pointings_df.iterrows():
+ mjd_tai = float(pointing["observationMidpointMJD_TAI"])
+ desigs = orbits['ObjID']
+ for obj_id in sorted(desigs):
+ ephem_geom_params = EphemerisGeometryParameters()
+ ephem_geom_params.obj_id = obj_id
+ ephem_geom_params.mjd_tai = mjd_tai
+ v = sim_dict[obj_id]
+ sim, ex = v["sim"], v["ex"]
+ simt = sim.t
+ (
+ ephem_geom_params.rho,
+ ephem_geom_params.rho_mag,
+ _,
+ ephem_geom_params.r_ast,
+ ephem_geom_params.v_ast,
+ ) = integrate_light_time(
+ sim, ex, pointing["JD_TDB"] - ephem.jd_ref, pointing["r_obs"], iter = 1, lt0=0
+ )
+ rho_hat = ephem_geom_params.rho / ephem_geom_params.rho_mag
+
+ state_vec = equatorial_to_ecliptic(np.hstack([ephem_geom_params.r_ast, ephem_geom_params.v_ast]))
+ states.append(state_vec)
+ names.append(obj_id)
+ radecs.append(hp.vec2ang(rho_hat, lonlat=True))
+
+ results = pd.DataFrame()
+ results['ObjID'] = names
+ results[['x', 'y', 'z', 'vx', 'vy', 'vz']] = states 
+ results[['ra', 'dec']] = radecs
+ return results
+

python/lsst/ap/association/mpSkyEphemerisQuery.py

@@ -0,0 +1,225 @@
+# This file is part of ap_association.
+#
+# Developed for the LSST Data Management System.
+# This product includes software developed by the LSST Project
+# (https://www.lsst.org).
+# See the COPYRIGHT file at the top-level directory of this distribution
+# for details of code ownership.
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see <https://www.gnu.org/licenses/>.
+
+"""Solar System Object Query to MPSky in place of a internal Rubin solar
+system object caching/retrieval code.
+
+Will compute the location for of known SSObjects within a visit. This code
+blocks on web requests, so should not be used as part of any real-time or
+time-sensitive system. Use in a larger pipeline at your own risk.
+"""
+
+__all__ = ["MPSkyEphemerisQueryConfig", "MPSkyEphemerisQueryTask"]
+
+
+import pandas as pd
+import numpy as np
+import pyarrow as pa
+import requests
+import sys
+
+import lsst.pipe.base as pipeBase
+import lsst.pex.config as pexConfig
+from lsst.utils.timer import timeMethod
+from lsst.geom import SpherePoint
+
+from lsst.pipe.base import PipelineTask, PipelineTaskConfig, PipelineTaskConnections
+import lsst.pipe.base.connectionTypes as connTypes
+
+# Enforce an error for unsafe column/array value setting in pandas.
+pd.options.mode.chained_assignment = 'raise'
+
+
+class MPSkyEphemerisQueryConnections(PipelineTaskConnections,
+ dimensions=("instrument",
+ "group", "detector")):
+
+ predictedRegionTime = pipeBase.connectionTypes.Input(
+ doc="The predicted exposure region and time",
+ name="predictedRegionTime",
+ storageClass="RegionTimeInfo",
+ dimensions={"instrument", "group", "detector"},
+ )
+
+ ssObjects = connTypes.Output(
+ doc="MPSky-provided Solar System objects observable in this detector-visit",
+ name="preloaded_ssObjects",
+ storageClass="DataFrame",
+ dimensions=("instrument", "group", "detector"),
+ )
+
+
+class MPSkyEphemerisQueryConfig(
+ PipelineTaskConfig,
+ pipelineConnections=MPSkyEphemerisQueryConnections):
+ observerCode = pexConfig.Field(
+ dtype=str,
+ doc="IAU Minor Planet Center observer code for queries "
+ "(Rubin Obs./LSST default is X05)",
+ default='X05'
+ )
+ queryRadiusDegrees = pexConfig.Field(
+ dtype=float,
+ doc="On sky radius for Ephemeris cone search. Defaults "
+ "to the radius of Rubin Obs FoV in degrees",
+ default=1.75
+ )
+ MPSkyURL = pexConfig.Field(
+ dtype=str,
+ doc="URL to query mpsky service",
+ default="http://sdfrome001.sdf.slac.stanford.edu:3666"
+ )
+
+
+class MPSkyEphemerisQueryTask(PipelineTask):
+ """Tasks to query the MPSky service and retrieve the solar system objects
+ that are observable within the input visit.
+ """
+ ConfigClass = MPSkyEphemerisQueryConfig
+ _DefaultName = "MPSkyEphemerisQuery"
+
+ def runQuantum(self, butlerQC, inputRefs, outputRefs):
+ inputs = butlerQC.get(inputRefs)
+ inputs["visit"] = butlerQC.quantum.dataId["visit"]
+
+ outputs = self.run(**inputs)
+
+ butlerQC.put(outputs, outputRefs)
+
+ @timeMethod
+ def run(self, predictedRegionTime):
+ """Parse the information on the current visit and retrieve the
+ observable solar system objects from MPSky.
+
+ Parameters
+ ----------
+ predictedRegionTime : `pipe.base.utils.RegionTimeInfo`
+ RegionTime of the predicted exposure
+
+ Returns
+ -------
+ result : `lsst.pipe.base.Struct`
+ Results struct with components:
+
+ - ``ssObjects``: `pandas.DataFrame`
+ DataFrame containing Solar System Objects in field of view as
+ retrieved by MPSky. The columns are as follows:
+ ``Name``
+ object name (`str`)
+ ``ra``
+ RA in decimal degrees (`float`)
+ ``dec``
+ DEC in decimal degrees (`float`)
+ ``obj_poly``
+ DO NOT USE until t_min issue is resolved
+ ``obs_poly``
+ DO NOT USE until t_min issue is resolved
+
+
+
+ """
+ # Grab the visitInfo from the raw to get the information needed on the
+ # full visit.
+
+ # Midpoint time of the exposure in JD
+ region = predictedRegionTime.region
+ timespan = predictedRegionTime.timespan
+ expCenter = SpherePoint(region.getBoundingCircle().getCenter())
+
+ # Make sure date is non-NaN.
+ expMidPointEPOCH = (timespan.begin.mjd + timespan.end.mjd)/2
+
+ # MPSky service query
+ MPSkySsObjects = self._MPSkyConeSearch(expCenter, expMidPointEPOCH, self.config.queryRadiusDegrees)
+ # Add the visit as an extra column.
+
+ return pipeBase.Struct(
+ ssObjects=MPSkySsObjects,
+ )
+
+ def _MPSkyConeSearch(self, expCenter, epochMJD, queryRadius):
+ """Query MPSky ephemeris service using the exposure boresight.
+
+ Parameters
+ ----------
+ expCenter : `lsst.geom.SpherePoint`
+ Center of Exposure RADEC [deg]
+ epochJD : `float`
+ Mid point JD of exposure, in UTC [EPOCH].
+ queryRadius : `float`
+ Radius of the cone search in degrees.
+
+ Returns
+ -------
+ MPSkySsObjects : `pandas.DataFrame`
+ DataFrame with Solar System Object information and RA/DEC position
+ within the visit.
+ """
+
+ fieldRA = expCenter.getRa().asDegrees()
+ fieldDec = expCenter.getDec().asDegrees()
+
+ params = {
+ "t": epochMJD,
+ "ra": fieldRA,
+ "dec": fieldDec,
+ "radius": queryRadius
+ }
+
+ try:
+ response = requests.get(self.config.MPSkyURL, params=params)
+ response.raise_for_status()
+ except requests.exceptions.ConnectionError as e:
+ print("failed to connect to the remote ephemerides service. details:", file=sys.stderr)
+ print(e, file=sys.stderr)
+ exit(-1)
+
+ with pa.input_stream(memoryview(response.content)) as fp:
+ fp.seek(0)
+ p = pa.ipc.read_tensor(fp)
+ op = pa.ipc.read_tensor(fp)
+ op, p = op, p
+ with pa.ipc.open_stream(fp) as reader:
+ r = next(reader)
+
+ ObjID = r["name"].to_numpy(zero_copy_only=False)
+ ra = r["ra"].to_numpy()
+ dec = r["dec"].to_numpy()
+
+ # ObjID, ra, dec, obj_poly, obs_poly = mpsky.query_service('https://sky.dirac.dev/ephemerides/',
+ # epochMJD, fieldRA, fieldDec, queryRadius)
+
+ MPSkySsObjects = pd.DataFrame()
+ MPSkySsObjects['ObjID'] = ObjID
+ MPSkySsObjects['ra'] = ra
+ MPSkySsObjects['dec'] = dec
+ MPSkySsObjects['obj_poly'] = list(np.zeros((len(MPSkySsObjects), 5))) # fix, eventually
+ MPSkySsObjects['obs_poly'] = list(np.zeros((len(MPSkySsObjects), 5))) # fix, eventually
+ MPSkySsObjects['Err(arcsec)'] = 2
+ MPSkySsObjects['ssObjectId'] = [abs(hash(v)) for v in MPSkySsObjects['ObjID'].values]
+ nFound = len(MPSkySsObjects)
+
+ if nFound == 0:
+ self.log.info("No Solar System objects found for visit.")
+
+ self.log.info("%d Solar System Objects in visit", nFound)
+
+ return MPSkySsObjects