Coverage for python/lsst/analysis/ap/templateInputs.py: 81%

245 statements  

« prev     ^ index     » next       coverage.py v7.16.2, created at 2026-10-08 00:05 +0000

1# This file is part of analysis_ap. 

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"""Find the images that make up the template at the position of a diaSource. 

23 

24A template written as an `lsst.images.DifferenceImage` records the coadds it 

25was built from in its ``templates`` component. Each coadd in turn records 

26its input images and their weights, either as a legacy 

27`lsst.afw.image.CoaddInputs` (``coaddInputs`` component) or as an 

28`lsst.images.cells.CoaddProvenance` (``provenance`` component). The tools 

29here follow that chain from a diaSource back to the observation times and 

30fractional weights of the template's input images. 

31""" 

32 

33__all__ = ["get_template_inputs", "summarize_template_inputs"] 

34 

35import logging 

36 

37import astropy.table 

38import astropy.units as u 

39import numpy as np 

40import pandas as pd 

41from astropy.coordinates import SkyCoord 

42 

43import lsst.geom 

44from lsst.daf.base import DateTime 

45from lsst.images import Box, Interval 

46 

47_LOG = logging.getLogger(__name__) 

48 

49_INPUT_COLUMNS = ["input_visit", "input_detector", "coadd_weight", "mid_mjd_tai", "exposure_time"] 

50"""Columns of the per-coadd input tables returned by the readers 

51(`list` [`str`]). 

52""" 

53 

54 

55def _empty_inputs(): 

56 """Return a per-coadd input table with no rows.""" 

57 return astropy.table.Table( 

58 names=["point"] + _INPUT_COLUMNS, 

59 dtype=[np.int64, np.int64, np.int64, np.float64, np.float64, np.float64], 

60 ) 

61 

62 

63def _nearest_pixel(value): 

64 """Return the integer index of the pixel containing a position.""" 

65 return np.floor(np.asarray(value) + 0.5).astype(np.int64) 

66 

67 

68def _clipped_box(ix, iy, bbox): 

69 """Return the box enclosing pixel indices, clipped to a bounding box. 

70 

71 Parameters 

72 ---------- 

73 ix, iy : `numpy.ndarray` 

74 Integer pixel indices to enclose. 

75 bbox : `lsst.images.Box` 

76 Bounding box of the image to read from. 

77 

78 Returns 

79 ------- 

80 box : `lsst.images.Box` or `None` 

81 Enclosing box, or `None` if no pixel is inside ``bbox``. 

82 """ 

83 inside = bbox.contains(x=ix, y=iy) 

84 if not inside.any(): 84 ↛ 85line 84 didn't jump to line 85 because the condition on line 84 was never true

85 return None 

86 return Box( 

87 y=Interval(iy[inside].min(), iy[inside].max() + 1), 

88 x=Interval(ix[inside].min(), ix[inside].max() + 1), 

89 ) 

90 

91 

92class _LegacyCoaddReader: 

93 """Inputs of a coadd whose provenance is a legacy 

94 `lsst.afw.image.CoaddInputs`. 

95 

96 Parameters 

97 ---------- 

98 ccds : `lsst.afw.table.ExposureCatalog` 

99 The ``ccds`` catalog of the coadd's ``coaddInputs``. 

100 butler : `lsst.daf.butler.Butler`, optional 

101 Butler to read the coadd's variance plane from. 

102 ref : `lsst.daf.butler.DatasetRef`, optional 

103 Reference to the coadd. 

104 """ 

105 

106 def __init__(self, ccds, butler=None, ref=None): 

107 self._ccds = ccds 

108 self._butler = butler 

109 self._ref = ref 

110 self._mid_mjd = np.full(len(ccds), np.nan) 

111 self._exposure_time = np.full(len(ccds), np.nan) 

112 for i, record in enumerate(ccds): 

113 visit_info = record.getVisitInfo() 

114 if visit_info is not None: 114 ↛ 112line 114 didn't jump to line 112 because the condition on line 114 was always true

