Coverage for python / lsst / pipe / tasks / ssoAssociation.py: 10%

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/>. 

21 

22"""Spatial association for Solar System Objects.""" 

23 

24__all__ = ["SolarSystemAssociationConfig", "SolarSystemAssociationTask"] 

25 

26from astropy import units as u 

27from astropy.coordinates import SkyCoord 

28from astropy.table import join, Column, Table 

29 

30import healpy as hp 

31import numpy as np 

32from numpy.polynomial.chebyshev import Chebyshev, chebval 

33from scipy.spatial import cKDTree 

34 

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 

40 

41from .ssp.ssobject import DIA_COLUMNS, DIA_DTYPES 

42 

43 

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 ) 

62 

63 

64class SolarSystemAssociationTask(pipeBase.Task): 

65 """Associate DIASources into existing SolarSystem Objects. 

66 

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" 

72 

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. 

77 

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 

92 

93 Returns 

94 ------- 

95 resultsStruct : `lsst.pipe.base.Struct` 

96 

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 """ 

107 

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) 

117 

118 if nSolarSystemObjects <= 0: 

119 return self._return_empty(diaSourceCatalog, ssObjects, source_column=source_column) 

120 

121 for c in ssObjects.columns: 

122 if ssObjects[c][0] is np.ma.core.MaskedConstant: 

123 ssObjects[c] = None 

124 

125 exposure_midpoint = visitInfo.date.toAstropy() 

126 if 'obs_x_poly' in ssObjects.columns: # mpSky ephemeris 

127 ref_time = exposure_midpoint.tai.mjd - ssObjects["tmin"].value[0] # all tmin should be identical 

128 ssObjects['obs_position'] = [ 

129 np.array([chebval(ref_time, row['obs_x_poly']), 

130 chebval(ref_time, row['obs_y_poly']), 

131 chebval(ref_time, row['obs_z_poly'])]) 

132 for row in ssObjects] 

133 ssObjects['obs_velocity'] = [ 

134 np.array([chebval(ref_time, Chebyshev(row['obs_x_poly']).deriv().coef), 

135 chebval(ref_time, Chebyshev(row['obs_y_poly']).deriv().coef), 

136 chebval(ref_time, Chebyshev(row['obs_z_poly']).deriv().coef)]) 

137 for row in ssObjects] 

138 ssObjects['obj_position'] = [ 

139 np.array([chebval(ref_time, row['obj_x_poly']), 

140 chebval(ref_time, row['obj_y_poly']), 

141 chebval(ref_time, row['obj_z_poly'])]) 

142 for row in ssObjects] 

143 ssObjects['obj_velocity'] = [ 

144 np.array([chebval(ref_time, Chebyshev(row['obj_x_poly']).deriv().coef), 

145 chebval(ref_time, Chebyshev(row['obj_y_poly']).deriv().coef), 

146 chebval(ref_time, Chebyshev(row['obj_z_poly']).deriv().coef)]) 

147 for row in ssObjects] 

148 vector = np.vstack(ssObjects['obj_position'].value - ssObjects['obs_position'].value) 

149 ras, decs = np.vstack(hp.vec2ang(vector, lonlat=True)) 

150 ssObjects['ephRa'] = ras 

151 ssObjects['ephDec'] = decs 

152 ssObjects['obs_position_x'], ssObjects['obs_position_y'], \ 

153 ssObjects['obs_position_z'] = ssObjects['obs_position'].value.T 

154 ssObjects['helio_x'], ssObjects['helio_y'], \ 

155 ssObjects['helio_z'] = ssObjects['obj_position'].value.T 

156 ssObjects['obs_velocity_x'], ssObjects['obs_velocity_y'], \ 

157 ssObjects['obs_velocity_z'] = ssObjects['obs_velocity'].value.T 

158 ssObjects['helio_vx'], ssObjects['helio_vy'], \ 

159 ssObjects['helio_vz'] = ssObjects['obj_velocity'].value.T 

160 ssObjects['topocentric_position'], ssObjects['topocentric_velocity'] = ( 

161 ssObjects['obj_position'] - ssObjects['obs_position'], 

162 ssObjects['obj_velocity'] - ssObjects['obs_velocity'], 

163 ) 

