Coverage for python / lsst / pipe / tasks / ssoAssociation.py: 10%
211 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-26 09:20 +0000
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"""Spatial association for Solar System Objects."""
24__all__ = ["SolarSystemAssociationConfig", "SolarSystemAssociationTask"]
26from astropy import units as u
27from astropy.coordinates import SkyCoord
28from astropy.table import join, Column, Table
30import healpy as hp
31import numpy as np
32from numpy.polynomial.chebyshev import Chebyshev, chebval
33from scipy.spatial import cKDTree
35from lsst.afw.image.exposure.exposureUtils import bbox_contains_sky_coords
36import lsst.pex.config as pexConfig
37import lsst.pipe.base as pipeBase
38from lsst.utils.timer import timeMethod
39from lsst.pipe.tasks.associationUtils import obj_id_to_ss_object_id
41from .ssp.ssobject import DIA_COLUMNS, DIA_DTYPES
44class SolarSystemAssociationConfig(pexConfig.Config):
45 """Config class for SolarSystemAssociationTask.
46 """
47 maxDistArcSeconds = pexConfig.Field(
48 dtype=float,
49 doc='Maximum distance in arcseconds to test for a DIASource to be a '
50 'match to a SSObject.',
51 default=1.0,
52 )
53 maxPixelMargin = pexConfig.RangeField(
54 doc="Maximum padding to add to the ccd bounding box before masking "
55 "SolarSystem objects to the ccd footprint. The bounding box will "
56 "be padded by the minimum of this number or the max uncertainty "
57 "of the SolarSystemObjects in pixels.",
58 dtype=int,
59 default=100,
60 min=0,
61 )
64class SolarSystemAssociationTask(pipeBase.Task):
65 """Associate DIASources into existing SolarSystem Objects.
67 This task performs the association of detected DIASources in a visit
68 with known solar system objects.
69 """
70 ConfigClass = SolarSystemAssociationConfig
71 _DefaultName = "ssoAssociation"
73 @timeMethod
74 def run(self, diaSourceCatalog, ssObjects, visitInfo, bbox, wcs):
75 """Create a searchable tree of unassociated DiaSources and match
76 to the nearest ssoObject.
78 Parameters
79 ----------
80 diaSourceCatalog : `astropy.table.Table`
81 Catalog of DiaSources. Modified in place to add ssObjectId to
82 successfully associated DiaSources.
83 ssObjects : `astropy.table.Table`
84 Set of solar system objects that should be within the footprint
85 of the current visit.
86 visitInfo : `lsst.afw.image.VisitInfo`
87 visitInfo of exposure used for exposure time
88 bbox : `lsst.geom.Box2I`
89 bbox of exposure used for masking
90 wcs : `lsst.afw.geom.SkyWcs`
91 wcs of exposure used for masking
93 Returns
94 -------
95 resultsStruct : `lsst.pipe.base.Struct`
97 - ``ssoAssocDiaSources`` : DiaSources that were associated with
98 solar system objects in this visit. (`astropy.table.Table`)
99 - ``unAssocDiaSources`` : Set of DiaSources that were not
100 associated with any solar system object. (`astropy.table.Table`)
101 - ``nTotalSsObjects`` : Total number of SolarSystemObjects
102 contained in the CCD footprint. (`int`)
103 - ``nAssociatedSsObjects`` : Number of SolarSystemObjects
104 that were associated with DiaSources. (`int`)
105 - ``ssSourceData`` : ssSource table data. (`Astropy.table.Table`)
106 """
108 # TODO DM-53699: the source_column should be made consistent, and
109 # should match the definition in sdm_schemas.
110 source_columns = DIA_COLUMNS.copy()
111 if 'diaSourceId' in diaSourceCatalog.columns:
112 source_column = 'diaSourceId'
113 else:
114 source_column = 'id'
115 source_columns[source_columns.index('diaSourceId')] = source_column
116 nSolarSystemObjects = len(ssObjects)
118 if nSolarSystemObjects <= 0:
119 return self._return_empty(diaSourceCatalog, ssObjects, source_column=source_column)
121 for c in ssObjects.columns:
122 if ssObjects[c][0] is np.ma.core.MaskedConstant:
123 ssObjects[c] = None
125 exposure_midpoint = visitInfo.date.toAstropy()
126 if 'obs_x_poly' in ssObjects.columns: # mpSky ephemeris
127 # tmin, tmax are mjd. All tmin should be identical, so just take the first.
128 ref_time = exposure_midpoint.tai.mjd - ssObjects["tmin"].quantity.to_value(u.d)[0]
129 # *_poly are Chebyshev polynomials encoding the observer and object positions.
130 # Positions are received in equatorial, cartesian au. Velocities are in au/d.
131 ssObjects['obs_position'] = [
132 np.array([chebval(ref_time, row['obs_x_poly']),
133 chebval(ref_time, row['obs_y_poly']),
134 chebval(ref_time, row['obs_z_poly'])])
135 for row in ssObjects] * u.au
136 ssObjects['obs_velocity'] = [
137 np.array([chebval(ref_time, Chebyshev(row['obs_x_poly']).deriv().coef),
138 chebval(ref_time, Chebyshev(row['obs_y_poly']).deriv().coef),
139 chebval(ref_time, Chebyshev(row['obs_z_poly']).deriv().coef)])
140 for row in ssObjects] * u.au/u.d
141 ssObjects['obj_position'] = [
142 np.array([chebval(ref_time, row['obj_x_poly']),
143 chebval(ref_time, row['obj_y_poly']),
144 chebval(ref_time, row['obj_z_poly'])])
145 for row in ssObjects] * u.au
146 ssObjects['obj_velocity'] = [
147 np.array([chebval(ref_time, Chebyshev(row['obj_x_poly']).deriv().coef),
148 chebval(ref_time, Chebyshev(row['obj_y_poly']).deriv().coef),
149 chebval(ref_time, Chebyshev(row['obj_z_poly']).deriv().coef)])
150 for row in ssObjects] * u.au/u.d
152 vector = np.vstack(ssObjects['obj_position'].quantity.to_value(u.au)
153 - ssObjects['obs_position'].quantity.to_value(u.au))
154 ras, decs = np.vstack(hp.vec2ang(vector, lonlat=True))
155 # Angles in degrees
156 ssObjects['ephRa'] = ras * u.deg
157 ssObjects['ephDec'] = decs * u.deg
158 # Positions in au, velocities in km/s.
159 ssObjects['obs_position_x'], ssObjects['obs_position_y'], \
160 ssObjects['obs_position_z'] = ssObjects['obs_position'].quantity.to_value(u.au).T
161 ssObjects['helio_x'], ssObjects['helio_y'], \
162 ssObjects['helio_z'] = ssObjects['obj_position'].quantity.to_value(u.au).T
163 ssObjects['obs_velocity_x'], ssObjects['obs_velocity_y'], \
164 ssObjects['obs_velocity_z'] = ssObjects['obs_velocity'].quantity.to_value(u.km/u.s).T
165 ssObjects['helio_vx'], ssObjects['helio_vy'], \
166 ssObjects['helio_vz'] = ssObjects['obj_velocity'].quantity.to_value(u.km/u.s).T
167 ssObjects['topocentric_position'], ssObjects['topocentric_velocity'] = (
168 ssObjects['obj_position'] - ssObjects['obs_position'],
169 ssObjects['obj_velocity'] - ssObjects['obs_velocity'],
170 )
171 ssObjects['topo_x'], ssObjects['topo_y'], ssObjects['topo_z'] = (
172 np.array(list(ssObjects['topocentric_position'].quantity.to_value(u.au))).T
173 )
174 ssObjects['topo_vx'], ssObjects['topo_vy'], ssObjects['topo_vz'] = (
175 np.array(list(ssObjects['topocentric_velocity'].quantity.to_value(u.km/u.s))).T
176 )
177 ssObjects['helioRange'] = np.linalg.norm(ssObjects['obj_position'], axis=1)
178 ssObjects['topoRange'] = np.linalg.norm(ssObjects['topocentric_position'], axis=1)
180 # Phase angle in degrees
181 ssObjects['phaseAngle'] = np.degrees(np.arccos(np.sum(
182 ssObjects['obj_position'].T * ssObjects['topocentric_position'].T
183 / ssObjects['helioRange'] / ssObjects['topoRange'], axis=0
184 )))
186 # ssObjects['trailedSourceMagTrue'] should already be filled
187 # Add other required columns with dummy values until we compute them properly.
188 # Fix in DM-53463
189 ssObjects['RARateCosDec_deg_day'] = 0
190 ssObjects['DecRate_deg_day'] = 0
191 ssObjects['RangeRate_LTC_km_s'] = 0
193 marginArcsec = ssObjects["Err(arcsec)"].max()
194 ssObjects['designation'] = ssObjects['MPCORB_unpacked_primary_provisional_designation']
196 columns_to_drop = [
197 "obs_position", "obs_velocity", "obj_position", "obj_velocity", "topocentric_position",
198 "topocentric_velocity", "obs_x_poly", "obs_y_poly", "obs_z_poly", "obj_x_poly", "obj_y_poly",
199 "obj_z_poly", "associated"
200 ]
202 else: # Sorcha ephemerides
203 if 'trailedSourceMagTrue' not in ssObjects.columns: # Only possible for historical CI ephemerides
204 ssObjects['trailedSourceMagTrue'] = 0 # Fix in DM-53462
205 ssObjects.rename_columns(
206 ['RATrue_deg', 'DecTrue_deg', 'phase_deg', 'Range_LTC_km', 'Obj_Sun_x_LTC_km',
207 'Obj_Sun_y_LTC_km', 'Obj_Sun_z_LTC_km', 'Obj_Sun_vx_LTC_km_s', 'Obj_Sun_vy_LTC_km_s',
208 'Obj_Sun_vz_LTC_km_s'],
209 ['ephRa', 'ephDec', 'phaseAngle', 'topoRange', 'helio_x', 'helio_y',
210 'helio_z', 'helio_vx', 'helio_vy', 'helio_vz'])
211 for col in ['topoRange', 'helio_x', 'helio_y', 'helio_z']:
212 ssObjects[col] = ssObjects[col].quantity.to_value(u.au)
213 for col in ['Obs_Sun_x_km', 'Obs_Sun_y_km', 'Obs_Sun_z_km']:
214 ssObjects[col.replace('km', 'au')] = ssObjects[col].quantity.to_value(u.au)
215 for col in ['helio_vx', 'helio_vy', 'helio_vz', 'Obs_Sun_vx_km_s', 'Obs_Sun_vy_km_s',
216 'Obs_Sun_vz_km_s']:
217 ssObjects[col] = ssObjects[col].quantity.to_value(u.km/u.s)
219 ssObjects['ssObjectId'] = [obj_id_to_ss_object_id(v) for v in ssObjects['packed_desig']]
220 for substring1, substring2 in [('x', 'x_au'), ('y', 'y_au'), ('z', 'z_au'),
221 ('vx', 'vx_km_s'), ('vy', 'vy_km_s'), ('vz', 'vz_km_s')]:
222 topoName = 'topo_' + substring1
223 helioName = 'helio_' + substring1
224 obsName = 'Obs_Sun_' + substring2
225 ssObjects[topoName] = ssObjects[helioName] - ssObjects[obsName]
227 ssObjects['helioRange'] = (
228 np.sqrt(ssObjects['helio_x']**2 + ssObjects['helio_y']**2
229 + ssObjects['helio_z']**2)
230 )
232 ssObjects['designation'] = ssObjects['ObjID']
234 marginArcsec = 1.0 # TODO: justify
236 columns_to_drop = ['FieldID', 'fieldMJD_TAI', 'fieldJD_TDB',
237 'Obs_Sun_x_km', 'Obs_Sun_y_km', 'Obs_Sun_z_km',
238 'Obs_Sun_vx_km_s', 'Obs_Sun_vy_km_s', 'Obs_Sun_vz_km_s',
239 'Obs_Sun_x_au', 'Obs_Sun_y_au', 'Obs_Sun_z_au']
241 stateVectorColumns = ['helio_x', 'helio_y', 'helio_z', 'helio_vx',
242 'helio_vy', 'helio_vz', 'topo_x', 'topo_y',
243 'topo_z', 'topo_vx', 'topo_vy', 'topo_vz']
245 mpcorbColumns = [col for col in ssObjects.columns if col[:7] == 'MPCORB_']
246 maskedObjects = self._maskToCcdRegion(
247 ssObjects,
248 bbox,
249 wcs,
250 marginArcsec).copy()
251 nSolarSystemObjects = len(maskedObjects)
252 if nSolarSystemObjects <= 0:
253 return self._return_empty(diaSourceCatalog, maskedObjects, source_column=source_column)
255 maxRadius = np.deg2rad(self.config.maxDistArcSeconds / 3600)
257 # Transform DIA RADEC coordinates to unit sphere xyz for tree building.
258 vectors = self._radec_to_xyz(diaSourceCatalog["ra"],
259 diaSourceCatalog["dec"])
261 # Create KDTree of DIA sources
262 tree = cKDTree(vectors)
264 nFound = 0
265 # Query the KDtree for DIA nearest neighbors to SSOs. Currently only
266 # picks the DiaSource with the shortest distance. We can do something
267 # fancier later.
268 ssSourceData, ssObjectIds, prov_ids = [], [], []
269 ras, decs, residual_ras, residual_decs, dia_ids = [], [], [], [], []
270 diaSourceCatalog["ssObjectId"] = 0
271 maskedObjects['associated'] = False
273 # Find all pairs of a source and an object within maxRadius
274 nearby_obj_source_pairs = []
275 for obj_idx, (ra, dec) in enumerate(zip(maskedObjects["ephRa"].data, maskedObjects["ephDec"].data)):
276 ssoVect = self._radec_to_xyz(ra, dec)
277 dist, idx = tree.query(ssoVect, distance_upper_bound=maxRadius)
278 if not np.isfinite(dist):
279 continue
280 for i in range(len(dist)):
281 nearby_obj_source_pairs.append((dist, obj_idx, idx[i]))
282 nearby_obj_source_pairs = sorted(nearby_obj_source_pairs)
284 # From closest to farthest, associate diaSources to SSOs.
285 # Skipping already-associated sources and objects.
286 all_cols = (
287 ["designation", "phaseAngle", "helioRange", "topoRange"] + stateVectorColumns + mpcorbColumns
288 + ["ephRa", "ephDec", "RARateCosDec_deg_day",
289 "DecRate_deg_day", "trailedSourceMagTrue", "RangeRate_LTC_km_s"]
290 )
292 used_src_indices, used_obj_indices = set(), set()
293 maskedObjects['associated'] = False
294 for dist, obj_idx, src_idx in nearby_obj_source_pairs:
295 if src_idx in used_src_indices or obj_idx in used_obj_indices:
296 continue
297 maskedObject = maskedObjects[obj_idx]
298 used_src_indices.add(src_idx)
299 used_obj_indices.add(obj_idx)
300 diaSourceCatalog[src_idx]["ssObjectId"] = maskedObject["ssObjectId"]
301 ssObjectIds.append(maskedObject["ssObjectId"])
302 ssSourceData.append(list(maskedObject[all_cols].values()))
303 dia_ra = diaSourceCatalog[src_idx]["ra"]
304 dia_dec = diaSourceCatalog[src_idx]["dec"]
305 dia_id = diaSourceCatalog[src_idx][source_column]
306 ras.append(dia_ra)
307 decs.append(dia_dec)
308 dia_ids.append(dia_id)
309 residual_ras.append(dia_ra - maskedObject["ephRa"])
310 residual_decs.append(dia_dec - maskedObject["ephDec"])
311 prov_ids.append(maskedObject['ObjID'])
312 maskedObjects['associated'][obj_idx] = True
313 nFound = len(ras)
315 self.log.info("Successfully associated %d / %d SolarSystemObjects.", nFound, nSolarSystemObjects)
316 self.metadata['nAssociatedSsObjects'] = nFound
317 self.metadata['nExpectedSsObjects'] = nSolarSystemObjects
318 assocSourceMask = diaSourceCatalog["ssObjectId"] != 0
319 unAssocObjectMask = np.logical_not(maskedObjects['associated'].value)
320 dtypes = [float if d is np.ma.core.MaskedConstant else type(d)
321 for d in maskedObjects[0][all_cols].values()]
322 ssSourceData = np.ma.filled(np.array(ssSourceData), None)
323 colnames = ["designation", "phaseAngle", "helioRange", "topoRange"]
324 colnames += stateVectorColumns + mpcorbColumns
325 colnames += ["ephRa", "ephDec", "ephRateRa", "ephRateDec", "ephVmag", "topoRangeRate"]
326 ssSourceData = Table(ssSourceData, names=colnames, dtype=dtypes)
327 if 'MPCORB_created_at' in ssSourceData.columns:
328 ssSourceData['MPCORB_created_at'] = 0
329 ssSourceData['MPCORB_updated_at'] = 0
330 ssSourceData['MPCORB_fitting_datetime'] = 0
331 for c in ['MPCORB_created_at', 'MPCORB_updated_at', 'MPCORB_fitting_datetime',
332 'MPCORB_u_param']:
333 ssSourceData[c] = ssSourceData[c].astype(np.float64).astype(np.int64)
334 ssSourceData['MPCORB_orbit_type_int'] = -1 # TODO: DM-54214. Mixed string and nan types.
335 ssSourceData['MPCORB_orbit_type_int'] = ssSourceData['MPCORB_orbit_type_int'].astype(np.int64)
336 ssSourceData['ssObjectId'] = Column(data=ssObjectIds, dtype=int)
337 ssSourceData["ra"] = ras
338 ssSourceData["dec"] = decs
339 ephOffsetRa = np.array(residual_ras * np.cos(np.radians(ssSourceData["dec"]))) * 3600 # in arcsec
340 ephOffsetDec = np.array(residual_decs) * 3600 # in arcsec
341 ssSourceData["ephOffsetRa"] = ephOffsetRa
342 ssSourceData["ephOffsetDec"] = ephOffsetDec
343 ephOffsetVec = np.array([ephOffsetRa, ephOffsetDec])
344 ssSourceData[source_column] = dia_ids
345 coords = SkyCoord(ra=ssSourceData['ra'] * u.deg, dec=ssSourceData['dec'] * u.deg)
346 ssSourceData['galLon'] = coords.galactic.l.deg
347 ssSourceData['galLat'] = coords.galactic.b.deg
348 ssSourceData['eclLambda'] = coords.barycentrictrueecliptic.lon.deg
349 ssSourceData['eclBeta'] = coords.barycentrictrueecliptic.lat.deg
350 ssSourceData['ephRate'] = np.sqrt((ssSourceData['ephRateRa']) ** 2
351 + (ssSourceData['ephRateDec']) ** 2)
352 ssSourceData['ephOffset'] = np.sqrt((ssSourceData['ephOffsetRa']) ** 2
353 + (ssSourceData['ephOffsetDec']) ** 2)
354 ssSourceData['topo_vtot'] = np.sqrt(ssSourceData['topo_vx'] ** 2
355 + ssSourceData['topo_vy'] ** 2
356 + ssSourceData['topo_vz'] ** 2)
357 ssSourceData['helio_vtot'] = np.sqrt(ssSourceData['helio_vx'] ** 2
358 + ssSourceData['helio_vy'] ** 2
359 + ssSourceData['helio_vz'] ** 2)
360 skyMotionNormal0 = (ssSourceData['ephRateRa']/ssSourceData['ephRate']).data
361 skyMotionNormal1 = (ssSourceData['ephRateDec']/ssSourceData['ephRate']).data
362 skyMotionNormal = np.array([skyMotionNormal0, skyMotionNormal1])
363 skyMotionOrthogonal = np.array([[0, -1], [1, 0]]) @ skyMotionNormal
364 ssSourceData['ephOffsetAlongTrack'] = (ephOffsetVec * skyMotionNormal).sum(axis=0)
365 ssSourceData['ephOffsetCrossTrack'] = (ephOffsetVec * skyMotionOrthogonal).sum(axis=0)
366 ssSourceData['diaDistanceRank'] = 1
367 sun_obs_x = ssSourceData['topo_x'] - ssSourceData['helio_x']
368 sun_obs_y = ssSourceData['topo_y'] - ssSourceData['helio_y']
369 sun_obs_z = ssSourceData['topo_z'] - ssSourceData['helio_z']
370 sun_obs_range = np.sqrt(sun_obs_x**2 + sun_obs_y**2 + sun_obs_z**2)
371 sun_obj_dot = (sun_obs_x * ssSourceData['topo_x'] + sun_obs_y * ssSourceData['topo_y']
372 + sun_obs_z * ssSourceData['topo_z'])
373 ssSourceData["elongation"] = np.degrees(np.arccos(sun_obj_dot
374 / (sun_obs_range * ssSourceData['topoRange'])))
375 ssSourceData["helioRangeRate"] = ((ssSourceData["helio_vx"] * ssSourceData["helio_x"]
376 + ssSourceData["helio_vy"] * ssSourceData["helio_y"]
377 + ssSourceData["helio_vz"] * ssSourceData["helio_z"])
378 / ssSourceData["helioRange"])
380 unassociatedObjects = maskedObjects[unAssocObjectMask]
381 unassociatedObjects.remove_columns(columns_to_drop)
382 unassociatedObjects['ra'] = unassociatedObjects['ephRa']
383 unassociatedObjects['dec'] = unassociatedObjects['ephDec']
385 # Add diaSource columns we care about when producing the SSObject table
386 if len(ssSourceData):
388 # Extract only DIA_COLUMNS
389 dia = diaSourceCatalog[source_columns].copy()
391 # Prefix all except diaSourceId
392 for c in source_columns:
393 if c != source_column:
394 dia.rename_column(c, f"DIA_{c}")
396 # Join on diaSourceId, keeping all rows of ssSourceData
397 ssSourceData = join(ssSourceData, dia, keys=source_column, join_type="left", uniq_col_name="")
399 return pipeBase.Struct(
400 ssoAssocDiaSources=diaSourceCatalog[assocSourceMask],
401 unAssocDiaSources=diaSourceCatalog[~assocSourceMask],
402 nTotalSsObjects=nSolarSystemObjects,
403 nAssociatedSsObjects=nFound,
404 associatedSsSources=ssSourceData,
405 unassociatedSsObjects=unassociatedObjects)
407 def _maskToCcdRegion(self, ssObjects, bbox, wcs, marginArcsec):
408 """Mask the input SolarSystemObjects to only those in the exposure
409 bounding box.
411 Parameters
412 ----------
413 ssObjects : `astropy.table.Table`
414 SolarSystemObjects to mask to ``exposure``.
415 bbox :
416 Exposure bbox used for masking
417 wcs :
418 Exposure wcs used for masking
419 marginArcsec : `float`
420 Maximum possible matching radius to pad onto the exposure bounding
421 box. If greater than ``maxPixelMargin``, ``maxPixelMargin`` will
422 be used.
424 Returns
425 -------
426 maskedSolarSystemObjects : `astropy.table.Table`
427 Set of SolarSystemObjects contained within the exposure bounds.
428 """
429 if len(ssObjects) == 0:
430 return ssObjects
431 padding = min(
432 int(np.ceil(marginArcsec / wcs.getPixelScale(bbox.getCenter()).asArcseconds())),
433 self.config.maxPixelMargin)
435 return ssObjects[bbox_contains_sky_coords(
436 bbox,
437 wcs,
438 ssObjects['ephRa'],
439 ssObjects['ephDec'],
440 padding)]
442 def _radec_to_xyz(self, ras, decs):
443 """Convert input ra/dec coordinates to spherical unit-vectors.
445 Parameters
446 ----------
447 ras : `array-like`
448 RA coordinates of objects in degrees.
449 decs : `array-like`
450 DEC coordinates of objects in degrees.
452 Returns
453 -------
454 vectors : `numpy.ndarray`, (N, 3)
455 Output unit-vectors
456 """
457 ras = np.radians(ras)
458 decs = np.radians(decs)
459 try:
460 vectors = np.empty((len(ras), 3))
461 except TypeError:
462 vectors = np.empty((1, 3))
464 sin_dec = np.sin(np.pi / 2 - decs)
465 vectors[:, 0] = sin_dec * np.cos(ras)
466 vectors[:, 1] = sin_dec * np.sin(ras)
467 vectors[:, 2] = np.cos(np.pi / 2 - decs)
469 return vectors
471 def _return_empty(self, diaSourceCatalog, emptySolarSystemObjects, source_column):
472 """Return a struct with all appropriate empty values for no SSO associations.
474 Parameters
475 ----------
476 diaSourceCatalog : `astropy.table.Table`
477 Used for column names
478 emptySolarSystemObjects : `astropy.table.Table`
479 Used for column names.
480 Returns
481 -------
482 results : `lsst.pipe.base.Struct`
483 Results struct with components.
484 - ``ssoAssocDiaSources`` : Empty. (`astropy.table.Table`)
485 - ``unAssocDiaSources`` : Input DiaSources. (`astropy.table.Table`)
486 - ``nTotalSsObjects`` : Zero. (`int`)
487 - ``nAssociatedSsObjects`` : Zero.
488 - ``associatedSsSources`` : Empty. (`Astropy.table.Table`)
489 - ``unassociatedSsObjects`` : Empty. (`Astropy.table.Table`)
492 Raises
493 ------
494 RuntimeError
495 Raised if duplicate DiaObjects or duplicate DiaSources are found.
496 """
497 self.log.info("No SolarSystemObjects found in detector bounding box.")
498 # TODO DM-53699: these column names and dtypes should match the
499 # definitions in sdm_schemas.
500 dia_columns = [col for col in DIA_COLUMNS if not col == source_column]
501 dia_dtypes = [d[0] for d in zip(DIA_DTYPES, DIA_COLUMNS) if not d[1] == source_column]
502 names = ['designation', 'eclBeta', 'eclLambda', 'ephDec', 'ephOffsetDec', 'ephOffsetRa', 'ephRa',
503 'galLat', 'galLon', 'elongation', 'ephOffset', 'ephOffsetAlongTrack', 'ephOffsetCrossTrack',
504 'ephRate', 'ephRateDec', 'ephRateRa', 'ephVmag', 'helio_vtot', 'helio_vx', 'helio_vy',
505 'helio_vz', 'helio_x', 'helio_y', 'helio_z', 'helioRange', 'helioRangeRate', 'phaseAngle',
506 'topo_vtot', 'topo_vx', 'topo_vy', 'topo_vz', 'topo_x', 'topo_y', 'topo_z', 'topoRange',
507 'topoRangeRate', source_column, 'ssObjectId', 'diaDistanceRank'] + dia_columns
508 dtypes = [str] + [float] * 35 + [int] * 3 + dia_dtypes
509 return pipeBase.Struct(
510 ssoAssocDiaSources=Table(names=diaSourceCatalog.columns),
511 unAssocDiaSources=diaSourceCatalog,
512 nTotalSsObjects=0,
513 nAssociatedSsObjects=0,
514 associatedSsSources=Table(names=names,
515 dtype=dtypes),
516 unassociatedSsObjects=Table(names=emptySolarSystemObjects.columns)
517 )