Coverage for python / lsst / pipe / tasks / generateEphemerides.py: 28%
117 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-09 01:51 -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"""Generate per-visit solar system ephemerides using Sorcha.
24Sorcha is an open-source solar system survey simulation tool, which
25efficiently generates ephemerides for solar system objects from input orbits
26and observation pointings. It is currently the fastest robust and maintained
27ephemeris generator, and can fit our use case well.
29Sorcha is a command line tool, and is not designed for direct use in python.
30This task creates a temporary directory in which to run Sorcha as designed,
31providing it with all required input files and reading the Sorcha-generated
32ephemerides from a csv. While awkward and un-pipetask-like, it works and takes
33advantage of Sorcha's as-designed speed.
35Eventually, this should be replaced with adam_core's ephemeris generation
36tools which propagate forward orbital uncertainty into on-sky ellipses.
37Doing so will require re-implementing the healpix binning method described
38in https://arxiv.org/abs/2506.02140. Doing so will not only improve this code
39but also allow on-sky uncertainties to be used during association. When this is
40done, mpsky should be modified to do the same.
41"""
43__all__ = ["GenerateEphemeridesConfig", "GenerateEphemeridesTask"]
46import numpy as np
47import glob
48from importlib import resources
49import os
50import pandas as pd
51from sqlite3 import connect as sqlite_connect
52from subprocess import Popen, PIPE
53from tempfile import TemporaryDirectory
54from textwrap import dedent
55from time import sleep
57from lsst.pex.config import Field
58from lsst.pipe.base import connectionTypes, NoWorkFound, PipelineTask, \
59 PipelineTaskConfig, PipelineTaskConnections, Struct
60import lsst.pipe.tasks
61from lsst.utils.timer import timeMethod
64class GenerateEphemeridesConnections(PipelineTaskConnections,
65 dimensions=("instrument",)):
67 visitTable = connectionTypes.Input(
68 doc="Table summarizing each visit, including ra, dec, and time",
69 name="preliminary_visit_table",
70 storageClass="ArrowAstropy",
71 dimensions={"instrument"},
72 )
73 mpcorb = connectionTypes.Input(
74 doc="Minor Planet Center orbit table used for association",
75 name="mpcorb",
76 storageClass="DataFrame",
77 dimensions={},
78 )
80 # The following 9 prequisite inputs are Sorcha's required auxiliary files.
81 de440s = connectionTypes.PrerequisiteInput(
82 doc="NAIF DE440 ephemeris file (de440s.bsp)",
83 name="de440s",
84 storageClass="SSPAuxiliaryFile",
85 dimensions={},
86 )
87 sb441_n16 = connectionTypes.PrerequisiteInput(
88 doc="NAIF DE440 ephemeris file (sb441_n16.bsp)",
89 name="sb441_n16",
90 storageClass="SSPAuxiliaryFile",
91 dimensions={},
92 )
93 obsCodes = connectionTypes.PrerequisiteInput(
94 doc="MPC observatory code file (ObsCodes.json)",
95 name="obscodes",
96 storageClass="SSPAuxiliaryFile",
97 dimensions={},
98 )
99 linux_p1550p2650 = connectionTypes.PrerequisiteInput(
100 doc="TODO (linux_p1550p2650.440)",
101 name="linux_p1550p2650",
102 storageClass="SSPAuxiliaryFile",
103 dimensions={},
104 )
105 pck00010 = connectionTypes.PrerequisiteInput(
106 doc="orientation of planets, moons, the Sun, and selected asteroids. (pck00010.pck)",
107 name="pck00010",
108 storageClass="SSPAuxiliaryFile",
109 dimensions={},
110 )
111 earth_latest_high_prec = connectionTypes.PrerequisiteInput(
112 doc="High-precision Earth orientation parameters (EOP) kernel",
113 name="earth_latest_high_prec",
114 storageClass="SSPAuxiliaryFile",
115 dimensions={},
116 )
117 earth_620120_250826 = connectionTypes.PrerequisiteInput(
118 doc="Historical EOP",
119 name="earth_620120_250826",
120 storageClass="SSPAuxiliaryFile",
121 dimensions={},
122 )
123 earth_2025_250826_2125_predict = connectionTypes.PrerequisiteInput(
124 doc="Longterm EOP predictions",
125 name="earth_2025_250826_2125_predict",
126 storageClass="SSPAuxiliaryFile",
127 dimensions={},
128 )
129 naif0012 = connectionTypes.PrerequisiteInput(
130 doc="Leapsecond tls file",
131 name="naif0012",
132 storageClass="SSPAuxiliaryFile",
133 dimensions={},
134 )
136 ssObjects = connectionTypes.Output(
137 doc="Sorcha-provided Solar System objects observable in this detector-visit",
138 name="preloaded_ss_object_visit",
139 storageClass="DataFrame",
140 dimensions=("instrument", "visit"),
141 multiple=True,
142 )
145class GenerateEphemeridesConfig(
146 PipelineTaskConfig,
147 pipelineConnections=GenerateEphemeridesConnections):
148 observatoryCode = Field(
149 dtype=str,
150 doc="IAU Minor Planet Center observer code for queries "
151 "(default is X05 for Rubin Obs./LSST)",
152 default='X05'
153 )
154 observatoryFOVRadius = Field(
155 dtype=float,
156 doc="The field of view of the observatory (degrees)",
157 default=2.06,
158 )
159 loggingSleepTime = Field(
160 dtype=int,
161 doc="Time (seconds) to sleep between checking Sorcha's progress",
162 default=60,
163 )
166class GenerateEphemeridesTask(PipelineTask):
167 """Task to query the Sorcha service and retrieve the solar system objects
168 that are observable within the input visit.
169 """
170 ConfigClass = GenerateEphemeridesConfig
171 _DefaultName = "GenerateEphemerides"
173 def runQuantum(
174 self,
175 butlerQC,
176 inputRefs,
177 outputRefs,
178 ):
179 """Do butler IO and transform to provide in memory
180 objects for tasks `~Task.run` method.
182 Parameters
183 ----------
184 butlerQC : `QuantumContext`
185 A butler which is specialized to operate in the context of a
186 `lsst.daf.butler.Quantum`.
187 inputRefs : `InputQuantizedConnection`
188 Datastructure whose attribute names are the names that identify
189 connections defined in corresponding `PipelineTaskConnections`
190 class. The values of these attributes are the
191 `lsst.daf.butler.DatasetRef` objects associated with the defined
192 input/prerequisite connections.
193 outputRefs : `OutputQuantizedConnection`
194 Datastructure whose attribute names are the names that identify
195 connections defined in corresponding `PipelineTaskConnections`
196 class. The values of these attributes are the
197 `lsst.daf.butler.DatasetRef` objects associated with the defined
198 output connections.
199 """
200 inputs = butlerQC.get(inputRefs)
201 outputs = self.run(**inputs)
202 for ref in outputRefs.ssObjects:
203 dataId = ref.dataId
204 ephemeris_visit = outputs.ssObjects[dataId['visit']]
205 butlerQC.put(ephemeris_visit, ref)
207 @timeMethod
208 def run(self, visitTable, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,
209 earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict, naif0012):
210 """Parse the information on the current visit and retrieve the
211 observable solar system objects from Sorcha.
213 Parameters
214 ----------
215 visitSummary : `lsst.afw.table.ExposureCatalog`
216 Has rows with .getVisitInfo, which give the center and time of exposure
218 mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,
219 earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict,
220 naif0012 : `lsst.pipe.tasks.sspAuxiliaryFile.SSPAuxiliaryFile`s
221 Minor Planet Center orbit table used for association
223 Returns
224 -------
225 result : `lsst.pipe.base.Struct`
226 Results struct with components:
228 - ``ssObjects`` : `pandas.DataFrame`
229 DataFrame containing Solar System Objects near the detector
230 footprint as retrieved by Sorcha. The columns are as follows:
232 ``Name``
233 object name (`str`)
234 ``ra``
235 RA in decimal degrees (`float`)
236 ``dec``
237 DEC in decimal degrees (`float`)
239 """
240 if len(visitTable) == 0:
241 raise NoWorkFound("No visits input!")
242 visitTable = visitTable.to_pandas()
243 (
244 visits, epochMJD, fieldRA, fieldDec, visitTime
245 ) = visitTable[
246 ['visitId', 'expMidptMJD', 'ra', 'dec', 'expTime']
247 ].values.T
249 # Confused about seconds/days units here
250 n = len(epochMJD)
251 inputVisits = pd.DataFrame({
252 "observationMidpointMJD": epochMJD,
253 "fieldRA": fieldRA,
254 "fieldDec": fieldDec,
255 "observationId": visits,
256 "visitTime": np.ones(n) * visitTime,
257 "observationStartMJD": epochMJD - (visitTime / 2) / 86400.0,
258 "visitExposureTime": visitTime,
260 # The following columns are required by Socha but only used after ephemeris generation
261 "filter": ["r"] * n, # Let's call V-band "r-band" to play nicely with Sorcha...
262 "seeingFwhmGeom": [0.01] * n,
263 "seeingFwhmEff": [0.01] * n,
264 "fiveSigmaDepth": [9999] * n,
265 "rotSkyPos": [0] * n,
266 })
268 desig_dict = {un: pack for pack, un in mpcorb[['packed_primary_provisional_designation',
269 'unpacked_primary_provisional_designation']].values}
270 inputOrbits = mpcorb[
271 ['unpacked_primary_provisional_designation', 'q', 'e', 'i',
272 'node', 'argperi', 'peri_time', 'epoch_mjd']
273 ].rename(columns={'unpacked_primary_provisional_designation': 'ObjID', 'epoch_mjd': 'epochMJD_TDB',
274 'i': 'inc', 'argperi': 'argPeri', 'peri_time': 't_p_MJD_TDB'})
275 inputOrbits['FORMAT'] = 'COM'
276 # Colors exactly like jake's prep_input_colors
277 inputColors = inputOrbits[["ObjID"]].copy()
278 inputColors["H_r"] = mpcorb['h']
279 inputColors["GS"] = 0.15 # Traditional
281 eph_str = resources.files(lsst.pipe.tasks).parents[3].joinpath("data/eph.ini").read_text()
282 eph_str = eph_str.replace("{OBSCODE}", self.config.observatoryCode)
283 eph_str = eph_str.replace("{FOV}", str(self.config.observatoryFOVRadius))
285 with TemporaryDirectory() as tmpdirname:
286 self.log.info(f'temp dir: {tmpdirname}')
288 # Observations SQLite
289 conn = sqlite_connect(f'{tmpdirname}/pointings.db')
290 inputVisits.to_sql('observations', conn, if_exists='replace', index=False)
291 conn.close()
293 with open(f'{tmpdirname}/eph.ini', 'w') as ephFile:
294 ephFile.write(eph_str)
296 cache = f'{tmpdirname}/sorcha_cache/'
297 self.log.info('making cache')
298 os.mkdir(cache)
299 for filename, fileref in [
300 ('de440s.bsp', de440s),
301 ('sb441-n16.bsp', sb441_n16),
302 ('ObsCodes.json', obsCodes),
303 ('linux_p1550p2650.440', linux_p1550p2650),
304 ('pck00010.pck', pck00010),
305 ('earth_latest_high_prec.bpc', earth_latest_high_prec),
306 ('earth_620120_250826.bpc', earth_620120_250826),
307 ('earth_2025_250826_2125_predict.bpc', earth_2025_250826_2125_predict),
308 ('naif0012.tls', naif0012),
309 ]:
310 self.log.info(f'writing {filename}')
311 with open(cache + filename, 'wb') as file:
312 file.write(fileref.fileContents.read())
314 abspath = f'{tmpdirname}/sorcha_cache/'
315 split = 79
316 n_iter = int(len(str(abspath)) / split) # Number of splits required
317 for n in range(n_iter, 0, -1):
318 abspath = abspath[:split * n] + "+' '" + abspath[split * n:]
320 meta_kernel_text = dedent(f"""\
321 \\begindata
323 PATH_VALUES = ('{abspath}')
325 PATH_SYMBOLS = ('A')
327 KERNELS_TO_LOAD=(
328 '$A/naif0012.tls',
329 '$A/earth_620120_250826.bpc',
330 '$A/earth_2025_250826_2125_predict.bpc',
331 '$A/pck00010.pck',
332 '$A/de440s.bsp',
333 '$A/earth_latest_high_prec.bpc',
334 )
336 \\begintext
337 """)
338 with open(f'{tmpdirname}/sorcha_cache/meta_kernel.txt', 'w') as meta_kernel_file:
339 meta_kernel_file.write(meta_kernel_text)
340 self.log.info('Sorcha process begun')
342 # Orbits and colors
343 inputOrbits.to_csv(f'{tmpdirname}/orbits.csv', index=False)
344 inputColors.to_csv(f'{tmpdirname}/colors.csv', index=False)
346 # FIXME: this is a workaround for breakage in sbpy 0.5.0 caused
347 # by astropy 7.2.0. DM-53467 once sbpy is updated.
348 try:
349 import sbpy.photometry # noqa: F401
350 sorcha_run = ["sorcha", "run"]
351 except ImportError:
352 import textwrap
353 self.log.info("WARNING: Running sorcha-run patched for sbpy issues.")
354 sorcha_run = ["python", "-c", textwrap.dedent("""
355 import re, sys
356 print("WARNING: Running sorcha-run patched for sbpy issues.", file=sys.stderr)
358 import astropy.units.utils
359 import astropy.units.docgen
360 astropy.units.utils.generate_unit_summary = \
361 astropy.units.docgen.generate_unit_summary
363 from sorcha_cmdline.run import main
364 sys.exit(main())
365 """)]
367 proc = Popen(
368 sorcha_run
369 + [
370 "-c", f"{tmpdirname}/eph.ini",
371 "-o", f"{tmpdirname}/",
372 "--ob", f"{tmpdirname}/orbits.csv",
373 "-p", f"{tmpdirname}/colors.csv",
374 "--pd", f"{tmpdirname}/pointings.db",
375 "--ew", f"{tmpdirname}/ephemeris",
376 "--ar", f"{tmpdirname}/sorcha_cache/",
377 "-t", "sorcha_output",
378 ],
379 stdout=PIPE,
380 stderr=PIPE,
381 text=True
382 )
383 # As long as Sorcha is running, print its output or 'still running' messages
384 outfile = None
385 n_iter = 0
386 while proc.poll() is None:
387 sleep(self.config.loggingSleepTime)
388 n_iter += 1
389 if outfile is None:
390 outfile = open(glob.glob(f'{tmpdirname}/sorcha_output*.log')[0], 'r')
391 new_text = outfile.read()
392 if len(new_text) > 0:
393 self.log.info(new_text)
394 else:
395 seconds_elapsed = n_iter * self.config.loggingSleepTime
396 self.log.info(f"Sorcha is still running ({seconds_elapsed} seconds)")
398 stdout, stderr = proc.communicate()
399 self.log.info(f"Sorcha STDOUT:\n {stdout}")
400 self.log.info(f"Sorcha STDERR:\n {stderr}")
402 eph_path = f'{tmpdirname}/sorcha_output.csv'
403 if not os.path.exists(eph_path):
404 raise FileNotFoundError(
405 " Sorcha did not create sorcha_output. Check STDOUT/STDERR above. "
406 f"Directory contents:\n{os.listdir(tmpdirname)}"
407 )
409 # Return Sorcha output directly
410 ephemeris = pd.read_csv(eph_path)
411 ephemeris['packed_desig'] = ephemeris['ObjID'].map(desig_dict)
412 perVisitSsObjects = {FieldID: group for FieldID, group in ephemeris.groupby("FieldID")}
413 for v in visits:
414 if v not in perVisitSsObjects:
415 perVisitSsObjects[v] = ephemeris.iloc[:0]
416 return Struct(ssObjects=perVisitSsObjects)