164 ssObjects['topo_x'], ssObjects['topo_y'], ssObjects['topo_z'] = ( 

165 np.array(list(ssObjects['topocentric_position'].value)).T 

166 ) 

167 ssObjects['topo_vx'], ssObjects['topo_vy'], ssObjects['topo_vz'] = ( 

168 np.array(list(ssObjects['topocentric_velocity'].value)).T 

169 ) 

170 ssObjects['helio_vx'], ssObjects['helio_vy'], \ 

171 ssObjects['helio_vz'] = np.array(list(ssObjects['obj_velocity'].value)).T 

172 ssObjects['helioRange'] = np.linalg.norm(ssObjects['obj_position'], axis=1) 

173 ssObjects['topoRange'] = np.linalg.norm(ssObjects['topocentric_position'], axis=1) 

174 ssObjects['phaseAngle'] = np.degrees(np.arccos(np.sum( 

175 ssObjects['obj_position'].T * ssObjects['topocentric_position'].T 

176 / ssObjects['helioRange'] / ssObjects['topoRange'], axis=0 

177 ))) 

178 # ssObjects['trailedSourceMagTrue'] should already be filled 

179 

180 # Add other required columns with dummy values until we compute them properly. 

181 # Fix in DM-53463 

182 ssObjects['RARateCosDec_deg_day'] = 0 

183 ssObjects['DecRate_deg_day'] = 0 

184 ssObjects['RangeRate_LTC_km_s'] = 0 

185 

186 marginArcsec = ssObjects["Err(arcsec)"].max() 

187 ssObjects['designation'] = ssObjects['MPCORB_unpacked_primary_provisional_designation'] 

188 

189 columns_to_drop = [ 

190 "obs_position", "obs_velocity", "obj_position", "obj_velocity", "topocentric_position", 

191 "topocentric_velocity", "obs_x_poly", "obs_y_poly", "obs_z_poly", "obj_x_poly", "obj_y_poly", 

192 "obj_z_poly", "associated" 

193 ] 

194 

195 else: # Sorcha ephemerides 

196 if 'trailedSourceMagTrue' not in ssObjects.columns: # Only possible for historical CI ephemerides 

197 ssObjects['trailedSourceMagTrue'] = 0 # Fix in DM-53462 

198 ssObjects.rename_columns( 

199 ['RATrue_deg', 'DecTrue_deg', 'phase_deg', 'Range_LTC_km', 'Obj_Sun_x_LTC_km', 

200 'Obj_Sun_y_LTC_km', 'Obj_Sun_z_LTC_km', 'Obj_Sun_vx_LTC_km_s', 'Obj_Sun_vy_LTC_km_s', 

201 'Obj_Sun_vz_LTC_km_s'], 

202 ['ephRa', 'ephDec', 'phaseAngle', 'topoRange', 'helio_x', 'helio_y', 

203 'helio_z', 'helio_vx', 'helio_vy', 'helio_vz']) 

204 ssObjects['ssObjectId'] = [obj_id_to_ss_object_id(v) for v in ssObjects['packed_desig']] 

205 ssObjects['helioRange'] = ( 

206 np.sqrt(ssObjects['helio_x']**2 + ssObjects['helio_y']**2 

207 + ssObjects['helio_z']**2) 

208 ) 

209 for substring1, substring2 in [('x', 'x_km'), ('y', 'y_km'), ('z', 'z_km'), 

210 ('vx', 'vx_km_s'), ('vy', 'vy_km_s'), ('vz', 'vz_km_s')]: 

211 topoName = 'topo_' + substring1 

212 helioName = 'helio_' + substring1 

213 obsName = 'Obs_Sun_' + substring2 

214 ssObjects[topoName] = ssObjects[helioName] - ssObjects[obsName] 

215 ssObjects['designation'] = ssObjects['ObjID'] 

216 

217 marginArcsec = 1.0 # TODO: justify 

218 

219 columns_to_drop = ['FieldID', 'fieldMJD_TAI', 'fieldJD_TDB', 

220 'Obs_Sun_x_km', 'Obs_Sun_y_km', 'Obs_Sun_z_km', 

221 'Obs_Sun_vx_km_s', 'Obs_Sun_vy_km_s', 'Obs_Sun_vz_km_s', 

222 '__index_level_0__'] 

223 

