Coverage for python / lsst / pipe / tasks / generateEphemerides.py: 27%
120 statements
« prev ^ index » next coverage.py v7.13.4, created at 2026-03-01 02:17 -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 )
79 # QG generation checks input visits to decide what output visits to expect.
80 # This dummy per-visit input helps QG generation decide to only expect output
81 # visits for visits actually present in the run, which this task gets in the
82 # not-per-visit visitTable input. See discussion on DM-54062 for more details.
83 dummyPerVisitInput = connectionTypes.Input(
84 doc="Dummy per-visit input to help trim QG generation",
85 name="visit_summary",
86 storageClass="ExposureCatalog",
87 dimensions={"instrument", "visit"},
88 multiple=True,
89 deferLoad=True,
90 )
92 # The following 9 prequisite inputs are Sorcha's required auxiliary files.
93 de440s = connectionTypes.PrerequisiteInput(
94 doc="NAIF DE440 ephemeris file (de440s.bsp)",
95 name="de440s",
96 storageClass="SSPAuxiliaryFile",
97 dimensions={},
98 )
99 sb441_n16 = connectionTypes.PrerequisiteInput(
100 doc="NAIF DE440 ephemeris file (sb441_n16.bsp)",
101 name="sb441_n16",
102 storageClass="SSPAuxiliaryFile",
103 dimensions={},
104 )
105 obsCodes = connectionTypes.PrerequisiteInput(
106 doc="MPC observatory code file (ObsCodes.json)",
107 name="obscodes",
108 storageClass="SSPAuxiliaryFile",
109 dimensions={},
110 )
111 linux_p1550p2650 = connectionTypes.PrerequisiteInput(
112 doc="TODO (linux_p1550p2650.440)",
113 name="linux_p1550p2650",
114 storageClass="SSPAuxiliaryFile",
115 dimensions={},
116 )
117 pck00010 = connectionTypes.PrerequisiteInput(
118 doc="orientation of planets, moons, the Sun, and selected asteroids. (pck00010.pck)",
119 name="pck00010",
120 storageClass="SSPAuxiliaryFile",
121 dimensions={},
122 )
123 earth_latest_high_prec = connectionTypes.PrerequisiteInput(
124 doc="High-precision Earth orientation parameters (EOP) kernel",
125 name="earth_latest_high_prec",
126 storageClass="SSPAuxiliaryFile",
127 dimensions={},
128 )
129 earth_620120_250826 = connectionTypes.PrerequisiteInput(
130 doc="Historical EOP",
131 name="earth_620120_250826",
132 storageClass="SSPAuxiliaryFile",
133 dimensions={},
134 )
135 earth_2025_250826_2125_predict = connectionTypes.PrerequisiteInput(
136 doc="Longterm EOP predictions",
137 name="earth_2025_250826_2125_predict",
138 storageClass="SSPAuxiliaryFile",
139 dimensions={},
140 )
141 naif0012 = connectionTypes.PrerequisiteInput(
142 doc="Leapsecond tls file",
143 name="naif0012",
144 storageClass="SSPAuxiliaryFile",
145 dimensions={},
146 )
148 ssObjects = connectionTypes.Output(
149 doc="Sorcha-provided Solar System objects observable in this detector-visit",
150 name="preloaded_ss_object_visit",
151 storageClass="DataFrame",
152 dimensions=("instrument", "visit"),
153 multiple=True,
154 )
157class GenerateEphemeridesConfig(
158 PipelineTaskConfig,
159 pipelineConnections=GenerateEphemeridesConnections):
160 observatoryCode = Field(
161 dtype=str,
162 doc="IAU Minor Planet Center observer code for queries "
163 "(default is X05 for Rubin Obs./LSST)",
164 default='X05'
165 )
166 observatoryFOVRadius = Field(
167 dtype=float,
168 doc="The field of view of the observatory (degrees)",
169 default=2.06,
170 )
171 loggingSleepTime = Field(
172 dtype=int,
173 doc="Time (seconds) to sleep between checking Sorcha's progress",
174 default=60,
175 )
178class GenerateEphemeridesTask(PipelineTask):
179 """Task to query the Sorcha service and retrieve the solar system objects
180 that are observable within the input visit.
181 """
182 ConfigClass = GenerateEphemeridesConfig
183 _DefaultName = "GenerateEphemerides"
185 def runQuantum(
186 self,
187 butlerQC,
188 inputRefs,
189 outputRefs,
190 ):
191 """Do butler IO and transform to provide in memory
192 objects for tasks `~Task.run` method.
194 Parameters
195 ----------
196 butlerQC : `QuantumContext`
197 A butler which is specialized to operate in the context of a
198 `lsst.daf.butler.Quantum`.
199 inputRefs : `InputQuantizedConnection`
200 Datastructure whose attribute names are the names that identify
201 connections defined in corresponding `PipelineTaskConnections`
202 class. The values of these attributes are the
203 `lsst.daf.butler.DatasetRef` objects associated with the defined
204 input/prerequisite connections.
205 outputRefs : `OutputQuantizedConnection`
206 Datastructure whose attribute names are the names that identify
207 connections defined in corresponding `PipelineTaskConnections`
208 class. The values of these attributes are the
209 `lsst.daf.butler.DatasetRef` objects associated with the defined
210 output connections.
211 """
212 inputs = butlerQC.get(inputRefs)
213 del inputs["dummyPerVisitInput"]
214 outputs = self.run(**inputs)
215 for ref in outputRefs.ssObjects:
216 dataId = ref.dataId
217 if dataId['visit'] in outputs.ssObjects:
218 ephemeris_visit = outputs.ssObjects[dataId['visit']]
219 butlerQC.put(ephemeris_visit, ref)
221 @timeMethod
222 def run(self, visitTable, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,
223 earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict, naif0012):
224 """Parse the information on the current visit and retrieve the
225 observable solar system objects from Sorcha.
227 Parameters
228 ----------
229 visitSummary : `lsst.afw.table.ExposureCatalog`
230 Has rows with .getVisitInfo, which give the center and time of exposure
232 mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,
233 earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict,
234 naif0012 : `lsst.pipe.tasks.sspAuxiliaryFile.SSPAuxiliaryFile`s
235 Minor Planet Center orbit table used for association
237 Returns
238 -------
239 result : `lsst.pipe.base.Struct`
240 Results struct with components:
242 - ``ssObjects`` : `pandas.DataFrame`
243 DataFrame containing Solar System Objects near the detector
244 footprint as retrieved by Sorcha. The columns are as follows:
246 ``Name``
247 object name (`str`)
248 ``ra``
249 RA in decimal degrees (`float`)
250 ``dec``
251 DEC in decimal degrees (`float`)
253 """
254 if len(visitTable) == 0:
255 raise NoWorkFound("No visits input!")
256 visitTable = visitTable.to_pandas()
257 (
258 visits, epochMJD, fieldRA, fieldDec, visitTime
259 ) = visitTable[
260 ['visitId', 'expMidptMJD', 'ra', 'dec', 'expTime']
261 ].values.T
263 # Confused about seconds/days units here
264 n = len(epochMJD)
265 inputVisits = pd.DataFrame({
266 "observationMidpointMJD": epochMJD,
267 "fieldRA": fieldRA,
268 "fieldDec": fieldDec,
269 "observationId": visits,
270 "visitTime": np.ones(n) * visitTime,
271 "observationStartMJD": epochMJD - (visitTime / 2) / 86400.0,
272 "visitExposureTime": visitTime,
274 # The following columns are required by Socha but only used after ephemeris generation
275 "filter": ["r"] * n, # Let's call V-band "r-band" to play nicely with Sorcha...
276 "seeingFwhmGeom": [0.01] * n,
277 "seeingFwhmEff": [0.01] * n,
278 "fiveSigmaDepth": [9999] * n,
279 "rotSkyPos": [0] * n,
280 })
282 desig_dict = {un: pack for pack, un in mpcorb[['packed_primary_provisional_designation',
283 'unpacked_primary_provisional_designation']].values}
284 inputOrbits = mpcorb[
285 ['unpacked_primary_provisional_designation', 'q', 'e', 'i',
286 'node', 'argperi', 'peri_time', 'epoch_mjd']
287 ].rename(columns={'unpacked_primary_provisional_designation': 'ObjID', 'epoch_mjd': 'epochMJD_TDB',
288 'i': 'inc', 'argperi': 'argPeri', 'peri_time': 't_p_MJD_TDB'})
289 inputOrbits['FORMAT'] = 'COM'
290 # Colors exactly like jake's prep_input_colors
291 inputColors = inputOrbits[["ObjID"]].copy()
292 inputColors["H_r"] = mpcorb['h']
293 inputColors["GS"] = 0.15 # Traditional
295 eph_str = resources.files(lsst.pipe.tasks).parents[3].joinpath("data/eph.ini").read_text()
296 eph_str = eph_str.replace("{OBSCODE}", self.config.observatoryCode)
297 eph_str = eph_str.replace("{FOV}", str(self.config.observatoryFOVRadius))
299 with TemporaryDirectory() as tmpdirname:
300 self.log.info(f'temp dir: {tmpdirname}')
302 # Observations SQLite
303 conn = sqlite_connect(f'{tmpdirname}/pointings.db')
304 inputVisits.to_sql('observations', conn, if_exists='replace', index=False)
305 conn.close()
307 with open(f'{tmpdirname}/eph.ini', 'w') as ephFile:
308 ephFile.write(eph_str)
310 cache = f'{tmpdirname}/sorcha_cache/'
311 self.log.info('making cache')
312 os.mkdir(cache)
313 for filename, fileref in [
314 ('de440s.bsp', de440s),
315 ('sb441-n16.bsp', sb441_n16),
316 ('ObsCodes.json', obsCodes),
317 ('linux_p1550p2650.440', linux_p1550p2650),
318 ('pck00010.pck', pck00010),
319 ('earth_latest_high_prec.bpc', earth_latest_high_prec),
320 ('earth_620120_250826.bpc', earth_620120_250826),
321 ('earth_2025_250826_2125_predict.bpc', earth_2025_250826_2125_predict),
322 ('naif0012.tls', naif0012),
323 ]:
324 self.log.info(f'writing {filename}')
325 with open(cache + filename, 'wb') as file:
326 file.write(fileref.fileContents.read())
328 abspath = f'{tmpdirname}/sorcha_cache/'
329 split = 79
330 n_iter = int(len(str(abspath)) / split) # Number of splits required
331 for n in range(n_iter, 0, -1):
332 abspath = abspath[:split * n] + "+' '" + abspath[split * n:]
334 meta_kernel_text = dedent(f"""\
335 \\begindata
337 PATH_VALUES = ('{abspath}')
339 PATH_SYMBOLS = ('A')
341 KERNELS_TO_LOAD=(
342 '$A/naif0012.tls',
343 '$A/earth_620120_250826.bpc',
344 '$A/earth_2025_250826_2125_predict.bpc',
345 '$A/pck00010.pck',
346 '$A/de440s.bsp',
347 '$A/earth_latest_high_prec.bpc',
348 )
350 \\begintext
351 """)
352 with open(f'{tmpdirname}/sorcha_cache/meta_kernel.txt', 'w') as meta_kernel_file:
353 meta_kernel_file.write(meta_kernel_text)
354 self.log.info('Sorcha process begun')
356 # Orbits and colors
357 inputOrbits.to_csv(f'{tmpdirname}/orbits.csv', index=False)
358 inputColors.to_csv(f'{tmpdirname}/colors.csv', index=False)
360 # FIXME: this is a workaround for breakage in sbpy 0.5.0 caused
361 # by astropy 7.2.0. DM-53467 once sbpy is updated.
362 try:
363 import sbpy.photometry # noqa: F401
364 sorcha_run = ["sorcha", "run"]
365 except ImportError:
366 import textwrap
367 self.log.info("WARNING: Running sorcha-run patched for sbpy issues.")
368 sorcha_run = ["python", "-c", textwrap.dedent("""
369 import re, sys
370 print("WARNING: Running sorcha-run patched for sbpy issues.", file=sys.stderr)
372 import astropy.units.utils
373 import astropy.units.docgen
374 astropy.units.utils.generate_unit_summary = \
375 astropy.units.docgen.generate_unit_summary
377 from sorcha_cmdline.run import main
378 sys.exit(main())
379 """)]
380 proc = Popen(
381 sorcha_run
382 + [
383 "-c", f"{tmpdirname}/eph.ini",
384 "-o", f"{tmpdirname}/",
385 "--ob", f"{tmpdirname}/orbits.csv",
386 "-p", f"{tmpdirname}/colors.csv",
387 "--pd", f"{tmpdirname}/pointings.db",
388 "--ew", f"{tmpdirname}/ephemeris",
389 "--ar", f"{tmpdirname}/sorcha_cache/",
390 "-t", "sorcha_output",
391 ],
392 stdout=PIPE,
393 stderr=PIPE,
394 text=True
395 )
396 # As long as Sorcha is running, print its output or 'still running' messages
397 outfile = None
398 n_iter = 0
399 while proc.poll() is None:
400 sleep(self.config.loggingSleepTime)
401 n_iter += 1
402 if outfile is None:
403 outfile = open(glob.glob(f'{tmpdirname}/sorcha_output*.log')[0], 'r')
404 new_text = outfile.read()
405 if len(new_text) > 0:
406 self.log.info(new_text)
407 else:
408 seconds_elapsed = n_iter * self.config.loggingSleepTime
409 self.log.info(f"Sorcha is still running ({seconds_elapsed} seconds)")
411 stdout, stderr = proc.communicate()
412 self.log.info(f"Sorcha STDOUT:\n {stdout}")
413 self.log.info(f"Sorcha STDERR:\n {stderr}")
415 eph_path = f'{tmpdirname}/sorcha_output.csv'
416 if not os.path.exists(eph_path):
417 raise FileNotFoundError(
418 " Sorcha did not create sorcha_output. Check STDOUT/STDERR above. "
419 f"Directory contents:\n{os.listdir(tmpdirname)}"
420 )
422 # Return Sorcha output directly
423 ephemeris = pd.read_csv(eph_path)
424 ephemeris['packed_desig'] = ephemeris['ObjID'].map(desig_dict)
425 perVisitSsObjects = {FieldID: group for FieldID, group in ephemeris.groupby("FieldID")}
426 for v in visits:
427 if v not in perVisitSsObjects:
428 perVisitSsObjects[v] = ephemeris.iloc[:0]
429 return Struct(ssObjects=perVisitSsObjects)