115 self._mid_mjd[i] = visit_info.getDate().get(DateTime.MJD, DateTime.TAI) 

116 self._exposure_time[i] = visit_info.getExposureTime() 

117 

118 @classmethod 

119 def from_butler(cls, butler, ref): 

120 """Read the provenance of a coadd. 

121 

122 Parameters 

123 ---------- 

124 butler : `lsst.daf.butler.Butler` 

125 Butler to read from. 

126 ref : `lsst.daf.butler.DatasetRef` 

127 Reference to the coadd. 

128 

129 Returns 

130 ------- 

131 reader : `_LegacyCoaddReader` 

132 Reader for that coadd. 

133 """ 

134 coadd_inputs = butler.get(ref.makeComponentRef("coaddInputs")) 

135 return cls(coadd_inputs.ccds, butler=butler, ref=ref) 

136 

137 def find_inputs(self, ra, dec): 

138 """Find the input images that cover sky positions. 

139 

140 Parameters 

141 ---------- 

142 ra, dec : `numpy.ndarray` 

143 Sky positions, in degrees. 

144 

145 Returns 

146 ------- 

147 inputs : `astropy.table.Table` 

148 One row per (position, input image), with ``point`` giving the 

149 index of the position. 

150 

151 Notes 

152 ----- 

153 An input covers a position if the position is inside both the bounding 

154 box and the valid polygon of that input, the same rule 

155 `lsst.meas.algorithms.CoaddPsf` uses. 

156 """ 

157 points = [lsst.geom.SpherePoint(r, d, lsst.geom.degrees) for r, d in zip(ra, dec)] 

158 rows = [] 

159 for i, record in enumerate(self._ccds): 

160 bbox = lsst.geom.Box2D(record.getBBox()) 

161 polygon = record.getValidPolygon() 

162 for point, pixel in enumerate(record.getWcs().skyToPixel(points)): 

163 if bbox.contains(pixel) and (polygon is None or polygon.contains(pixel)): 

164 rows.append((point, record["visit"], record["ccd"], record["weight"], 

165 self._mid_mjd[i], self._exposure_time[i])) 

166 if not rows: 166 ↛ 167line 166 didn't jump to line 167 because the condition on line 166 was never true

167 return _empty_inputs() 

168 return astropy.table.Table(rows=rows, names=_empty_inputs().colnames) 

169 

170 def variance_at(self, ra, dec): 

171 """Read the coadd variance at sky positions. 

172 

173 Parameters 

174 ---------- 

175 ra, dec : `numpy.ndarray` 

176 Sky positions, in degrees. 

177 

178 Returns 

179 ------- 

180 variance : `numpy.ndarray` 

181 Variance of the nearest coadd pixel; ``nan`` outside the coadd. 

182 """ 

183 wcs = self._butler.get(self._ref.makeComponentRef("wcs")) 

184 points = [lsst.geom.SpherePoint(r, d, lsst.geom.degrees) for r, d in zip(ra, dec)] 

185 pixels = wcs.skyToPixel(points) 

186 ix = _nearest_pixel([p.getX() for p in pixels]) 

187 iy = _nearest_pixel([p.getY() for p in pixels]) 

188 coadd_bbox = Box.from_legacy(self._butler.get(self._ref.makeComponentRef("bbox"))) 

189 box = _clipped_box(ix, iy, coadd_bbox) 

190 result = np.full(len(ix), np.nan) 

191 if box is None: 191 ↛ 192line 191 didn't jump to line 192 because the condition on line 191 was never true

192 return result 

193 variance = self._butler.get(self._ref.makeComponentRef("variance"), 

194 parameters={"bbox": box.to_legacy()}) 

195 inside = box.contains(x=ix, y=iy) 

196 result[inside] = variance.array[iy[inside] - box.y.start, ix[inside] - box.x.start] 

197 return result 

198 

199 

200class _CellCoaddReader: 