224 stateVectorColumns = ['helio_x', 'helio_y', 'helio_z', 'helio_vx', 

225 'helio_vy', 'helio_vz', 'topo_x', 'topo_y', 

226 'topo_z', 'topo_vx', 'topo_vy', 'topo_vz'] 

227 

228 mpcorbColumns = [col for col in ssObjects.columns if col[:7] == 'MPCORB_'] 

229 maskedObjects = self._maskToCcdRegion( 

230 ssObjects, 

231 bbox, 

232 wcs, 

233 marginArcsec).copy() 

234 nSolarSystemObjects = len(maskedObjects) 

235 if nSolarSystemObjects <= 0: 

236 return self._return_empty(diaSourceCatalog, maskedObjects, source_column=source_column) 

237 

238 maxRadius = np.deg2rad(self.config.maxDistArcSeconds / 3600) 

239 

240 # Transform DIA RADEC coordinates to unit sphere xyz for tree building. 

241 vectors = self._radec_to_xyz(diaSourceCatalog["ra"], 

242 diaSourceCatalog["dec"]) 

243 

244 # Create KDTree of DIA sources 

245 tree = cKDTree(vectors) 

246 

247 nFound = 0 

248 # Query the KDtree for DIA nearest neighbors to SSOs. Currently only 

249 # picks the DiaSource with the shortest distance. We can do something 

250 # fancier later. 

251 ssSourceData, ssObjectIds, prov_ids = [], [], [] 

252 ras, decs, residual_ras, residual_decs, dia_ids = [], [], [], [], [] 

253 diaSourceCatalog["ssObjectId"] = 0 

254 maskedObjects['associated'] = False 

255 

256 # Find all pairs of a source and an object within maxRadius 

257 nearby_obj_source_pairs = [] 

258 for obj_idx, (ra, dec) in enumerate(zip(maskedObjects["ephRa"].data, maskedObjects["ephDec"].data)): 

259 ssoVect = self._radec_to_xyz(ra, dec) 

260 dist, idx = tree.query(ssoVect, distance_upper_bound=maxRadius) 

261 if not np.isfinite(dist): 

262 continue 

263 for i in range(len(dist)): 

264 nearby_obj_source_pairs.append((dist, obj_idx, idx[i])) 

265 nearby_obj_source_pairs = sorted(nearby_obj_source_pairs) 

266 

267 # From closest to farthest, associate diaSources to SSOs. 

268 # Skipping already-associated sources and objects. 

269 all_cols = ( 

270 ["designation", "phaseAngle", "helioRange", "topoRange"] + stateVectorColumns + mpcorbColumns 

271 + ["ephRa", "ephDec", "RARateCosDec_deg_day", 

272 "DecRate_deg_day", "trailedSourceMagTrue", "RangeRate_LTC_km_s"] 

273 ) 

274 

275 used_src_indices, used_obj_indices = set(), set() 

276 maskedObjects['associated'] = False 

277 for dist, obj_idx, src_idx in nearby_obj_source_pairs: 

278 if src_idx in used_src_indices or obj_idx in used_obj_indices: 

279 continue 

280 maskedObject = maskedObjects[obj_idx] 

281 used_src_indices.add(src_idx) 

282 used_obj_indices.add(obj_idx) 

283 diaSourceCatalog[src_idx]["ssObjectId"] = maskedObject["ssObjectId"] 

284 ssObjectIds.append(maskedObject["ssObjectId"]) 

285 ssSourceData.append(list(maskedObject[all_cols].values())) 

286 dia_ra = diaSourceCatalog[src_idx]["ra"] 

287 dia_dec = diaSourceCatalog[src_idx]["dec"] 

288 dia_id = diaSourceCatalog[src_idx][source_column] 

289 ras.append(dia_ra) 

290 decs.append(dia_dec) 

291 dia_ids.append(dia_id) 

292 residual_ras.append(dia_ra - maskedObject["ephRa"]) 

293 residual_decs.append(dia_dec - maskedObject["ephDec"]) 

294 prov_ids.append(maskedObject['ObjID']) 

295 maskedObjects['associated'][obj_idx] = True 

296 nFound = len(ras) 

297 

298 self.log.info("Successfully associated %d / %d SolarSystemObjects.", nFound, nSolarSystemObjects) 

