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