201 """Inputs of a coadd whose provenance is an 

202 `lsst.images.cells.CoaddProvenance`. 

203 

204 Parameters 

205 ---------- 

206 provenance : `lsst.images.cells.CoaddProvenance` 

207 Provenance of the coadd. 

208 sky_projection : `lsst.images.SkyProjection` 

209 Mapping from the coadd's pixels to the sky. 

210 grid : `lsst.images.CellGrid` 

211 Cell grid of the coadd's patch. 

212 times : `dict` [`tuple` [`str`, `int`], `tuple` [`float`, `float`]] 

213 Mid-exposure time (MJD, TAI) and exposure time (s), keyed by 

214 ``(instrument, visit)``. 

215 butler : `lsst.daf.butler.Butler`, optional 

216 Butler to read the coadd's variance plane from. 

217 ref : `lsst.daf.butler.DatasetRef`, optional 

218 Reference to the coadd. 

219 """ 

220 

221 def __init__(self, provenance, sky_projection, grid, times, butler=None, ref=None): 

222 self._sky_projection = sky_projection 

223 self._grid = grid 

224 self._butler = butler 

225 self._ref = ref 

226 self._polygons = { 

227 (str(row["instrument"]), int(row["visit"]), int(row["detector"])): row["polygon"] 

228 for row in provenance.inputs 

229 } 

230 self._contributions = provenance.contributions 

231 self._rows_by_cell = {} 

232 for n, row in enumerate(self._contributions): 

233 self._rows_by_cell.setdefault((int(row["cell_i"]), int(row["cell_j"])), []).append(n) 

234 self._times = times 

235 

236 @classmethod 

237 def from_butler(cls, butler, ref): 

238 """Read the provenance of a coadd and the times of its inputs. 

239 

240 Parameters 

241 ---------- 

242 butler : `lsst.daf.butler.Butler` 

243 Butler to read from. 

244 ref : `lsst.daf.butler.DatasetRef` 

245 Reference to the coadd. 

246 

247 Returns 

248 ------- 

249 reader : `_CellCoaddReader` 

250 Reader for that coadd. 

251 """ 

252 parts = butler.get(ref.makeComponentRef("components"), 

253 parameters={"components": ["provenance", "sky_projection", "psf"]}) 

254 provenance = parts["provenance"] 

255 times = {} 

256 for instrument in np.unique(provenance.inputs["instrument"]): 

257 visits = provenance.inputs["visit"][provenance.inputs["instrument"] == instrument] 

258 records = butler.query_dimension_records( 

259 "visit", instrument=str(instrument), where="visit IN (:visits)", 

260 bind={"visits": [int(v) for v in np.unique(visits)]}, limit=None, 

261 ) 

262 for record in records: 

263 timespan = record.timespan 

264 mid = timespan.begin + (timespan.end - timespan.begin)/2 

265 times[str(instrument), record.id] = (mid.tai.mjd, record.exposure_time) 

266 return cls(provenance, parts["sky_projection"], parts["psf"].bounds.grid, times, 

267 butler=butler, ref=ref) 

268 

269 def _pixels(self, ra, dec): 

270 """Return the coadd pixel positions of sky positions.""" 

271 xy = self._sky_projection.sky_to_pixel(SkyCoord(ra*u.deg, dec*u.deg)) 

272 return np.atleast_1d(xy.x), np.atleast_1d(xy.y) 

273 

274 def find_inputs(self, ra, dec): 

275 """Find the input images that cover sky positions. 

276 

277 Parameters 

278 ---------- 

279 ra, dec : `numpy.ndarray` 

280 Sky positions, in degrees. 

281 

282 Returns 

283 ------- 

284 inputs : `astropy.table.Table` 

285 One row per (position, input image), with ``point`` giving the 

286 index of the position. 

287 

288 Notes 

289 ----- 

290 An input covers a position if it contributes to the cell containing 

291 the position and its overlap polygon contains the position. The 

292 weight is the one recorded for that cell. 

293 """ 

294 x, y = self._pixels(ra, dec) 

295 rows = [] 

296 for point, (px, py) in enumerate(zip(x, y)): 

