Coverage for python/lsst/pipe/tasks/generateEphemerides.py: 34%

1# This file is part of ap_association.

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.

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 os

48import pandas as pd

49from importlib import resources

50from sqlite3 import connect as sqlite_connect

51from subprocess import run, PIPE

52from tempfile import TemporaryDirectory

53from textwrap import dedent

56from lsst.pex.config import Field

57from lsst.pipe.base import connectionTypes, NoWorkFound, PipelineTask, \

58 PipelineTaskConfig, PipelineTaskConnections, Struct

59import lsst.pipe.tasks

60from lsst.utils.timer import timeMethod

63class GenerateEphemeridesConnections(PipelineTaskConnections,

64 dimensions=("instrument",)):

66 visitSummaries = connectionTypes.Input(

67 doc="Summary of visit information including ra, dec, and time",

68 name="preliminary_visit_summary",

69 storageClass="ExposureCatalog",

70 dimensions={"instrument", "visit"},

71 multiple=True

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 )

135

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 )

143

144

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

160

161class GenerateEphemeridesTask(PipelineTask):

162 """Task to query the Sorcha service and retrieve the solar system objects

163 that are observable within the input visit.

164 """

165 ConfigClass = GenerateEphemeridesConfig

166 _DefaultName = "GenerateEphemerides"

167

168 def runQuantum(

169 self,

170 butlerQC,

171 inputRefs,

172 outputRefs,

173 ):

174 """Do butler IO and transform to provide in memory

175 objects for tasks `~Task.run` method.

176

177 Parameters

178 ----------

179 butlerQC : `QuantumContext`

180 A butler which is specialized to operate in the context of a

181 `lsst.daf.butler.Quantum`.

182 inputRefs : `InputQuantizedConnection`

183 Datastructure whose attribute names are the names that identify

184 connections defined in corresponding `PipelineTaskConnections`

185 class. The values of these attributes are the

186 `lsst.daf.butler.DatasetRef` objects associated with the defined

187 input/prerequisite connections.

188 outputRefs : `OutputQuantizedConnection`

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 output connections.

194 """

195 inputs = butlerQC.get(inputRefs)

196 outputs = self.run(**inputs)

197 for ref in outputRefs.ssObjects:

198 dataId = ref.dataId

199 ephemeris_visit = outputs.ssObjects[dataId['visit']]

200 butlerQC.put(ephemeris_visit, ref)

201

202 @timeMethod

203 def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,

204 earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict, naif0012):

205 """Parse the information on the current visit and retrieve the

206 observable solar system objects from Sorcha.

207

208 Parameters

209 ----------

210 visitSummary : `lsst.afw.table.ExposureCatalog`

211 Has rows with .getVisitInfo, which give the center and time of exposure

212

213 mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,

214 earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict,

215 naif0012 : `lsst.pipe.tasks.sspAuxiliaryFile.SSPAuxiliaryFile`s

216 Minor Planet Center orbit table used for association

217

218 Returns

219 -------

220 result : `lsst.pipe.base.Struct`

221 Results struct with components:

222

223 - ``ssObjects`` : `pandas.DataFrame`

224 DataFrame containing Solar System Objects near the detector

225 footprint as retrieved by Sorcha. The columns are as follows:

226

227 ``Name``

228 object name (`str`)

229 ``ra``

230 RA in decimal degrees (`float`)

231 ``dec``

232 DEC in decimal degrees (`float`)

233

234 """

235 if len(visitSummaries) == 0:

236 raise NoWorkFound("No visits input!")

237 visitInfos = [vs[0].getVisitInfo() for vs in visitSummaries]

238 fieldRA = np.array([vi.getBoresightRaDec().getRa().asDegrees() for vi in visitInfos])

239 fieldDec = np.array([vi.getBoresightRaDec().getDec().asDegrees() for vi in visitInfos])

240 epochMJD = np.array([vi.date.toAstropy().tai.mjd for vi in visitInfos])

241 visits = [vi.id for vi in visitInfos]

242

243 # Confused about seconds/days units here

244 n = len(epochMJD)

245 visitTime = np.ones(n) * 30.0 # seconds

246 inputVisits = pd.DataFrame({

247 "observationMidpointMJD": epochMJD,

248 "fieldRA": fieldRA,

249 "fieldDec": fieldDec,

250 "observationId": visits,

251 "visitTime": np.ones(n) * visitTime,

252 "observationStartMJD": epochMJD - (visitTime / 2) / 86400.0,

253 "visitExposureTime": visitTime,

254

255 # The following columns are required by Socha but only used after ephemeris generation

256 "filter": ["r"] * n,

257 "seeingFwhmGeom": [-1] * n,

258 "seeingFwhmEff": [-1] * n,

259 "fiveSigmaDepth": [-1] * n,

260 "rotSkyPos": [-1] * n,

261 })

262

