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
6059import lsst .pipe .tasks
61- from lsst .pipe .tasks .associationUtils import obj_id_to_ss_object_id
62- from lsst .sphgeom import ConvexPolygon
6360from lsst .utils .timer import timeMethod
6461
62+
6563class GenerateEphemeridesConnections (PipelineTaskConnections ,
66- dimensions = ("instrument" ,)):
64+ dimensions = ("instrument" ,)):
6765
6866 visitSummaries = connectionTypes .Input (
6967 doc = "Summary of visit information including ra, dec, and time" ,
7068 name = "preliminary_visit_summary" ,
7169 storageClass = "ExposureCatalog" ,
7270 dimensions = {"instrument" , "visit" },
73- multiple = True
71+ multiple = True
7472 )
7573 mpcorb = connectionTypes .Input (
7674 doc = "Minor Planet Center orbit table used for association" ,
@@ -142,19 +140,18 @@ class GenerateEphemeridesConnections(PipelineTaskConnections,
142140 dimensions = ("instrument" , "visit" ),
143141 multiple = True ,
144142 )
145-
146143
147144
148145class GenerateEphemeridesConfig (
149146 PipelineTaskConfig ,
150147 pipelineConnections = GenerateEphemeridesConnections ):
151- observatoryCode = pexConfig . Field (
148+ observatoryCode = Field (
152149 dtype = str ,
153150 doc = "IAU Minor Planet Center observer code for queries "
154151 "(default is X05 for Rubin Obs./LSST)" ,
155152 default = 'X05'
156153 )
157- observatoryFOVRadius = pexConfig . Field (
154+ observatoryFOVRadius = Field (
158155 dtype = float ,
159156 doc = "The field of view of the observatory (degrees)" ,
160157 default = 2.06 ,
@@ -202,7 +199,6 @@ def runQuantum(
202199 ephemeris_visit = outputs .ssObjects [dataId ['visit' ]]
203200 butlerQC .put (ephemeris_visit , ref )
204201
205-
206202 @timeMethod
207203 def run (self , visitSummaries , mpcorb , de440s , sb441_n16 , obsCodes , linux_p1550p2650 , pck00010 ,
208204 earth_latest_high_prec , earth_620120_250826 , earth_2025_250826_2125_predict , naif0012 ):
@@ -214,6 +210,11 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
214210 visitSummary : `lsst.afw.table.ExposureCatalog`
215211 Has rows with .getVisitInfo, which give the center and time of exposure
216212
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+
217218 Returns
218219 -------
219220 result : `lsst.pipe.base.Struct`
@@ -229,7 +230,7 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
229230 RA in decimal degrees (`float`)
230231 ``dec``
231232 DEC in decimal degrees (`float`)
232-
233+
233234 """
234235 if len (visitSummaries ) == 0 :
235236 raise NoWorkFound ("No visits input!" )
@@ -247,7 +248,7 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
247248 "fieldRA" : fieldRA ,
248249 "fieldDec" : fieldDec ,
249250 "observationId" : visits ,
250- "visitTime" : np .ones (n ) * visitTime ,
251+ "visitTime" : np .ones (n ) * visitTime ,
251252 "observationStartMJD" : epochMJD - (visitTime / 2 ) / 86400.0 ,
252253 "visitExposureTime" : visitTime ,
253254
@@ -258,38 +259,39 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
258259 "fiveSigmaDepth" : [- 1 ] * n ,
259260 "rotSkyPos" : [- 1 ] * n ,
260261 })
261- mpcorb = mpcorb .iloc [:10_000 ] # TODO: remove
262262
263263 inputOrbits = mpcorb [
264- ['packed_primary_provisional_designation' , 'q' , 'e' , 'i' , 'node' ,'argperi' ,'peri_time' ,'epoch_mjd' ]
264+ ['packed_primary_provisional_designation' , 'q' , 'e' , 'i' ,
265+ 'node' , 'argperi' , 'peri_time' , 'epoch_mjd' ]
265266 ].rename (columns = {'packed_primary_provisional_designation' : 'ObjID' , 'epoch_mjd' : 'epochMJD_TDB' ,
266267 'i' : 'inc' , 'argperi' : 'argPeri' , 'peri_time' : 't_p_MJD_TDB' })
267268 inputOrbits ['FORMAT' ] = 'COM'
268269 # Colors exactly like jake's prep_input_colors
269270 inputColors = inputOrbits [["ObjID" ]].copy ()
270271 inputColors ["H_r" ] = mpcorb ['h' ]
271- inputColors ["GS" ] = 0.15 # Traditional
272+ inputColors ["GS" ] = 0.15 # Traditional
272273
273274 eph_str = resources .files (lsst .pipe .tasks ).parents [3 ].joinpath ("data/eph.ini" ).read_text ()
274275 eph_str = eph_str .replace ("{OBSCODE}" , self .config .observatoryCode )
275276 eph_str = eph_str .replace ("{FOV}" , str (self .config .observatoryFOVRadius ))
276277
277- with tempfile . TemporaryDirectory () as tmpdirname :
278+ with TemporaryDirectory () as tmpdirname :
278279 self .log .info (f'temp dir: { tmpdirname } ' )
279-
280+
280281 # Orbits
281282 inputOrbits .to_csv (f'{ tmpdirname } /orbits.csv' , index = False )
282283 # Observations SQLite
283- conn = sqlite3 . connect (f'{ tmpdirname } /pointings.db' )
284+ conn = sqlite_connect (f'{ tmpdirname } /pointings.db' )
284285 inputVisits .to_sql ('observations' , conn , if_exists = 'replace' , index = False )
285286 conn .close ()
286287
287- open (f'{ tmpdirname } /eph.ini' , 'w' ).write (eph_str )
288+ with open (f'{ tmpdirname } /eph.ini' , 'w' ) as ephFile :
289+ ephFile .write (eph_str )
290+
288291 inputColors .to_csv (f'{ tmpdirname } /colors.csv' , index = False )
289292
290-
291293 cache = f'{ tmpdirname } /sorcha_cache/'
292- self .log .info (f 'making cache' )
294+ self .log .info ('making cache' )
293295 os .mkdir (cache )
294296 # DONE
295297 for filename , fileref in [
@@ -304,28 +306,32 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
304306 ('naif0012.tls' , naif0012 ),
305307 ]:
306308 self .log .info (f'writing { filename } ' )
307- open (cache + filename , 'wb' ). write ( fileref . fileContents . read ())
308- meta_kernel_text = f""" \\ begindata
309+ with open (cache + filename , 'wb' ) as file :
310+ file . write ( fileref . fileContents . read ())
309311
310- PATH_VALUES = ('{ tmpdirname } /sorcha_cache/')
312+ meta_kernel_text = dedent (f"""\
313+ \\ begindata
311314
312- PATH_SYMBOLS = ('A ')
315+ PATH_VALUES = ('{ tmpdirname } /sorcha_cache/ ')
313316
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- )
317+ PATH_SYMBOLS = ('A')
322318
323- \\ begintext
324- """
325- open (f'{ tmpdirname } /sorcha_cache/meta_kernel.txt' , 'w' ).write (meta_kernel_text )
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 )
326332 self .log .info ('Sorcha process begun' )
327333
328- result = subprocess . run (
334+ result = run (
329335 [
330336 "sorcha" ,
331337 "run" ,
@@ -337,18 +343,19 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
337343 "--ew" , f"{ tmpdirname } /ephemeris" ,
338344 "--ar" , f"{ tmpdirname } /sorcha_cache/"
339345 ],
340- stdout = subprocess . PIPE ,
341- stderr = subprocess . PIPE ,
346+ stdout = PIPE ,
347+ stderr = PIPE ,
342348 text = True
343349 )
344350
345351 self .log .info (f"Sorcha STDOUT:\n { result .stdout } " )
346352 self .log .info (f"Sorcha STDERR:\n { result .stderr } " )
347-
353+
348354 eph_path = f'{ tmpdirname } /ephemeris.csv'
349355 if not os .path .exists (eph_path ):
350356 raise FileNotFoundError (
351- f" Sorcha did not create ephemeris. Check STDOUT/STDERR above. Directory contents:\n { os .listdir (tmpdirname )} "
357+ " Sorcha did not create ephemeris. Check STDOUT/STDERR above. "
358+ f"Directory contents:\n { os .listdir (tmpdirname )} "
352359 )
353360
354361 # Return Sorcha output directly
@@ -358,7 +365,3 @@ def run(self, visitSummaries, mpcorb, de440s, sb441_n16, obsCodes, linux_p1550p2
358365 if v not in perVisitSsObjects :
359366 perVisitSsObjects [v ] = ephemeris .iloc [:0 ]
360367 return Struct (ssObjects = perVisitSsObjects )
361-
362-
363-
364-
0 commit comments