297 cell = self._grid.index_of(y=int(_nearest_pixel(py)), x=int(_nearest_pixel(px))) 

298 for n in self._rows_by_cell.get((cell.i, cell.j), []): 

299 row = self._contributions[n] 

300 instrument = str(row["instrument"]) 

301 key = (instrument, int(row["visit"]), int(row["detector"])) 

302 if not self._polygons[key].contains(x=px, y=py): 

303 continue 

304 mid_mjd, exposure_time = self._times.get(key[:2], (np.nan, np.nan)) 

305 rows.append((point, key[1], key[2], row["weight"], mid_mjd, exposure_time)) 

306 if not rows: 306 ↛ 307line 306 didn't jump to line 307 because the condition on line 306 was never true

307 return _empty_inputs() 

308 return astropy.table.Table(rows=rows, names=_empty_inputs().colnames) 

309 

310 def variance_at(self, ra, dec): 

311 """Read the coadd variance at sky positions. 

312 

313 Parameters 

314 ---------- 

315 ra, dec : `numpy.ndarray` 

316 Sky positions, in degrees. 

317 

318 Returns 

319 ------- 

320 variance : `numpy.ndarray` 

321 Variance of the nearest coadd pixel; ``nan`` outside the coadd. 

322 """ 

323 x, y = self._pixels(ra, dec) 

324 ix = _nearest_pixel(x) 

325 iy = _nearest_pixel(y) 

326 box = _clipped_box(ix, iy, self._grid.bbox) 

327 result = np.full(len(ix), np.nan) 

328 if box is None: 

329 return result 

330 variance = self._butler.get(self._ref.makeComponentRef("variance"), parameters={"bbox": box}) 

331 inside = box.contains(x=ix, y=iy) 

332 result[inside] = variance.array[iy[inside] - box.y.start, ix[inside] - box.x.start] 

333 return result 

334 

335 

336_READERS = {"coaddInputs": _LegacyCoaddReader, "provenance": _CellCoaddReader} 

337"""Coadd reader for each provenance component, in order of preference 

338(`dict` [`str`, `type`]). 

339""" 

340 

341 

342def _get_reader(butler, info, readers): 

343 """Return the reader for one template coadd, reading it if needed. 

344 

345 Parameters 

346 ---------- 

347 butler : `lsst.daf.butler.Butler` 

348 Butler to read from. 

349 info : `lsst.images.DifferenceImageTemplateInfo` 

350 Record of the coadd. 

351 readers : `dict` 

352 Readers already made, keyed by dataset ID; updated in place. 

353 

354 Returns 

355 ------- 

356 reader : `_LegacyCoaddReader` or `_CellCoaddReader` 

357 Reader for the coadd. 

358 

359 Raises 

360 ------ 

361 LookupError 

362 Raised if the butler has no dataset with the coadd's ID. 

363 TypeError 

364 Raised if the coadd's storage class has no known provenance 

365 component. 

366 """ 

367 if info.dataset_id in readers: 367 ↛ 368line 367 didn't jump to line 368 because the condition on line 367 was never true

368 return readers[info.dataset_id] 

369 ref = butler.get_dataset(info.dataset_id) 

370 if ref is None: 370 ↛ 371line 370 didn't jump to line 371 because the condition on line 370 was never true

371 raise LookupError(f"Template coadd tract={info.tract}, patch={info.patch} " 

372 f"(id={info.dataset_id}, run={info.dataset_run!r}) is not in this butler.") 

373 components = ref.datasetType.storageClass.allComponents() 

374 for component, reader_type in _READERS.items(): 374 ↛ 378line 374 didn't jump to line 378 because the loop on line 374 didn't complete

375 if component in components: 375 ↛ 374line 375 didn't jump to line 374 because the condition on line 375 was always true

376 readers[info.dataset_id] = reader_type.from_butler(butler, ref) 

377 return readers[info.dataset_id] 

378 raise TypeError(f"Coadd {ref} has storage class {ref.datasetType.storageClass_name}, which has none " 

379 f"of the provenance components {list(_READERS)}.") 