299 self.metadata['nAssociatedSsObjects'] = nFound 

300 self.metadata['nExpectedSsObjects'] = nSolarSystemObjects 

301 assocSourceMask = diaSourceCatalog["ssObjectId"] != 0 

302 unAssocObjectMask = np.logical_not(maskedObjects['associated'].value) 

303 dtypes = [type(d) if d is not np.ma.core.MaskedConstant else float 

304 for d in maskedObjects[0][all_cols].values()] 

305 ssSourceData = np.ma.filled(np.array(ssSourceData), None) 

306 colnames = ["designation", "phaseAngle", "helioRange", "topoRange"] 

307 colnames += stateVectorColumns + mpcorbColumns 

308 colnames += ["ephRa", "ephDec", "ephRateRa", "ephRateDec", "ephVmag", "topoRangeRate"] 

309 ssSourceData = Table(ssSourceData, names=colnames, dtype=dtypes) 

310 if 'MPCORB_created_at' in ssSourceData.columns: 

311 ssSourceData['MPCORB_created_at'] = 0 

312 ssSourceData['MPCORB_updated_at'] = 0 

313 ssSourceData['MPCORB_fitting_datetime'] = 0 

314 for c in ['MPCORB_created_at', 'MPCORB_updated_at', 'MPCORB_fitting_datetime', 

315 'MPCORB_u_param']: 

316 ssSourceData[c] = ssSourceData[c].astype(np.float64).astype(np.int64) 

317 ssSourceData['MPCORB_orbit_type_int'] = -1 # TODO: DM-54214. Mixed string and nan types. 

318 ssSourceData['MPCORB_orbit_type_int'] = ssSourceData['MPCORB_orbit_type_int'].astype(np.int64) 

319 ssSourceData['ssObjectId'] = Column(data=ssObjectIds, dtype=int) 

320 ssSourceData["ra"] = ras 

321 ssSourceData["dec"] = decs 

322 ephOffsetRa = np.array(residual_ras * np.cos(np.radians(ssSourceData["dec"]))) * 3600 # in arcsec 

323 ephOffsetDec = np.array(residual_decs) * 3600 # in arcsec 

324 ssSourceData["ephOffsetRa"] = ephOffsetRa 

325 ssSourceData["ephOffsetDec"] = ephOffsetDec 

326 ephOffsetVec = np.array([ephOffsetRa, ephOffsetDec]) 

327 ssSourceData[source_column] = dia_ids 

328 coords = SkyCoord(ra=ssSourceData['ra'].value * u.deg, dec=ssSourceData['dec'].value * u.deg) 

329 ssSourceData['galLon'] = coords.galactic.l.deg 

330 ssSourceData['galLat'] = coords.galactic.b.deg 

331 ssSourceData['eclLambda'] = coords.barycentrictrueecliptic.lon.deg 

332 ssSourceData['eclBeta'] = coords.barycentrictrueecliptic.lat.deg 

333 ssSourceData['ephRate'] = np.sqrt((ssSourceData['ephRateRa']) ** 2 

334 + (ssSourceData['ephRateDec']) ** 2) 

335 ssSourceData['ephOffset'] = np.sqrt((ssSourceData['ephOffsetRa']) ** 2 

336 + (ssSourceData['ephOffsetDec']) ** 2) 

337 ssSourceData['topo_vtot'] = np.sqrt(ssSourceData['topo_vx'] ** 2 

338 + ssSourceData['topo_vy'] ** 2 

339 + ssSourceData['topo_vy'] ** 2) 

340 ssSourceData['helio_vtot'] = np.sqrt(ssSourceData['helio_vx'] ** 2 

341 + ssSourceData['helio_vy'] ** 2 

342 + ssSourceData['helio_vy'] ** 2) 

343 skyMotionNormal0 = (ssSourceData['ephRateRa']/ssSourceData['ephRate']).data 

344 skyMotionNormal1 = (ssSourceData['ephRateDec']/ssSourceData['ephRate']).data 

345 skyMotionNormal = np.array([skyMotionNormal0, skyMotionNormal1]) 

346 skyMotionOrthogonal = np.array([[0, -1], [1, 0]]) @ skyMotionNormal 

