Coverage for python / lsst / ap / association / mpSkyEphemerisQuery.py: 31%
90 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-07 01:56 -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 os
32import pandas as pd
33import pyarrow as pa
34import requests
36from lsst.ap.association.utils import getMidpointFromTimespan
37from lsst.geom import SpherePoint
38import lsst.pex.config as pexConfig
39from lsst.pipe.tasks.associationUtils import obj_id_to_ss_object_id
40from lsst.utils.timer import timeMethod
42from lsst.pipe.base import connectionTypes, NoWorkFound, PipelineTask, \
43 PipelineTaskConfig, PipelineTaskConnections, Struct
45ELEMENTS_LIST = "q,e,inc,node,argPeri,t_p_MJD_TDB,epochMJD_TDB".split(',')
46ELEMENTS_FORMAT = " ".join(["{:> 11f}"] * len(ELEMENTS_LIST))
47NAMING_MAP = {'inc': 'incl', 'argPeri': 'peri', 't_p_MJD_TDB': 't_p', 'epochMJD_TDB': 'epoch',
48 'q': 'q', 'e': 'e', 'node': 'node'}
51class MPSkyEphemerisQueryConnections(PipelineTaskConnections,
52 dimensions=("instrument",
53 "group", "detector")):
55 predictedRegionTime = connectionTypes.Input(
56 doc="The predicted exposure region and time",
57 name="regionTimeInfo",
58 storageClass="RegionTimeInfo",
59 dimensions={"instrument", "group", "detector"},
60 )
62 ssObjects = connectionTypes.Output(
63 doc="MPSky-provided Solar System objects observable in this detector-visit",
64 name="preloaded_SsObjects",
65 storageClass="ArrowAstropy",
66 dimensions=("instrument", "group", "detector"),
67 )
70class MPSkyEphemerisQueryConfig(
71 PipelineTaskConfig,
72 pipelineConnections=MPSkyEphemerisQueryConnections):
73 observerCode = pexConfig.Field(
74 dtype=str,
75 doc="IAU Minor Planet Center observer code for queries "
76 "(default is X05 for Rubin Obs./LSST)",
77 default='X05'
78 )
79 queryBufferRadiusDegrees = pexConfig.Field(
80 dtype=float,
81 doc="Buffer radius in degrees added to detector bounding circle for ephemeris "
82 "cone search. Defaults to 10 deg/day * 30 minutes",
83 default=0.208
84 )
85 mpSkyRequestTimeoutSeconds = pexConfig.Field(
86 dtype=float,
87 doc="Time in seconds to wait for mpSky request before failing ",
88 default=1.0
89 )
90 mpSkyFallbackURL = pexConfig.Field(
91 dtype=str,
92 doc="mpSky default URL if MP_SKY_URL environment variable unset",
93 default="http://sdfiana014.sdf.slac.stanford.edu:3666/ephemerides/"
94 )
95 mpcorb_prefix = pexConfig.Field(
96 dtype=str,
97 doc="Column prefix to indicate MPCORB columns for downstream tasks.",
98 default="MPCORB_"
99 )
102class MPSkyEphemerisQueryTask(PipelineTask):
103 """Task to query the MPSky service and retrieve the solar system objects
104 that are observable within the input visit.
105 """
106 ConfigClass = MPSkyEphemerisQueryConfig
107 _DefaultName = "mpSkyEphemerisQuery"
109 @timeMethod
110 def run(self, predictedRegionTime):
111 """Parse the information on the current visit and retrieve the
112 observable solar system objects from MPSky.
114 Parameters
115 ----------
116 predictedRegionTime : `lsst.pipe.base.utils.RegionTimeInfo`
117 Predicted footprint and timespan of the exposure
119 Returns
120 -------
121 result : `lsst.pipe.base.Struct`
122 Results struct with components:
124 - ``ssObjects`` : `pandas.DataFrame`
125 DataFrame containing Solar System Objects near the detector
126 footprint as retrieved by MPSky. The columns are as follows:
128 ``Name``
129 object name (`str`)
130 ``ra``
131 RA in decimal degrees (`float`)
132 ``dec``
133 DEC in decimal degrees (`float`)
134 ``obj_X_poly``, ``obj_Y_poly``, ``obj_Z_poly``
135 Chebyshev coefficients for object path
136 ``obs_X_poly``, ``obs_Y_poly``, ``obs_Z_poly``
137 Chebyshev coefficients for observer path
138 ``t_min``
139 Lower time bound for polynomials
140 ``t_max``
141 Upper time bound for polynomials
142 """
143 # Get detector center and timespan midpoint from predictedRegionTime.
144 region = predictedRegionTime.region
145 timespan = predictedRegionTime.timespan
146 expCenter = SpherePoint(region.getBoundingCircle().getCenter())
147 expRadius = region.getBoundingCircle().getOpeningAngle().asDegrees()
148 expMidPointEPOCH = getMidpointFromTimespan(timespan, allowUnbounded=False).mjd
150 # MPSky service query
151 try:
152 mpSkyURL = os.environ['MP_SKY_URL']
153 except KeyError:
154 self.log.warning("MP_SKY_URL is not defined. Falling back to %s.", self.config.mpSkyFallbackURL)
155 mpSkyURL = self.config.mpSkyFallbackURL
156 mpSkySsObjects, elements = self._mpSkyConeSearch(expCenter, expMidPointEPOCH,
157 expRadius + self.config.queryBufferRadiusDegrees,
158 mpSkyURL)
159 if elements is not None:
160 for element in ELEMENTS_LIST:
161 mpSkySsObjects[self.config.mpcorb_prefix + NAMING_MAP[element]] = elements[element]
162 # ssObjectId and mpcDesignation are also in the MPCORB schema so let's add them here.
163 return Struct(
164 ssObjects=mpSkySsObjects,
165 )
167 def read_mp_sky_response(self, response):
168 """Extract ephemerides from an MPSky web request
170 Parameters
171 ----------
172 response : `requests.Response`
173 MPSky message
175 Returns
176 -------
177 objID : `np.ndarray`
178 Designations of nearby objects
179 ra : `np.ndarray`
180 Array of object right ascensions
181 dec : `np.ndarray`
182 Array of object declinations
183 object_polynomial : `np.ndarray`, (N,M)
184 Array of object cartesian position polynomials
185 observer_polynomial : `np.ndarray`, (N,M)
186 Array of observer cartesian position polynomials
187 t_min : `np.ndarray`
188 Lower time bound for polynomials
189 t_max : `np.ndarray`
190 Upper time bound for polynomials
191 elements : `pandas.DataFrame` or `NoneType`
192 Orbital elements of ssObjects (to pass to alerts)
193 """
194 with pa.input_stream(memoryview(response.content)) as fp:
195 fp.seek(0)
196 p = pa.ipc.read_tensor(fp)
197 op = pa.ipc.read_tensor(fp)
198 with pa.ipc.open_stream(fp) as reader:
199 schema = reader.schema
200 r = next(reader)
202 # construct an elements pandas dataframe
203 if "q" in schema.names:
204 cols = {col: r[col].to_numpy() for col in ELEMENTS_LIST}
205 elements = pd.DataFrame(cols)
206 else:
207 elements = None
209 return (
210 r["name"].to_numpy(zero_copy_only=False), r["ra"].to_numpy(), r["dec"].to_numpy(),
211 p.to_numpy(), op.to_numpy(), r["tmin"].to_numpy(), r["tmax"].to_numpy(), elements
212 )
214 def _mpSkyConeSearch(self, expCenter, epochMJD, queryRadius, mpSkyURL):
215 """Query MPSky ephemeris service for objects near the expected detector position
217 Parameters
218 ----------
219 expCenter : `lsst.geom.SpherePoint`
220 Center of search cone
221 epochMJD : `float`
222 Epoch of cone search, (MJD in UTC).
223 queryRadius : `float`
224 Radius of the cone search in degrees.
225 mpSkyURL : `str`
226 URL to query for MPSky.
228 Returns
229 -------
230 mpSkySsObjects : `pandas.DataFrame`
231 DataFrame with Solar System Object information and RA/DEC position
232 within the visit.
233 elements : `pandas.DataFrame` or `None`
234 Orbital elements of ssObjects (to pass to alerts)
235 """
236 fieldRA = expCenter.getRa().asDegrees()
237 fieldDec = expCenter.getDec().asDegrees()
239 params = {
240 "t": epochMJD,
241 "ra": fieldRA,
242 "dec": fieldDec,
243 "radius": queryRadius,
244 "return_elements": True
245 }
247 response = None
248 try:
249 response = requests.get(mpSkyURL, params=params, timeout=self.config.mpSkyRequestTimeoutSeconds)
250 response.raise_for_status()
251 response = self.read_mp_sky_response(response)
252 objID, ra, dec, object_polynomial, observer_polynomial, tmin, tmax, elements = response
254 mpSkySsObjects = pd.DataFrame()
255 mpSkySsObjects['ObjID'] = objID
256 mpSkySsObjects['ra'] = ra
257 mpSkySsObjects['obj_x_poly'] = [poly[0] for poly in object_polynomial.T]
258 mpSkySsObjects['obj_y_poly'] = [poly[1] for poly in object_polynomial.T]
259 mpSkySsObjects['obj_z_poly'] = [poly[2] for poly in object_polynomial.T]
260 mpSkySsObjects['obs_x_poly'] = [observer_polynomial.T[0] for
261 i in range(len(mpSkySsObjects))]
262 mpSkySsObjects['obs_y_poly'] = [observer_polynomial.T[1] for
263 i in range(len(mpSkySsObjects))]
264 mpSkySsObjects['obs_z_poly'] = [observer_polynomial.T[2] for
265 i in range(len(mpSkySsObjects))]
266 mpSkySsObjects['tmin'] = tmin
267 mpSkySsObjects['tmax'] = tmax
268 mpSkySsObjects['dec'] = dec
269 mpSkySsObjects['Err(arcsec)'] = 2
270 mpSkySsObjects['ssObjectId'] = [obj_id_to_ss_object_id(v) for v in mpSkySsObjects['ObjID'].values]
271 nFound = len(mpSkySsObjects)
273 if nFound == 0:
274 self.log.info("No Solar System objects found for visit.")
275 else:
276 self.log.info("%d Solar System Objects in visit", nFound)
277 except requests.RequestException as e:
278 body = getattr(response, "text", "<No response>")
279 message = f"Query to the remote ephemerides service failed. Got response {body}"
280 self.log.error(message)
281 raise NoWorkFound(f"{message}: {e}") from e
283 return mpSkySsObjects, elements