380 

381 

382def _combine_weights(coadds, variances): 

383 """Compute the fraction of the template flux at one position that came 

384 from each input image. 

385 

386 Parameters 

387 ---------- 

388 coadds : `list` [`tuple`] 

389 One ``(key, tract, inputs)`` entry per coadd whose template bounds 

390 contain the position, where ``inputs`` is that coadd's input table 

391 for the position. 

392 variances : `dict` 

393 Coadd variance at the position, keyed by ``key``. Needed only when 

394 there is more than one coadd. 

395 

396 Returns 

397 ------- 

398 fractions : `dict` 

399 Weight to apply to each coadd's ``coadd_weight`` column, keyed by 

400 ``key``. Coadds that contribute nothing are omitted. 

401 

402 Notes 

403 ----- 

404 This follows `lsst.ip.diffim.GetTemplateTask`: the patches of each tract 

405 are averaged with per-pixel weights ``variance**-0.5``, and the tracts 

406 are then averaged the same way using the averaged variance of each tract. 

407 Within one coadd, the fraction from an input is its coadd weight over the 

408 sum of the weights of the inputs covering the position. 

409 """ 

410 coadds = [(key, tract, inputs) for key, tract, inputs in coadds if inputs["coadd_weight"].sum() > 0] 

411 if len(coadds) == 1: 

412 key, _, inputs = coadds[0] 

413 return {key: 1.0/inputs["coadd_weight"].sum()} 

414 patch_weights = {} 

415 tract_variances = {} 

416 for key, tract, inputs in coadds: 

417 variance = variances.get(key, np.nan) 

418 if np.isfinite(variance) and variance > 0: 

419 patch_weights.setdefault(tract, {})[key] = variance**-0.5 

420 for tract, weights in patch_weights.items(): 

421 total = sum(weights.values()) 

422 for key in weights: 

423 weights[key] /= total 

424 tract_variances[tract] = sum(weights[key]*variances[key] for key in weights) 

425 tract_weights = {tract: variance**-0.5 for tract, variance in tract_variances.items()} 

426 tract_total = sum(tract_weights.values()) 

427 fractions = {} 

428 for key, tract, inputs in coadds: 

429 if key in patch_weights.get(tract, {}): 

430 fractions[key] = (tract_weights[tract]/tract_total*patch_weights[tract][key] 

431 / inputs["coadd_weight"].sum()) 

432 return fractions 

433 

434 

435def _template_inputs_for_detector(butler, data_id, sources, template_dataset, readers): 

436 """Find the template inputs for the diaSources on one detector. 

437 

438 Parameters 

439 ---------- 

440 butler : `lsst.daf.butler.Butler` 

441 Butler to read from. 

442 data_id : `dict` 

443 Data ID of the detector. 

444 sources : `astropy.table.Table` 

445 DiaSources on the detector. 

446 template_dataset : `str` 

447 Name of the template dataset type. 

448 readers : `dict` 

449 Coadd readers already made, keyed by dataset ID; updated in place. 

450 

451 Returns 

452 ------- 

453 inputs : `list` [`astropy.table.Table`] 

454 Input tables, one per contributing coadd. 

455 """ 

456 templates = butler.get(f"{template_dataset}.templates", data_id) 

457 ra = np.asarray(sources["ra"], dtype=float) 

458 dec = np.asarray(sources["dec"], dtype=float) 

459 x = np.asarray(sources["x"], dtype=float) 

460 y = np.asarray(sources["y"], dtype=float) 

461 per_coadd = [] 

462 for info in templates: 

463 contains = np.atleast_1d(info.bounds.contains(x=x, y=y)) 

464 if not contains.any(): 464 ↛ 465line 464 didn't jump to line 465 because the condition on line 464 was never true

465 continue 

466 indices = np.flatnonzero(contains) 

467 reader = _get_reader(butler, info, readers) 

468 inputs = reader.find_inputs(ra[indices], dec[indices]) 