263 inputOrbits = mpcorb[

264 ['packed_primary_provisional_designation', 'q', 'e', 'i',

265 'node', 'argperi', 'peri_time', 'epoch_mjd']

266 ].rename(columns={'packed_primary_provisional_designation': 'ObjID', 'epoch_mjd': 'epochMJD_TDB',

267 'i': 'inc', 'argperi': 'argPeri', 'peri_time': 't_p_MJD_TDB'})

268 inputOrbits['FORMAT'] = 'COM'

269 # Colors exactly like jake's prep_input_colors

270 inputColors = inputOrbits[["ObjID"]].copy()

271 inputColors["H_r"] = mpcorb['h']

272 inputColors["GS"] = 0.15 # Traditional

273

274 eph_str = resources.files(lsst.pipe.tasks).parents[3].joinpath("data/eph.ini").read_text()

275 eph_str = eph_str.replace("{OBSCODE}", self.config.observatoryCode)

276 eph_str = eph_str.replace("{FOV}", str(self.config.observatoryFOVRadius))

277

278 with TemporaryDirectory() as tmpdirname:

279 self.log.info(f'temp dir: {tmpdirname}')

280

281 # Orbits

282 inputOrbits.to_csv(f'{tmpdirname}/orbits.csv', index=False)

283 # Observations SQLite

284 conn = sqlite_connect(f'{tmpdirname}/pointings.db')

285 inputVisits.to_sql('observations', conn, if_exists='replace', index=False)

286 conn.close()

287

288 with open(f'{tmpdirname}/eph.ini', 'w') as ephFile:

289 ephFile.write(eph_str)

290

291 inputColors.to_csv(f'{tmpdirname}/colors.csv', index=False)

292

293 cache = f'{tmpdirname}/sorcha_cache/'

294 self.log.info('making cache')

295 os.mkdir(cache)

296 # DONE

297 for filename, fileref in [

298 ('de440s.bsp', de440s),

299 ('sb441-n16.bsp', sb441_n16),

300 ('ObsCodes.json', obsCodes),

301 ('linux_p1550p2650.440', linux_p1550p2650),

302 ('pck00010.pck', pck00010),

303 ('earth_latest_high_prec.bpc', earth_latest_high_prec),

304 ('earth_620120_250826.bpc', earth_620120_250826),

305 ('earth_2025_250826_2125_predict.bpc', earth_2025_250826_2125_predict),

306 ('naif0012.tls', naif0012),

307 ]:

308 self.log.info(f'writing {filename}')

309 with open(cache + filename, 'wb') as file:

310 file.write(fileref.fileContents.read())

311

312 meta_kernel_text = dedent(f"""\

313 \\begindata

314

315 PATH_VALUES = ('{tmpdirname}/sorcha_cache/')

316

317 PATH_SYMBOLS = ('A')

318

319 KERNELS_TO_LOAD=(

320 '$A/naif0012.tls',

321 '$A/earth_620120_250826.bpc',

322 '$A/earth_2025_250826_2125_predict.bpc',

323 '$A/pck00010.pck',

324 '$A/de440s.bsp',

325 '$A/earth_latest_high_prec.bpc',

326 )

327

328 \\begintext

329 """)

330 with open(f'{tmpdirname}/sorcha_cache/meta_kernel.txt', 'w') as meta_kernel_file:

331 meta_kernel_file.write(meta_kernel_text)

332 self.log.info('Sorcha process begun')

333

334 result = run(

335 [

336 "sorcha",

337 "run",

338 "-c", f"{tmpdirname}/eph.ini",

339 "-o", f"{tmpdirname}/",

340 "--ob", f"{tmpdirname}/orbits.csv",

341 "-p", f"{tmpdirname}/colors.csv",

342 "--pd", f"{tmpdirname}/pointings.db",

343 "--ew", f"{tmpdirname}/ephemeris",

344 "--ar", f"{tmpdirname}/sorcha_cache/"

345 ],

346 stdout=PIPE,

347 stderr=PIPE,

348 text=True

349 )

350

351 self.log.info(f"Sorcha STDOUT:\n {result.stdout}")

352 self.log.info(f"Sorcha STDERR:\n {result.stderr}")

353

354 eph_path = f'{tmpdirname}/ephemeris.csv'

355 if not os.path.exists(eph_path):

356 raise FileNotFoundError(

357 " Sorcha did not create ephemeris. Check STDOUT/STDERR above. "

358 f"Directory contents:\n{os.listdir(tmpdirname)}"

359 )

360

361 # Return Sorcha output directly

362 ephemeris = pd.read_csv(eph_path)

363 perVisitSsObjects = {FieldID: group for FieldID, group in ephemeris.groupby("FieldID")}

364 for v in visits:

365 if v not in perVisitSsObjects:

366 perVisitSsObjects[v] = ephemeris.iloc[:0]

367 return Struct(ssObjects=perVisitSsObjects)