3232ephemerides from a csv. While awkward and un-pipetask-like, it works and takes
3333advantage of Sorcha's as-designed speed.
3434
35- Eventually, this should be replaced with adam_core's ephemeris generation
35+ Eventually, this should be replaced with adam_core's ephemeris generation
3636tools which propagate forward orbital uncertainty into on-sky ellipses.
3737Doing so will require re-implementing the healpix binning method described
3838in https://arxiv.org/abs/2506.02140. Doing so will not only improve this code
3939but also allow on-sky uncertainties to be used during association. When this is
40- done, mpsky should be modified to do the same.
40+ done, mpsky should be modified to do the same.
4141"""
4242
4343__all__ = ["GenerateEphemeridesConfig" , "GenerateEphemeridesTask" ]
4747import os
4848import pandas as pd
4949from importlib import resources
50- import sqlite3
51- import subprocess
52- import tempfile
50+ from sqlite3 import connect as sqlite_connect
51+ from subprocess import run , PIPE
52+ from tempfile import TemporaryDirectory
53+ from textwrap import dedent
5354
5455
55- import lsst .daf .butler as dafButler
56- from lsst .geom import SpherePoint , degrees
57- import lsst .pex .config as pexConfig
56+ from lsst .pex .config import Field
5857from lsst .pipe .base import connectionTypes , NoWorkFound , PipelineTask , \
5958 PipelineTaskConfig , PipelineTaskConnections , Struct
60- import lsst .pipe .tasks
61- from lsst .pipe .tasks .associationUtils import obj_id_to_ss_object_id
62- from lsst .sphgeom import ConvexPolygon
6359from lsst .utils .timer import timeMethod
6460
61+
6562class GenerateEphemeridesConnections (PipelineTaskConnections ,
66- dimensions = ("instrument" ,)):
63+ dimensions = ("instrument" ,)):
6764
6865 visitSummaries = connectionTypes .Input (
6966 doc = "Summary of visit information including ra, dec, and time" ,
7067 name = "preliminary_visit_summary" ,
7168 storageClass = "ExposureCatalog" ,
7269 dimensions = {"instrument" , "visit" },
73- multiple = True
70+ multiple = True
7471 )
7572 mpcorb = connectionTypes .Input (
7673 doc = "Minor Planet Center orbit table used for association" ,
@@ -142,19 +139,18 @@ class GenerateEphemeridesConnections(PipelineTaskConnections,
142139 dimensions = ("instrument" , "visit" ),
143140 multiple = True ,
144141 )
145-
146142
147143
148144class GenerateEphemeridesConfig (
149145 PipelineTaskConfig ,
150146 pipelineConnections = GenerateEphemeridesConnections ):
151- observatoryCode = pexConfig . Field (
147+ observatoryCode = Field (
152148 dtype = str ,
153149 doc = "IAU Minor Planet Center observer code for queries "
154150 "(default is X05 for Rubin Obs./LSST)" ,
155151 default = 'X05'
156152 )
157- observatoryFOVRadius = pexConfig . Field (
153+ observatoryFOVRadius = Field (
158154 dtype = float ,
159155 doc = "The field of view of the observatory (degrees)" ,
160156 default = 2.06 ,
@@ -202,7 +198,6 @@ def runQuantum(
202198 ephemeris_visit = outputs .ssObjects [dataId ['visit' ]]
203199 butlerQC .put (ephemeris_visit , ref )
204200
205-
206201 @timeMethod
207202 def run (self , visitSummaries , mpcorb , de440s , sb441_n16 , obsCodes , linux_p1550p2650 , pck00010 ,
208203 earth_latest_high_prec , earth_620120_250826 , earth_2025_250826_2125_predict , naif0012 ):
@@ -214,6 +209,11 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
214209 visitSummary : `lsst.afw.table.ExposureCatalog`
215210 Has rows with .getVisitInfo, which give the center and time of exposure
216211
212+ mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2650, pck00010,
213+ earth_latest_high_prec, earth_620120_250826, earth_2025_250826_2125_predict,
214+ naif0012 : `lsst.pipe.tasks.sspAuxiliaryFile.SSPAuxiliaryFile`s
215+ Minor Planet Center orbit table used for association
216+
217217 Returns
218218 -------
219219 result : `lsst.pipe.base.Struct`
@@ -229,7 +229,7 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
229229 RA in decimal degrees (`float`)
230230 ``dec``
231231 DEC in decimal degrees (`float`)
232-
232+
233233 """
234234 if len (visitSummaries ) == 0 :
235235 raise NoWorkFound ("No visits input!" )
@@ -247,7 +247,7 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
247247 "fieldRA" : fieldRA ,
248248 "fieldDec" : fieldDec ,
249249 "observationId" : visits ,
250- "visitTime" : np .ones (n ) * visitTime ,
250+ "visitTime" : np .ones (n ) * visitTime ,
251251 "observationStartMJD" : epochMJD - (visitTime / 2 ) / 86400.0 ,
252252 "visitExposureTime" : visitTime ,
253253
@@ -258,38 +258,39 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
258258 "fiveSigmaDepth" : [- 1 ] * n ,
259259 "rotSkyPos" : [- 1 ] * n ,
260260 })
261- mpcorb = mpcorb .iloc [:10_000 ] # TODO: remove
262261
263262 inputOrbits = mpcorb [
264- ['packed_primary_provisional_designation' , 'q' , 'e' , 'i' , 'node' ,'argperi' ,'peri_time' ,'epoch_mjd' ]
263+ ['packed_primary_provisional_designation' , 'q' , 'e' , 'i' ,
264+ 'node' , 'argperi' , 'peri_time' , 'epoch_mjd' ]
265265 ].rename (columns = {'packed_primary_provisional_designation' : 'ObjID' , 'epoch_mjd' : 'epochMJD_TDB' ,
266266 'i' : 'inc' , 'argperi' : 'argPeri' , 'peri_time' : 't_p_MJD_TDB' })
267267 inputOrbits ['FORMAT' ] = 'COM'
268268 # Colors exactly like jake's prep_input_colors
269269 inputColors = inputOrbits [["ObjID" ]].copy ()
270270 inputColors ["H_r" ] = mpcorb ['h' ]
271- inputColors ["GS" ] = 0.15 # Traditional
271+ inputColors ["GS" ] = 0.15 # Traditional
272272
273273 eph_str = resources .files (lsst .pipe .tasks ).parents [3 ].joinpath ("data/eph.ini" ).read_text ()
274274 eph_str = eph_str .replace ("{OBSCODE}" , self .config .observatoryCode )
275275 eph_str = eph_str .replace ("{FOV}" , str (self .config .observatoryFOVRadius ))
276276
277- with tempfile . TemporaryDirectory () as tmpdirname :
277+ with TemporaryDirectory () as tmpdirname :
278278 self .log .info (f'temp dir: { tmpdirname } ' )
279-
279+
280280 # Orbits
281281 inputOrbits .to_csv (f'{ tmpdirname } /orbits.csv' , index = False )
282282 # Observations SQLite
283- conn = sqlite3 . connect (f'{ tmpdirname } /pointings.db' )
283+ conn = sqlite_connect (f'{ tmpdirname } /pointings.db' )
284284 inputVisits .to_sql ('observations' , conn , if_exists = 'replace' , index = False )
285285 conn .close ()
286286
287- open (f'{ tmpdirname } /eph.ini' , 'w' ).write (eph_str )
287+ with open (f'{ tmpdirname } /eph.ini' , 'w' ) as ephFile :
288+ ephFile .write (eph_str )
289+
288290 inputColors .to_csv (f'{ tmpdirname } /colors.csv' , index = False )
289291
290-
291292 cache = f'{ tmpdirname } /sorcha_cache/'
292- self .log .info (f 'making cache' )
293+ self .log .info ('making cache' )
293294 os .mkdir (cache )
294295 # DONE
295296 for filename , fileref in [
@@ -304,28 +305,32 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
304305 ('naif0012.tls' , naif0012 ),
305306 ]:
306307 self .log .info (f'writing { filename } ' )
307- open (cache + filename , 'wb' ). write ( fileref . fileContents . read ())
308- meta_kernel_text = f""" \\ begindata
308+ with open (cache + filename , 'wb' ) as file :
309+ file . write ( fileref . fileContents . read ())
309310
310- PATH_VALUES = ('{ tmpdirname } /sorcha_cache/')
311+ meta_kernel_text = dedent (f"""\
312+ \\ begindata
311313
312- PATH_SYMBOLS = ('A ')
314+ PATH_VALUES = ('{ tmpdirname } /sorcha_cache/ ')
313315
314- KERNELS_TO_LOAD=(
315- '$A/naif0012.tls',
316- '$A/earth_620120_250826.bpc',
317- '$A/earth_2025_250826_2125_predict.bpc',
318- '$A/pck00010.pck',
319- '$A/de440s.bsp',
320- '$A/earth_latest_high_prec.bpc',
321- )
316+ PATH_SYMBOLS = ('A')
322317
323- \\ begintext
324- """
325- open (f'{ tmpdirname } /sorcha_cache/meta_kernel.txt' , 'w' ).write (meta_kernel_text )
318+ KERNELS_TO_LOAD=(
319+ '$A/naif0012.tls',
320+ '$A/earth_620120_250826.bpc',
321+ '$A/earth_2025_250826_2125_predict.bpc',
322+ '$A/pck00010.pck',
323+ '$A/de440s.bsp',
324+ '$A/earth_latest_high_prec.bpc',
325+ )
326+
327+ \\ begintext
328+ """ )
329+ with open (f'{ tmpdirname } /sorcha_cache/meta_kernel.txt' , 'w' ) as meta_kernel_file :
330+ meta_kernel_file .write (meta_kernel_text )
326331 self .log .info ('Sorcha process begun' )
327332
328- result = subprocess . run (
333+ result = run (
329334 [
330335 "sorcha" ,
331336 "run" ,
@@ -337,18 +342,19 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
337342 "--ew" , f"{ tmpdirname } /ephemeris" ,
338343 "--ar" , f"{ tmpdirname } /sorcha_cache/"
339344 ],
340- stdout = subprocess . PIPE ,
341- stderr = subprocess . PIPE ,
345+ stdout = PIPE ,
346+ stderr = PIPE ,
342347 text = True
343348 )
344349
345350 self .log .info (f"Sorcha STDOUT:\n { result .stdout } " )
346351 self .log .info (f"Sorcha STDERR:\n { result .stderr } " )
347-
352+
348353 eph_path = f'{ tmpdirname } /ephemeris.csv'
349354 if not os .path .exists (eph_path ):
350355 raise FileNotFoundError (
351- f" Sorcha did not create ephemeris. Check STDOUT/STDERR above. Directory contents:\n { os .listdir (tmpdirname )} "
356+ " Sorcha did not create ephemeris. Check STDOUT/STDERR above. "
357+ f"Directory contents:\n { os .listdir (tmpdirname )} "
352358 )
353359
354360 # Return Sorcha output directly
@@ -358,7 +364,3 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
358364 if v not in perVisitSsObjects :
359365 perVisitSsObjects [v ] = ephemeris .iloc [:0 ]
360366 return Struct (ssObjects = perVisitSsObjects )
361-
362-
363-
364-
0 commit comments