469 inputs["point"] = indices[inputs["point"]] 

470 per_coadd.append((info, reader, indices, inputs)) 

471 n_coadds = np.zeros(len(sources), dtype=int) 

472 for _, _, indices, _ in per_coadd: 

473 n_coadds[indices] += 1 

474 variances = {} 

475 for info, reader, indices, _ in per_coadd: 

476 overlap = indices[n_coadds[indices] > 1] 

477 if len(overlap): 477 ↛ 475line 477 didn't jump to line 475 because the condition on line 477 was always true

478 for point, variance in zip(overlap, reader.variance_at(ra[overlap], dec[overlap])): 

479 variances[point, info.dataset_id] = variance 

480 results = [] 

481 for point in range(len(sources)): 

482 coadds = [] 

483 infos = {} 

484 for info, _, _, inputs in per_coadd: 

485 mine = inputs[inputs["point"] == point] 

486 if len(mine): 

487 coadds.append(((point, info.dataset_id), info.tract, mine)) 

488 infos[point, info.dataset_id] = info 

489 fractions = _combine_weights(coadds, variances) if coadds else {} 

490 if not fractions: 490 ↛ 491line 490 didn't jump to line 491 because the condition on line 490 was never true

491 _LOG.warning("No template coadd input contributes to diaSource %s.", 

492 sources["diaSourceId"][point]) 

493 continue 

494 for key, _, inputs in coadds: 

495 if key not in fractions: 495 ↛ 496line 495 didn't jump to line 496 because the condition on line 495 was never true

496 continue 

497 table = inputs.copy() 

498 table["tract"] = infos[key].tract 

499 table["patch"] = infos[key].patch 

500 table["weight"] = table["coadd_weight"]*fractions[key] 

501 results.append(table) 

502 return results 

503 

504 

505def get_template_inputs(butler, dia_sources, *, instrument=None, template_dataset="template_detector"): 

506 """Find the observation times and weights of the images that make up 

507 the template at the position of each diaSource. 

508 

509 Parameters 

510 ---------- 

511 butler : `lsst.daf.butler.Butler` 

512 Butler with the templates and the coadds they were built from. 

513 Coadds are found by dataset ID, so their collections need not be 

514 searched. 

515 dia_sources : `astropy.table.Table` or `pandas.DataFrame` 

516 DiaSources, with columns ``diaSourceId``, ``visit``, ``detector``, 

517 ``ra``, ``dec`` (degrees) and ``x``, ``y`` (science image pixels). 

518 instrument : `str`, optional 

519 Instrument name, if the butler has no default. 

520 template_dataset : `str`, optional 

521 Template dataset type; must be an `lsst.images.DifferenceImage`. 

522 

523 Returns 

524 ------- 

525 inputs : `astropy.table.Table` 

526 One row per (diaSource, template coadd, input image), with columns: 

527 

528 ``diaSourceId``, ``visit``, ``detector`` 

529 The diaSource and its science image. 

530 ``tract``, ``patch`` 

531 Template coadd the input contributed through. 

532 ``input_visit``, ``input_detector`` 

533 The input image. 

534 ``mid_mjd_tai`` 

535 Mid-exposure time of the input (MJD, TAI); ``nan`` if unknown. 

536 ``exposure_time`` 

537 Exposure time of the input (s); ``nan`` if unknown. 

538 ``coadd_weight`` 

539 Weight of the input as recorded by the coadd. 

540 ``weight`` 

541 Fraction of the template flux at the diaSource that came from 

542 this input through this coadd. These sum to 1 for each 

543 diaSource; an input in several overlapping coadds has one row 

544 for each. 

545 

546 Raises 

547 ------ 

548 TypeError 

549 Raised if ``template_dataset`` has no ``templates`` component (a 

550 legacy Exposure), or if a coadd has no known provenance component. 

551 LookupError 

552 Raised if a template coadd is missing from the butler. 

553 

554 Notes 

555 ----- 

556 The weights are evaluated at the diaSource position only. Warping the 

557 template and PSF matching both mix in nearby template pixels, which can 

558 come from other coadds or inputs; that is not modeled. Per-pixel input 

559 rejection (masked pixels, or artifacts clipped by 

560 `lsst.drp.tasks.assemble_coadd.CompareWarpAssembleCoaddTask`) is not 

561 recorded in coadd provenance, so a rejected input is still reported. 

562 """ 

