Coverage for python / lsst / ap / association / mpSkyEphemerisQuery.py: 24%
106 statements
« prev ^ index » next coverage.py v7.13.2, created at 2026-02-10 01:41 -0800
1# This file is part of ap_association.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22"""Solar System Object Query to MPSky in place of a internal Rubin solar
23system object caching/retrieval code.
25Will compute the location for of known SSObjects within a visit-detector combination.
26"""
28__all__ = ["MPSkyEphemerisQueryConfig", "MPSkyEphemerisQueryTask"]
31import numpy as np
32import os
33import pandas as pd
34import pyarrow as pa
35import requests
37from lsst.ap.association.utils import getMidpointFromTimespan
38from lsst.geom import SpherePoint
39import lsst.pex.config as pexConfig
40from lsst.pipe.tasks.associationUtils import obj_id_to_ss_object_id
41from lsst.utils.timer import timeMethod
43from lsst.pipe.base import connectionTypes, NoWorkFound, PipelineTask, \
44 PipelineTaskConfig, PipelineTaskConnections, Struct
46ELEMENTS_LIST = "q,e,inc,node,argPeri,t_p_MJD_TDB,epochMJD_TDB".split(',')
47ELEMENTS_FORMAT = " ".join(["{:> 11f}"] * len(ELEMENTS_LIST))
48NAMING_MAP = {'inc': 'incl', 'argPeri': 'peri', 't_p_MJD_TDB': 't_p', 'epochMJD_TDB': 'epoch',
49 'q': 'q', 'e': 'e', 'node': 'node'}
51ELEMENTS_LIST_MPC_ORBITS = ['designation', 'id', 'packed_primary_provisional_designation',
52 'unpacked_primary_provisional_designation', 'mpc_orb_jsonb',
53 'created_at', 'updated_at', 'orbit_type_int', 'u_param', 'nopp',
54 'arc_length_total', 'arc_length_sel', 'nobs_total', 'nobs_total_sel',
55 'a', 'q', 'e', 'i', 'node', 'argperi', 'peri_time', 'yarkovsky', 'srp',
56 'a1', 'a2', 'a3', 'dt', 'mean_anomaly', 'period', 'mean_motion', 'a_unc',
57 'q_unc', 'e_unc', 'i_unc', 'node_unc', 'argperi_unc', 'peri_time_unc',
58 'yarkovsky_unc', 'srp_unc', 'a1_unc', 'a2_unc', 'a3_unc', 'dt_unc',
59 'mean_anomaly_unc', 'period_unc', 'mean_motion_unc', 'epoch_mjd', 'h', 'g',
60 'not_normalized_rms', 'normalized_rms', 'earth_moid', 'fitting_datetime']
63class MPSkyEphemerisQueryConnections(PipelineTaskConnections,
64 dimensions=("instrument",
65 "group", "detector")):
67 predictedRegionTime = connectionTypes.Input(
68 doc="The predicted exposure region and time",
69 name="regionTimeInfo",
70 storageClass="RegionTimeInfo",
71 dimensions={"instrument", "group", "detector"},
72 )
74 ssObjects = connectionTypes.Output(
75 doc="MPSky-provided Solar System objects observable in this detector-visit",
76 name="preloaded_SsObjects",
77 storageClass="ArrowAstropy",
78 dimensions=("instrument", "group", "detector"),
79 )
82class MPSkyEphemerisQueryConfig(
83 PipelineTaskConfig,
84 pipelineConnections=MPSkyEphemerisQueryConnections):
85 observerCode = pexConfig.Field(
86 dtype=str,
87 doc="IAU Minor Planet Center observer code for queries "
88 "(default is X05 for Rubin Obs./LSST)",
89 default='X05'
90 )
91 queryBufferRadiusDegrees = pexConfig.Field(
92 dtype=float,
93 doc="Buffer radius in degrees added to detector bounding circle for ephemeris "
94 "cone search. Defaults to 10 deg/day * 30 minutes",
95 default=0.208
96 )
97 mpSkyRequestTimeoutSeconds = pexConfig.Field(
98 dtype=float,
99 doc="Time in seconds to wait for mpSky request before failing ",
100 default=10.0
101 )
102 mpSkyFallbackURL = pexConfig.Field(
103 dtype=str,
104 doc="mpSky default URL if MP_SKY_URL environment variable unset",
105 default="http://sdfiana014.sdf.slac.stanford.edu:3666/ephemerides/"
106 )
107 mpcorb_prefix = pexConfig.Field(
108 dtype=str,
109 doc="Column prefix to indicate MPCORB columns for downstream tasks.",
110 default="MPCORB_"
111 )
114class MPSkyEphemerisQueryTask(PipelineTask):
115 """Task to query the MPSky service and retrieve the solar system objects
116 that are observable within the input visit.
117 """
118 ConfigClass = MPSkyEphemerisQueryConfig
119 _DefaultName = "mpSkyEphemerisQuery"
121 @timeMethod
122 def run(self, predictedRegionTime):
123 """Parse the information on the current visit and retrieve the
124 observable solar system objects from MPSky.
126 Parameters
127 ----------
128 predictedRegionTime : `lsst.pipe.base.utils.RegionTimeInfo`
129 Predicted footprint and timespan of the exposure
131 Returns
132 -------
133 result : `lsst.pipe.base.Struct`
134 Results struct with components:
136 - ``ssObjects`` : `pandas.DataFrame`
137 DataFrame containing Solar System Objects near the detector
138 footprint as retrieved by MPSky. The columns are as follows:
140 ``Name``
141 object name (`str`)
142 ``ra``
143 RA in decimal degrees (`float`)
144 ``dec``
145 DEC in decimal degrees (`float`)
146 ``obj_X_poly``, ``obj_Y_poly``, ``obj_Z_poly``
147 Chebyshev coefficients for object path
148 ``obs_X_poly``, ``obs_Y_poly``, ``obs_Z_poly``
149 Chebyshev coefficients for observer path
150 ``t_min``
151 Lower time bound for polynomials
152 ``t_max``
153 Upper time bound for polynomials
154 """
155 # Get detector center and timespan midpoint from predictedRegionTime.
156 region = predictedRegionTime.region
157 timespan = predictedRegionTime.timespan
158 expCenter = SpherePoint(region.getBoundingCircle().getCenter())
159 expRadius = region.getBoundingCircle().getOpeningAngle().asDegrees()
160 expMidPointEPOCH = getMidpointFromTimespan(timespan, allowUnbounded=False).mjd
162 # MPSky service query
163 try:
164 mpSkyURL = os.environ['MP_SKY_URL']
165 except KeyError:
166 self.log.warning("MP_SKY_URL is not defined. Falling back to %s.", self.config.mpSkyFallbackURL)
167 mpSkyURL = self.config.mpSkyFallbackURL
168 mpSkySsObjects, elements = self._mpSkyConeSearch(expCenter, expMidPointEPOCH,
169 expRadius + self.config.queryBufferRadiusDegrees,
170 mpSkyURL)
171 if elements is not None:
172 if 'mpc_orb_jsonb' in elements.columns:
173 for element in ELEMENTS_LIST_MPC_ORBITS:
174 mpSkySsObjects[self.config.mpcorb_prefix + element] = elements[element]
175 else:
176 for element in ELEMENTS_LIST:
177 mpSkySsObjects[self.config.mpcorb_prefix + NAMING_MAP[element]] = elements[element]
178 mpSkySsObjects = mpSkySsObjects.fillna(value=np.nan)
179 # ssObjectId and mpcDesignation are also in the MPCORB schema so let's add them here.
180 return Struct(
181 ssObjects=mpSkySsObjects,
182 )
184 def _mpSkyConeSearch(self, expCenter, epochMJD, queryRadius, mpSkyURL):
185 """Query MPSky ephemeris service for objects near the expected detector position
187 Parameters
188 ----------
189 expCenter : `lsst.geom.SpherePoint`
190 Center of search cone
191 epochMJD : `float`
192 Epoch of cone search, (MJD in UTC).
193 queryRadius : `float`
194 Radius of the cone search in degrees.
195 mpSkyURL : `str`
196 URL to query for MPSky.
198 Returns
199 -------
200 mpSkySsObjects : `pandas.DataFrame`
201 DataFrame with Solar System Object information and RA/DEC position
202 within the visit.
203 elements : `pandas.DataFrame` or `None`
204 Orbital elements of ssObjects (to pass to alerts)
205 """
206 fieldRA = expCenter.getRa().asDegrees()
207 fieldDec = expCenter.getDec().asDegrees()
209 params = {
210 "t": epochMJD,
211 "ra": fieldRA,
212 "dec": fieldDec,
213 "radius": queryRadius,
214 "return_elements": "extended"
215 }
217 response = None
218 try:
219 response = requests.get(mpSkyURL, params=params, timeout=self.config.mpSkyRequestTimeoutSeconds)
220 response.raise_for_status()
221 with pa.input_stream(memoryview(response.content)) as fp:
222 fp.seek(0)
223 p = pa.ipc.read_tensor(fp)
224 op = pa.ipc.read_tensor(fp)
225 with pa.ipc.open_stream(fp) as reader:
226 schema = reader.schema
227 r = next(reader)
229 # construct an elements pandas dataframe
230 if "q" in schema.names:
231 if "mpc_orb_jsonb" in schema.names:
232 colnames = ELEMENTS_LIST_MPC_ORBITS
233 else:
234 colnames = ELEMENTS_LIST
235 cols = {col: r[col].to_numpy(zero_copy_only=False) for col in colnames}
236 elements = pd.DataFrame(cols)
237 for c in cols:
238 if c in ['created_at', 'fitting_datetime', 'updated_at']:
239 elements[c] = None
240 # TODO (DM-53952): Handle properly, probably as a datetime
241 else:
242 elements = None
244 mpSkySsObjects = pd.DataFrame()
245 mpSkySsObjects['ObjID'] = r["name"].to_numpy(zero_copy_only=False)
246 mpSkySsObjects['ra'] = r["ra"].to_numpy()
247 mpSkySsObjects['dec'] = r["dec"].to_numpy()
248 mpSkySsObjects['tmin'] = r["tmin"].to_numpy()
249 mpSkySsObjects['tmax'] = r["tmax"].to_numpy()
251 mpSkySsObjects['trailedSourceMagTrue'] = np.nan
252 if 'Vmag' in r:
253 mpSkySsObjects['trailedSourceMagTrue'] = r["Vmag"].to_numpy()
255 object_polynomial = p.to_numpy()
256 observer_polynomial = op.to_numpy()
258 mpSkySsObjects['obj_x_poly'] = [poly[0] for poly in object_polynomial.T]
259 mpSkySsObjects['obj_y_poly'] = [poly[1] for poly in object_polynomial.T]
260 mpSkySsObjects['obj_z_poly'] = [poly[2] for poly in object_polynomial.T]
261 mpSkySsObjects['obs_x_poly'] = [observer_polynomial.T[0] for
262 i in range(len(mpSkySsObjects))]
263 mpSkySsObjects['obs_y_poly'] = [observer_polynomial.T[1] for
264 i in range(len(mpSkySsObjects))]
265 mpSkySsObjects['obs_z_poly'] = [observer_polynomial.T[2] for
266 i in range(len(mpSkySsObjects))]
267 mpSkySsObjects['Err(arcsec)'] = 2
268 if 'packed_primary_provisional_designation' in elements.columns:
269 desigs = elements['packed_primary_provisional_designation'].values
270 else:
271 desigs = mpSkySsObjects['ObjID'].values
272 mpSkySsObjects['ssObjectId'] = [obj_id_to_ss_object_id(v) for v in desigs]
274 nFound = len(mpSkySsObjects)
276 if nFound == 0:
277 self.log.info("No Solar System objects found for visit.")
278 else:
279 self.log.info("%d Solar System Objects in visit", nFound)
280 except requests.RequestException as e:
281 body = getattr(response, "text", "<No response>")
282 message = f"Query to the remote ephemerides service failed. Got response {body}"
283 self.log.error(message)
284 raise NoWorkFound(f"{message}: {e}") from e
286 return mpSkySsObjects, elements