347 ssSourceData['ephOffsetAlongTrack'] = (ephOffsetVec * skyMotionNormal).sum(axis=0) 

348 ssSourceData['ephOffsetCrossTrack'] = (ephOffsetVec * skyMotionOrthogonal).sum(axis=0) 

349 ssSourceData['diaDistanceRank'] = 1 

350 sun_obs_x = ssSourceData['topo_x'] - ssSourceData['helio_x'] 

351 sun_obs_y = ssSourceData['topo_y'] - ssSourceData['helio_y'] 

352 sun_obs_z = ssSourceData['topo_z'] - ssSourceData['helio_z'] 

353 sun_obs_range = np.sqrt(sun_obs_x**2 + sun_obs_y**2 + sun_obs_z**2) 

354 sun_obj_dot = (sun_obs_x * ssSourceData['topo_x'] + sun_obs_y * ssSourceData['topo_y'] 

355 + sun_obs_z * ssSourceData['topo_z']) 

356 ssSourceData["elongation"] = np.degrees(np.arccos(sun_obj_dot 

357 / (sun_obs_range * ssSourceData['topoRange']))) 

358 ssSourceData["helioRangeRate"] = ((ssSourceData["helio_vx"] * ssSourceData["helio_x"] 

359 + ssSourceData["helio_vy"] * ssSourceData["helio_y"] 

360 + ssSourceData["helio_vz"] * ssSourceData["helio_z"]) 

361 / ssSourceData["helioRange"]) 

362 

363 for distanceName in ['helio_x', 'helio_y', 'helio_z', 'topo_x', 'topo_y', 'topo_z', 

364 'helioRange', 'topoRange']: 

365 ssSourceData[distanceName] = ssSourceData[distanceName].astype(np.float64) 

366 ssSourceData[distanceName] *= (u.km).to(u.AU) 

367 unassociatedObjects = maskedObjects[unAssocObjectMask] 

368 unassociatedObjects.remove_columns(columns_to_drop) 

369 unassociatedObjects['ra'] = unassociatedObjects['ephRa'] 

370 unassociatedObjects['dec'] = unassociatedObjects['ephDec'] 

371 

372 # Add diaSource columns we care about when producing the SSObject table 

373 if len(ssSourceData): 

374 

375 # Extract only DIA_COLUMNS 

376 dia = diaSourceCatalog[source_columns].copy() 

377 

378 # Prefix all except diaSourceId 

379 for c in source_columns: 

380 if c != source_column: 

381 dia.rename_column(c, f"DIA_{c}") 

382 

383 # Join on diaSourceId, keeping all rows of ssSourceData 

384 ssSourceData = join(ssSourceData, dia, keys=source_column, join_type="left", uniq_col_name="") 

385 

386 return pipeBase.Struct( 

387 ssoAssocDiaSources=diaSourceCatalog[assocSourceMask], 

388 unAssocDiaSources=diaSourceCatalog[~assocSourceMask], 

389 nTotalSsObjects=nSolarSystemObjects, 

390 nAssociatedSsObjects=nFound, 

391 associatedSsSources=ssSourceData, 

392 unassociatedSsObjects=unassociatedObjects) 

393 

394 def _maskToCcdRegion(self, ssObjects, bbox, wcs, marginArcsec): 

395 """Mask the input SolarSystemObjects to only those in the exposure 

396 bounding box. 

397 

398 Parameters 

399 ---------- 

400 ssObjects : `astropy.table.Table` 

401 SolarSystemObjects to mask to ``exposure``. 

402 bbox : 

403 Exposure bbox used for masking 

404 wcs : 

405 Exposure wcs used for masking 

406 marginArcsec : `float` 

407 Maximum possible matching radius to pad onto the exposure bounding 

408 box. If greater than ``maxPixelMargin``, ``maxPixelMargin`` will 

409 be used. 

410 

411 Returns 

412 ------- 

413 maskedSolarSystemObjects : `astropy.table.Table` 

414 Set of SolarSystemObjects contained within the exposure bounds. 

415 """ 

416 if len(ssObjects) == 0: 

417 return ssObjects 

418 padding = min( 

419 int(np.ceil(marginArcsec / wcs.getPixelScale(bbox.getCenter()).asArcseconds())), 

420 self.config.maxPixelMargin) 

421 