563 dataset_type = butler.get_dataset_type(template_dataset) 

564 if "templates" not in dataset_type.storageClass.allComponents(): 

565 raise TypeError(f"{template_dataset!r} has storage class {dataset_type.storageClass_name}, which " 

566 "has no 'templates' component; templates written as legacy Exposures are not " 

567 "supported.") 

568 if isinstance(dia_sources, pd.DataFrame): 568 ↛ 569line 568 didn't jump to line 569 because the condition on line 568 was never true

569 sources = astropy.table.Table.from_pandas(dia_sources) 

570 else: 

571 sources = astropy.table.Table(dia_sources) 

572 readers = {} 

573 results = [] 

574 pairs = np.unique(np.column_stack([sources["visit"], sources["detector"]]), axis=0) 

575 for visit, detector in pairs: 

576 selected = sources[(sources["visit"] == visit) & (sources["detector"] == detector)] 

577 data_id = {"visit": int(visit), "detector": int(detector)} 

578 if instrument is not None: 578 ↛ 579line 578 didn't jump to line 579 because the condition on line 578 was never true

579 data_id["instrument"] = instrument 

580 for table in _template_inputs_for_detector(butler, data_id, selected, template_dataset, readers): 

581 table["diaSourceId"] = selected["diaSourceId"][table["point"]] 

582 table["visit"] = visit 

583 table["detector"] = detector 

584 results.append(table) 

585 columns = ["diaSourceId", "visit", "detector", "tract", "patch"] + _INPUT_COLUMNS + ["weight"] 

586 if not results: 586 ↛ 587line 586 didn't jump to line 587 because the condition on line 586 was never true

587 empty = _empty_inputs() 

588 for name in ["diaSourceId", "visit", "detector", "tract", "patch"]: 

589 empty[name] = np.zeros(0, dtype=np.int64) 

590 empty["weight"] = np.zeros(0) 

591 return empty[columns] 

592 return astropy.table.vstack(results)[columns] 

593 

594 

595def summarize_template_inputs(inputs): 

596 """Summarize the template inputs of each diaSource. 

597 

598 Parameters 

599 ---------- 

600 inputs : `astropy.table.Table` 

601 Output of `get_template_inputs`. 

602 

603 Returns 

604 ------- 

605 summary : `astropy.table.Table` 

606 One row per diaSource, with columns ``diaSourceId``, ``n_inputs`` 

607 (distinct input images), ``mean_mjd_tai`` (weighted mean 

608 mid-exposure time), ``min_mjd_tai`` and ``max_mjd_tai``. Inputs with 

609 unknown times are left out of the time columns. 

610 """ 

611 rows = [] 

612 for group in inputs.group_by("diaSourceId").groups: 

613 known = np.isfinite(group["mid_mjd_tai"]) 

614 times = group["mid_mjd_tai"][known] 

615 weights = group["weight"][known] 

616 n_inputs = len(np.unique(np.column_stack([group["input_visit"], group["input_detector"]]), axis=0)) 

617 if len(times): 617 ↛ 621line 617 didn't jump to line 621 because the condition on line 617 was always true

618 row = (np.average(times, weights=weights) if weights.sum() > 0 else np.nan, 

619 times.min(), times.max()) 

620 else: 

621 row = (np.nan, np.nan, np.nan) 

622 rows.append((group["diaSourceId"][0], n_inputs) + row) 

623 return astropy.table.Table( 

624 rows=rows, 

625 names=["diaSourceId", "n_inputs", "mean_mjd_tai", "min_mjd_tai", "max_mjd_tai"], 

626 dtype=[np.int64, np.int64, np.float64, np.float64, np.float64], 

627 )