422 return ssObjects[bbox_contains_sky_coords( 

423 bbox, 

424 wcs, 

425 ssObjects['ephRa'].value * u.degree, 

426 ssObjects['ephDec'].value * u.degree, 

427 padding)] 

428 

429 def _radec_to_xyz(self, ras, decs): 

430 """Convert input ra/dec coordinates to spherical unit-vectors. 

431 

432 Parameters 

433 ---------- 

434 ras : `array-like` 

435 RA coordinates of objects in degrees. 

436 decs : `array-like` 

437 DEC coordinates of objects in degrees. 

438 

439 Returns 

440 ------- 

441 vectors : `numpy.ndarray`, (N, 3) 

442 Output unit-vectors 

443 """ 

444 ras = np.radians(ras) 

445 decs = np.radians(decs) 

446 try: 

447 vectors = np.empty((len(ras), 3)) 

448 except TypeError: 

449 vectors = np.empty((1, 3)) 

450 

451 sin_dec = np.sin(np.pi / 2 - decs) 

452 vectors[:, 0] = sin_dec * np.cos(ras) 

453 vectors[:, 1] = sin_dec * np.sin(ras) 

454 vectors[:, 2] = np.cos(np.pi / 2 - decs) 

455 

456 return vectors 

457 

458 def _return_empty(self, diaSourceCatalog, emptySolarSystemObjects, source_column): 

459 """Return a struct with all appropriate empty values for no SSO associations. 

460 

461 Parameters 

462 ---------- 

463 diaSourceCatalog : `astropy.table.Table` 

464 Used for column names 

465 emptySolarSystemObjects : `astropy.table.Table` 

466 Used for column names. 

467 Returns 

468 ------- 

469 results : `lsst.pipe.base.Struct` 

470 Results struct with components. 

471 - ``ssoAssocDiaSources`` : Empty. (`astropy.table.Table`) 

472 - ``unAssocDiaSources`` : Input DiaSources. (`astropy.table.Table`) 

473 - ``nTotalSsObjects`` : Zero. (`int`) 

474 - ``nAssociatedSsObjects`` : Zero. 

475 - ``associatedSsSources`` : Empty. (`Astropy.table.Table`) 

476 - ``unassociatedSsObjects`` : Empty. (`Astropy.table.Table`) 

477 

478 

479 Raises 

480 ------ 

481 RuntimeError 

482 Raised if duplicate DiaObjects or duplicate DiaSources are found. 

483 """ 

484 self.log.info("No SolarSystemObjects found in detector bounding box.") 

485 # TODO DM-53699: these column names and dtypes should match the 

486 # definitions in sdm_schemas. 

487 dia_columns = [col for col in DIA_COLUMNS if not col == source_column] 

488 dia_dtypes = [d[0] for d in zip(DIA_DTYPES, DIA_COLUMNS) if not d[1] == source_column] 

489 names = ['designation', 'eclBeta', 'eclLambda', 'ephDec', 'ephOffsetDec', 'ephOffsetRa', 'ephRa', 

490 'galLat', 'galLon', 'elongation', 'ephOffset', 'ephOffsetAlongTrack', 'ephOffsetCrossTrack', 

491 'ephRate', 'ephRateDec', 'ephRateRa', 'ephVmag', 'helio_vtot', 'helio_vx', 'helio_vy', 

492 'helio_vz', 'helio_x', 'helio_y', 'helio_z', 'helioRange', 'helioRangeRate', 'phaseAngle', 

493 'topo_vtot', 'topo_vx', 'topo_vy', 'topo_vz', 'topo_x', 'topo_y', 'topo_z', 'topoRange', 

494 'topoRangeRate', source_column, 'ssObjectId', 'diaDistanceRank'] + dia_columns 

495 dtypes = [str] + [float] * 35 + [int] * 3 + dia_dtypes 

496 return pipeBase.Struct( 

497 ssoAssocDiaSources=Table(names=diaSourceCatalog.columns), 

498 unAssocDiaSources=diaSourceCatalog, 

499 nTotalSsObjects=0, 

500 nAssociatedSsObjects=0, 

501 associatedSsSources=Table(names=names, 

502 dtype=dtypes), 

503 unassociatedSsObjects=Table(names=emptySolarSystemObjects.columns) 

504 )