Coverage for python / lsst / pipe / tasks / selectImages.py: 21%

1# This file is part of pipe_tasks. 

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__all__ = ["BaseSelectImagesTask", "BaseExposureInfo", "WcsSelectImagesTask", "PsfWcsSelectImagesTask", 

23 "DatabaseSelectImagesConfig", "BestSeeingSelectVisitsTask", 

24 "BestSeeingQuantileSelectVisitsTask"] 

25 

26import numpy as np 

27import lsst.sphgeom 

28import lsst.pex.config as pexConfig 

29import lsst.pex.exceptions as pexExceptions 

30import lsst.geom as geom 

31import lsst.pipe.base as pipeBase 

32from lsst.skymap import BaseSkyMap 

33from lsst.daf.base import DateTime 

34from lsst.utils.timer import timeMethod 

35 

36 

37def _meetsVisitSummaryMinValues(visit, visitSummary, visitSummaryMinValues, logger=None): 

38 """Check if the visitSummary meets minimum values in visitSummaryMinValues. 

39 

40 Parameters 

41 ---------- 

42 visit : `int` 

43 Visit number. 

44 visitSummary : `lsst.afw.table.ExposureCatalog` 

45 Exposure catalog with per-detector summary information. 

46 visitSummaryMinValues : `dict` 

47 Dictionary with column names as keys and minimum allowed values as values. 

48 logger : `lsst.log.Logger` 

49 Logger to log debug and warning messages. 

50 

51 Returns 

52 ------- 

53 result : `bool` 

54 True if all columns in visitSummary meet the minimum values specified 

55 in visitSummaryMinValues, False otherwise. 

56 """ 

57 for columnName, minThreshold in visitSummaryMinValues.items(): 

58 # Get values for this column across all detectors with a valid WCS 

59 values = np.asarray([vs[columnName] for vs in visitSummary if vs.getWcs()]) 

60 meanValue = np.nanmean(values) 

61 if meanValue < minThreshold: 

62 if logger is not None: 

63 logger.info(f'Rejecting visit {visit}: mean {columnName} ({meanValue:.3f}) ' 

64 f'is below the minimum threshold ({minThreshold:.3f}).') 

65 return False 

66 return True 

67 

68 

69class DatabaseSelectImagesConfig(pexConfig.Config): 

70 """Base configuration for subclasses of BaseSelectImagesTask that use a 

71 database. 

72 """ 

73 

74 host = pexConfig.Field( 

75 doc="Database server host name", 

76 dtype=str, 

77 ) 

78 port = pexConfig.Field( 

79 doc="Database server port", 

80 dtype=int, 

81 ) 

82 database = pexConfig.Field( 

83 doc="Name of database", 

84 dtype=str, 

85 ) 

86 maxExposures = pexConfig.Field( 

87 doc="maximum exposures to select; intended for debugging; ignored if None", 

88 dtype=int, 

89 optional=True, 

90 ) 

91 

92 

93class BaseExposureInfo(pipeBase.Struct): 

94 """Data about a selected exposure. 

95 

96 Parameters 

97 ---------- 

98 dataId : `dict` 

99 Data ID keys of exposure. 

100 coordList : `list` [`lsst.afw.geom.SpherePoint`] 

101 ICRS coordinates of the corners of the exposure 

102 plus any others items that are desired. 

103 """ 

104 

105 def __init__(self, dataId, coordList): 

106 super(BaseExposureInfo, self).__init__(dataId=dataId, coordList=coordList) 

107 

108 

109class BaseSelectImagesTask(pipeBase.Task): 

110 """Base task for selecting images suitable for coaddition. 

111 """ 

112 

113 ConfigClass = pexConfig.Config 

114 _DefaultName = "selectImages" 

115 

116 @timeMethod 

117 def run(self, coordList): 

118 """Select images suitable for coaddition in a particular region. 

119 

120 Parameters 

121 ---------- 

122 coordList : `list` [`lsst.geom.SpherePoint`] or `None` 

123 List of coordinates defining region of interest; if `None`, then 

124 select all images subclasses may add additional keyword arguments, 

125 as required. 

126 

127 Returns 

128 ------- 

129 result : `pipeBase.Struct` 

130 Results as a struct with attributes: 

131 

132 ``exposureInfoList`` 

133 A list of exposure information objects (subclasses of 

134 BaseExposureInfo), which have at least the following fields: 

135 - dataId: Data ID dictionary (`dict`). 

136 - coordList: ICRS coordinates of the corners of the exposure. 

137 (`list` [`lsst.geom.SpherePoint`]) 

138 """ 

139 raise NotImplementedError() 

140 

141 

142def _extractKeyValue(dataList, keys=None): 

143 """Extract the keys and values from a list of dataIds. 

144 

145 The input dataList is a list of objects that have 'dataId' members. 

146 This allows it to be used for both a list of data references and a 

147 list of ExposureInfo. 

148 

149 Parameters 

150 ---------- 

151 dataList : `Unknown` 

152 keys : `Unknown` 

153 

154 Returns 

155 ------- 

156 keys : `Unknown` 

157 values : `Unknown` 

158 

159 Raises 

160 ------ 

161 RuntimeError 

162 Raised if DataId keys are inconsistent. 

163 """ 

164 assert len(dataList) > 0 

165 if keys is None: 

166 keys = sorted(dataList[0].dataId.keys()) 

167 keySet = set(keys) 

168 values = list() 

169 for data in dataList: 

170 thisKeys = set(data.dataId.keys()) 

171 if thisKeys != keySet: 

172 raise RuntimeError("DataId keys inconsistent: %s vs %s" % (keySet, thisKeys)) 

173 values.append(tuple(data.dataId[k] for k in keys)) 

174 return keys, values 

175 

176 

177class SelectStruct(pipeBase.Struct): 

178 """A container for data to be passed to the WcsSelectImagesTask. 

179 

180 Parameters 

181 ---------- 

182 dataRef : `Unknown` 

183 Data reference. 

184 wcs : `lsst.afw.geom.SkyWcs` 

185 Coordinate system definition (wcs). 

186 bbox : `lsst.geom.box.Box2I` 

187 Integer bounding box for image. 

188 """ 

189 

190 def __init__(self, dataRef, wcs, bbox): 

191 super(SelectStruct, self).__init__(dataRef=dataRef, wcs=wcs, bbox=bbox) 

192 

193 

194class WcsSelectImagesConfig(pexConfig.Config): 

195 excludeDetectors = pexConfig.ListField( 

196 dtype=int, 

197 default=[], 

198 doc="Detectors to exclude from selection.", 

199 optional=True, 

200 ) 

201 

202 

203class WcsSelectImagesTask(BaseSelectImagesTask): 

204 """Select images using their Wcs. 

205 

206 We use the "convexHull" method of lsst.sphgeom.ConvexPolygon to define 

207 polygons on the celestial sphere, and test the polygon of the 

208 patch for overlap with the polygon of the image. 

209 

210 We use "convexHull" instead of generating a ConvexPolygon 

211 directly because the standard for the inputs to ConvexPolygon 

212 are pretty high and we don't want to be responsible for reaching them. 

213 """ 

214 

215 ConfigClass = WcsSelectImagesConfig 

216 _DefaultName = "WcsSelectImages" 

217 

218 def run(self, wcsList, bboxList, coordList, dataIds=None, **kwargs): 

219 """Return indices of provided lists that meet the selection criteria. 

220 

221 Parameters 

222 ---------- 

223 wcsList : `list` [`lsst.afw.geom.SkyWcs`] 

224 Specifying the WCS's of the input ccds to be selected. 

225 bboxList : `list` [`lsst.geom.Box2I`] 

226 Specifying the bounding boxes of the input ccds to be selected. 

227 coordList : `list` [`lsst.geom.SpherePoint`] 

228 ICRS coordinates specifying boundary of the patch. 

229 dataIds : iterable [`lsst.daf.butler.dataId`] or `None`, optional 

230 An iterable object of dataIds which point to reference catalogs. 

231 **kwargs 

232 Additional keyword arguments. 

233 

234 Returns 

235 ------- 

236 result : `list` [`int`] 

237 The indices of selected ccds. 

238 """ 

239 if dataIds is None: 

240 dataIds = [None] * len(wcsList) 

241 patchVertices = [coord.getVector() for coord in coordList] 

242 patchPoly = lsst.sphgeom.ConvexPolygon.convexHull(patchVertices) 

243 result = [] 

244 for i, (imageWcs, imageBox, dataId) in enumerate(zip(wcsList, bboxList, dataIds)): 

245 if dataId and ("detector" in dataId) and (dataId["detector"] in self.config.excludeDetectors): 

246 self.log.info("De-selecting exposure %s because detector %s is exluded from processing", 

247 dataId, dataId["detector"]) 

248 elif imageWcs is None: 

249 self.log.info("De-selecting exposure %s: Exposure has no WCS.", dataId) 

250 else: 

251 imageCorners = self.getValidImageCorners(imageWcs, imageBox, patchPoly, dataId) 

252 if imageCorners: 

253 result.append(i) 

254 return result 

255 

256 def getValidImageCorners(self, imageWcs, imageBox, patchPoly, dataId=None): 

257 """Return corners or `None` if bad. 

258 

259 Parameters 

260 ---------- 

261 imageWcs : `Unknown` 

262 imageBox : `Unknown` 

263 patchPoly : `Unknown` 

264 dataId : `Unknown` 

265 """ 

266 try: 

267 imageCorners = [imageWcs.pixelToSky(pix) for pix in geom.Box2D(imageBox).getCorners()] 

268 except (pexExceptions.DomainError, pexExceptions.RuntimeError) as e: 

269 # Protecting ourselves from awful Wcs solutions in input images 

270 self.log.debug("WCS error in testing calexp %s (%s): deselecting", dataId, e) 

271 return None 

272 

273 imagePoly = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in imageCorners]) 

274 if imagePoly is None: 

275 self.log.debug("Unable to create polygon from image %s: deselecting", dataId) 

276 return None 

277 

278 if patchPoly.intersects(imagePoly): 

279 # "intersects" also covers "contains" or "is contained by" 

280 self.log.info("Selecting calexp %s", dataId) 

281 return imageCorners 

282 

283 return None 

284 

285 

286class PsfWcsSelectImagesConnections(pipeBase.PipelineTaskConnections, 

287 dimensions=("tract", "patch", "skymap", "instrument", "visit"), 

288 defaultTemplates={"coaddName": "deep"}): 

289 pass 

290 

291 

292class PsfWcsSelectImagesConfig(pipeBase.PipelineTaskConfig, 

293 pipelineConnections=PsfWcsSelectImagesConnections): 

294 maxPsfFwhm = pexConfig.Field( 

295 doc="Maximum PSF FWHM (in arcseconds) to select. This can be used in conjunction with the " 

296 "per-visit maxPsfFwhm in BestSeeingSelectVisitsTask to ensure that no outlier CCD has a " 

297 "PSF FWHM larger than this value.", 

298 dtype=float, 

299 default=1.8, 

300 optional=True, 

301 ) 

302 maxEllipResidual = pexConfig.Field( 

303 doc="Maximum median ellipticity residual", 

304 dtype=float, 

305 default=0.007, 

306 optional=True, 

307 ) 

308 maxSizeScatter = pexConfig.Field( 

309 doc="Maximum scatter in the size residuals", 

310 dtype=float, 

311 optional=True, 

312 ) 

313 maxScaledSizeScatter = pexConfig.Field( 

314 doc="Maximum scatter in the size residuals, scaled by the median size", 

315 dtype=float, 

316 default=0.019, 

317 optional=True, 

318 ) 

319 maxPsfTraceRadiusDelta = pexConfig.Field( 

320 doc="Maximum delta (max - min) of model PSF trace radius values evaluated on a grid on " 

321 "the unmasked detector pixels (pixel).", 

322 dtype=float, 

323 default=0.7, 

324 optional=True, 

325 ) 

326 maxPsfApFluxDelta = pexConfig.Field( 

327 doc="Maximum delta (max - min) of model PSF aperture flux (with aperture radius of " 

328 "max(2, 3*psfSigma)) values evaluated on a grid on the unmasked detector pixels (based " 

329 "on a normalized-to-one flux).", 

330 dtype=float, 

331 default=0.24, 

332 optional=True, 

333 ) 

334 maxPsfApCorrSigmaScaledDelta = pexConfig.Field( 

335 doc="Maximum delta (max - min) of model PSF aperture correction values evaluated on a grid " 

336 "on the unmasked detector pixels scaled (divided) by the measured model psfSigma.", 

337 dtype=float, 

338 default=0.22, 

339 optional=True, 

340 ) 

341 minNPsfStarPerBand = pexConfig.DictField( 

342 keytype=str, 

343 itemtype=int, 

344 default={ 

345 "u": 6, 

346 "g": 6, 

347 "r": 6, 

348 "i": 6, 

349 "z": 6, 

350 "y": 6, 

351 "fallback": 6, 

352 }, 

353 doc="Minimum number of PSF stars for the final PSF model to be considered " 

354 "well-constrained and suitible for inclusion in the coadd. This number should " 

355 "take into consideration the spatial order used for the PSF model. If the current " 

356 "band for the exposure is not included as a key in this dict, the value associated " 

357 "with the \"fallback\" key will be used.", 

358 ) 

359 maxStarEPerBand = pexConfig.DictField( 

360 keytype=str, 

361 itemtype=float, 

362 default={ 

363 "u": 0.4, 

364 "g": 0.4, 

365 "r": 0.4, 

366 "i": 0.4, 

367 "z": 0.4, 

368 "y": 0.4, 

369 "fallback": 0.4, 

370 }, 

371 doc="Maximum median of the ellipticity (sqrt(starE1**2.0 + starE2**2.0)) " 

372 "distribution of the star sources used in the PSF model. If the current band " 

373 "for the exposure is not included as a key in this dict, the value associated " 

374 "with the \"fallback\" key will be used.", 

375 ) 

376 maxStarUnNormalizedEPerBand = pexConfig.DictField( 

377 keytype=str, 

378 itemtype=float, 

379 default={ 

380 "u": 2.8, 

381 "g": 2.8, 

382 "r": 2.8, 

383 "i": 2.8, 

384 "z": 2.8, 

385 "y": 2.8, 

386 "fallback": 2.8, 

387 }, 

388 doc="Maximum median of the unnormalized ellipticity " 

389 "(sqrt((starXX - starYY)**2.0 + (2.0*starXY)**2.0)) distribution of the star " 

390 "sources used in the PSF model. If the current band for the exposure is not " 

391 "included as a key in this dict, the value associated with the \"fallback\" key " 

392 "will be used.", 

393 ) 

394 excludeDetectors = pexConfig.ListField( 

395 dtype=int, 

396 default=[], 

397 doc="Detectors to exclude from selection.", 

398 optional=True, 

399 ) 

400 

401 def validate(self): 

402 super().validate() 

403 if "fallback" not in self.minNPsfStarPerBand: 

404 msg = ("Must include a \"fallback\" key in the config.minNPsfStarPerBand config dict. " 

405 f"It is currenly: {self.minNPsfStarPerBand}.") 

406 raise pexConfig.FieldValidationError(self.__class__.minNPsfStarPerBand, self, msg) 

407 if "fallback" not in self.maxStarEPerBand: 

408 msg = ("Must include a \"fallback\" key in the config.maxStarEPerBand config dict. " 

409 f"It is currenly: {self.maxStarEPerBand}.") 

410 raise pexConfig.FieldValidationError(self.__class__.maxStarEPerBand, self, msg) 

411 if "fallback" not in self.maxStarUnNormalizedEPerBand: 

412 msg = ("Must include a \"fallback\" key in the config.maxStarUnNormalizedEPerBand " 

413 f"config dict. It is currenly: {self.maxStarUnNormalizedEPerBand}.") 

414 raise pexConfig.FieldValidationError(self.__class__.maxStarUnNormalizedEPerBand, self, msg) 

415 

416 

417class PsfWcsSelectImagesTask(WcsSelectImagesTask): 

418 """Select images using their Wcs and cuts on the PSF properties. 

419 

420 The PSF quality criteria are based on the size and ellipticity 

421 residuals from the adaptive second moments of the star and the PSF. 

422 

423 The criteria are: 

424 - the median of the ellipticty residuals. 

425 - the robust scatter of the size residuals (using the median absolute 

426 deviation). 

427 - the robust scatter of the size residuals scaled by the square of 

428 the median size. 

429 """ 

430 

431 ConfigClass = PsfWcsSelectImagesConfig 

432 _DefaultName = "PsfWcsSelectImages" 

433 

434 def run(self, wcsList, bboxList, coordList, visitSummary, dataIds=None, **kwargs): 

435 """Return indices of provided lists that meet the selection criteria. 

436 

437 Parameters 

438 ---------- 

439 wcsList : `list` [`lsst.afw.geom.SkyWcs`] 

440 Specifying the WCS's of the input ccds to be selected. 

441 bboxList : `list` [`lsst.geom.Box2I`] 

442 Specifying the bounding boxes of the input ccds to be selected. 

443 coordList : `list` [`lsst.geom.SpherePoint`] 

444 ICRS coordinates specifying boundary of the patch. 

445 visitSummary : `list` [`lsst.afw.table.ExposureCatalog`] 

446 containing the PSF shape information for the input ccds to be 

447 selected. 

448 dataIds : iterable [`lsst.daf.butler.dataId`] or `None`, optional 

449 An iterable object of dataIds which point to reference catalogs. 

450 **kwargs 

451 Additional keyword arguments. 

452 

453 Returns 

454 ------- 

455 goodPsf : `list` [`int`] 

456 The indices of selected ccds. 

457 """ 

458 goodWcs = super(PsfWcsSelectImagesTask, self).run(wcsList=wcsList, bboxList=bboxList, 

459 coordList=coordList, dataIds=dataIds) 

460 

461 goodPsf = [] 

462 

463 for i, dataId in enumerate(dataIds): 

464 if i not in goodWcs: 

465 continue 

466 if self.isValid(visitSummary, dataId["detector"]): 

467 goodPsf.append(i) 

468 

469 return goodPsf 

470 

471 def isValid(self, visitSummary, detectorId): 

472 """Should this ccd be selected based on its PSF shape information. 

473 

474 Parameters 

475 ---------- 

476 visitSummary : `lsst.afw.table.ExposureCatalog` 

477 Exposure catalog with per-detector summary information. 

478 detectorId : `int` 

479 Detector identifier. 

480 

481 Returns 

482 ------- 

483 valid : `bool` 

484 True if selected. 

485 """ 

486 row = visitSummary.find(detectorId) 

487 if row is None: 

488 # This is not listed, so it must be bad. 

489 self.log.warning("Removing detector %d because summary stats not available.", detectorId) 

490 return False 

491 

492 band = row["band"] 

493 if band in self.config.minNPsfStarPerBand: 

494 minNPsfStar = self.config.minNPsfStarPerBand[band] 

495 else: 

496 minNPsfStar = self.config.minNPsfStarPerBand["fallback"] 

497 if band in self.config.maxStarEPerBand: 

498 maxStarE = self.config.maxStarEPerBand[band] 

499 else: 

500 maxStarE = self.config.maxStarEPerBand["fallback"] 

501 if band in self.config.maxStarUnNormalizedEPerBand: 

502 maxStarUnNormalizedE = self.config.maxStarUnNormalizedEPerBand[band] 

503 else: 

504 maxStarUnNormalizedE = self.config.maxStarUnNormalizedEPerBand["fallback"] 

505 nPsfStar = row["nPsfStar"] 

506 medianE = np.sqrt(row["psfStarDeltaE1Median"]**2. + row["psfStarDeltaE2Median"]**2.) 

507 scatterSize = row["psfStarDeltaSizeScatter"] 

508 scaledScatterSize = row["psfStarScaledDeltaSizeScatter"] 

509 psfTraceRadiusDelta = row["psfTraceRadiusDelta"] 

510 psfApFluxDelta = row["psfApFluxDelta"] 

511 psfApCorrSigmaScaledDelta = row["psfApCorrSigmaScaledDelta"] 

512 starEMedian = row["starEMedian"] 

513 starUnNormalizedEMedian = row["starUnNormalizedEMedian"] 

514 

515 valid = True 

516 if self.config.maxEllipResidual and not (medianE <= self.config.maxEllipResidual): 

517 self.log.info("Removing visit %d detector %d because median e residual too large: %f vs %f", 

518 row["visit"], detectorId, medianE, self.config.maxEllipResidual) 

519 valid = False 

520 elif self.config.maxSizeScatter and not (scatterSize <= self.config.maxSizeScatter): 

521 self.log.info("Removing visit %d detector %d because size scatter too large: %f vs %f", 

522 row["visit"], detectorId, scatterSize, self.config.maxSizeScatter) 

523 valid = False 

524 elif self.config.maxScaledSizeScatter and not (scaledScatterSize <= self.config.maxScaledSizeScatter): 

525 self.log.info("Removing visit %d detector %d because scaled size scatter too large: %f vs %f", 

526 row["visit"], detectorId, scaledScatterSize, self.config.maxScaledSizeScatter) 

527 valid = False 

528 elif ( 

529 self.config.maxPsfTraceRadiusDelta is not None 

530 and not (psfTraceRadiusDelta <= self.config.maxPsfTraceRadiusDelta) 

531 ): 

532 self.log.info( 

533 "Removing visit %d detector %d because max-min delta of model PSF trace radius values " 

534 "across the unmasked detector pixels is not finite or too large: %.3f vs %.3f (pixels)", 

535 row["visit"], detectorId, psfTraceRadiusDelta, self.config.maxPsfTraceRadiusDelta 

536 ) 

537 valid = False 

538 elif ( 

539 self.config.maxPsfApFluxDelta is not None 

540 and not (psfApFluxDelta <= self.config.maxPsfApFluxDelta) 

541 ): 

542 self.log.info( 

543 "Removing visit %d detector %d because max-min delta of model PSF aperture flux values " 

544 "across the unmasked detector pixels is not finite or too large: %.3f vs %.3f (pixels)", 

545 row["visit"], detectorId, psfApFluxDelta, self.config.maxPsfApFluxDelta 

546 ) 

547 valid = False 

548 elif ( 

549 self.config.maxPsfApCorrSigmaScaledDelta is not None 

550 and not (psfApCorrSigmaScaledDelta <= self.config.maxPsfApCorrSigmaScaledDelta) 

551 ): 

552 self.log.info( 

553 "Removing visit %d detector %d because max-min delta of the model PSF apterture correction " 

554 "values scaled (divided) by the psfSigma across the unmasked detector pixels is not " 

555 "finite or too large: %.3f vs %.3f (fatcor)", 

556 row["visit"], detectorId, psfApCorrSigmaScaledDelta, self.config.maxPsfApCorrSigmaScaledDelta 

557 ) 

558 valid = False 

559 elif minNPsfStar and (nPsfStar < minNPsfStar): 

560 self.log.info( 

561 "Removing visit %d detector %d because the PSF model had too few stars: %d vs %d", 

562 row["visit"], detectorId, nPsfStar, minNPsfStar 

563 ) 

564 valid = False 

565 elif self.config.maxPsfFwhm: 

566 pixelScale = row['pixelScale'] 

567 psfSigma = row['psfSigma'] 

568 fwhm = psfSigma * pixelScale * np.sqrt(8.*np.log(2.)) 

569 if not (fwhm <= self.config.maxPsfFwhm): 

570 self.log.info("Removing visit %d detector %d because FWHM too large: %f vs %f", 

571 row["visit"], detectorId, fwhm, self.config.maxPsfFwhm) 

572 valid = False 

573 elif maxStarE and not (starEMedian <= maxStarE): 

574 self.log.info( 

575 "Removing visit %d detector %d because the median star ellipticity is too large: %d vs %d", 

576 row["visit"], detectorId, starEMedian, maxStarE 

577 ) 

578 valid = False 

579 elif maxStarUnNormalizedE and not (starUnNormalizedEMedian <= maxStarUnNormalizedE): 

580 self.log.info( 

581 "Removing visit %d detector %d because the median star unnormalized ellipticity " 

582 "is too large: %d vs %d", 

583 row["visit"], detectorId, starUnNormalizedEMedian, maxStarUnNormalizedE 

584 ) 

585 valid = False 

586 

587 return valid 

588 

589 

590class BestSeeingSelectVisitsConnections(pipeBase.PipelineTaskConnections, 

591 dimensions=("tract", "patch", "skymap", "band", "instrument"), 

592 defaultTemplates={"coaddName": "goodSeeing", 

593 "calexpType": ""}): 

594 skyMap = pipeBase.connectionTypes.Input( 

595 doc="Input definition of geometry/bbox and projection/wcs for coadded exposures", 

596 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME, 

597 storageClass="SkyMap", 

598 dimensions=("skymap",), 

599 ) 

600 visitSummaries = pipeBase.connectionTypes.Input( 

601 doc="Per-visit consolidated exposure metadata", 

602 name="finalVisitSummary", 

603 storageClass="ExposureCatalog", 

604 dimensions=("instrument", "visit",), 

605 multiple=True, 

606 deferLoad=True 

607 ) 

608 goodVisits = pipeBase.connectionTypes.Output( 

609 doc="Selected visits to be coadded.", 

610 name="{coaddName}Visits", 

611 storageClass="StructuredDataDict", 

612 dimensions=("instrument", "tract", "patch", "skymap", "band"), 

613 ) 

614 

615 

616class BestSeeingSelectVisitsConfig(pipeBase.PipelineTaskConfig, 

617 pipelineConnections=BestSeeingSelectVisitsConnections): 

618 nVisitsMax = pexConfig.RangeField( 

619 dtype=int, 

620 doc="Maximum number of visits to select; use -1 to select all.", 

621 default=12, 

622 min=-1, 

623 ) 

624 maxPsfFwhm = pexConfig.Field( 

625 dtype=float, 

626 doc="Maximum PSF FWHM (in arcseconds) to select", 

627 default=1.5, 

628 optional=True 

629 ) 

630 minPsfFwhm = pexConfig.Field( 

631 dtype=float, 

632 doc="Minimum PSF FWHM (in arcseconds) to select", 

633 default=0., 

634 optional=True 

635 ) 

636 doConfirmOverlap = pexConfig.Field( 

637 dtype=bool, 

638 doc="Do remove visits that do not actually overlap the patch?", 

639 default=True, 

640 ) 

641 minMJD = pexConfig.Field( 

642 dtype=float, 

643 doc="Minimum visit MJD to select", 

644 default=None, 

645 optional=True 

646 ) 

647 maxMJD = pexConfig.Field( 

648 dtype=float, 

649 doc="Maximum visit MJD to select", 

650 default=None, 

651 optional=True 

652 ) 

653 visitSummaryMinValues = pexConfig.DictField( 

654 keytype=str, 

655 itemtype=float, 

656 doc=("Optional thresholds for visit selection. Keys are visit_summary column names" 

657 "(e.g. 'effTimeZeroPointScale'), and values are minimum allowed values."), 

658 default=None, 

659 optional=True 

660 ) 

661 

662 

663class BestSeeingSelectVisitsTask(pipeBase.PipelineTask): 

664 """Select up to a maximum number of the best-seeing visits. 

665 

666 Don't exceed the FWHM range specified by configs min(max)PsfFwhm. 

667 This Task is a port of the Gen2 image-selector used in the AP pipeline: 

668 BestSeeingSelectImagesTask. This Task selects full visits based on the 

669 average PSF of the entire visit. 

670 """ 

671 

672 ConfigClass = BestSeeingSelectVisitsConfig 

673 _DefaultName = 'bestSeeingSelectVisits' 

674 

675 def runQuantum(self, butlerQC, inputRefs, outputRefs): 

676 inputs = butlerQC.get(inputRefs) 

677 quantumDataId = butlerQC.quantum.dataId 

678 outputs = self.run(**inputs, dataId=quantumDataId) 

679 butlerQC.put(outputs, outputRefs) 

680 

681 def run(self, visitSummaries, skyMap, dataId): 

682 """Run task. 

683 

684 Parameters 

685 ---------- 

686 visitSummary : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`] 

687 List of `lsst.pipe.base.connections.DeferredDatasetRef` of 

688 visitSummary tables of type `lsst.afw.table.ExposureCatalog`. 

689 skyMap : `lsst.skyMap.SkyMap` 

690 SkyMap for checking visits overlap patch. 

691 dataId : `dict` of dataId keys 

692 For retrieving patch info for checking visits overlap patch. 

693 

694 Returns 

695 ------- 

696 result : `lsst.pipe.base.Struct` 

697 Results as a struct with attributes: 

698 

699 ``goodVisits`` 

700 A `dict` with selected visit ids as keys, 

701 so that it can be be saved as a StructuredDataDict. 

702 StructuredDataList's are currently limited. 

703 """ 

704 if self.config.doConfirmOverlap: 

705 patchPolygon = self.makePatchPolygon(skyMap, dataId) 

706 

707 inputVisits = [visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries] 

708 fwhmSizes = [] 

709 visits = [] 

710 for visit, visitSummary in zip(inputVisits, visitSummaries): 

711 # read in one-by-one and only once. There may be hundreds 

712 visitSummary = visitSummary.get() 

713 

714 # mjd is guaranteed to be the same for every detector in the 

715 # visitSummary. 

716 mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD) 

717 

718 pixelScales = np.array([vs['pixelScale'] for vs in visitSummary if vs.getWcs()]) 

719 # psfSigma is PSF model determinant radius at chip center in pixels 

720 psfSigmas = np.array([vs['psfSigma'] for vs in visitSummary if vs.getWcs()]) 

721 fwhm = np.nanmean(psfSigmas * pixelScales) * np.sqrt(8.*np.log(2.)) 

722 

723 if self.config.maxPsfFwhm and fwhm > self.config.maxPsfFwhm: 

724 continue 

725 if self.config.minPsfFwhm and fwhm < self.config.minPsfFwhm: 

726 continue 

727 if self.config.minMJD and mjd < self.config.minMJD: 

728 self.log.debug('MJD %f earlier than %.2f; rejecting', mjd, self.config.minMJD) 

729 continue 

730 if self.config.maxMJD and mjd > self.config.maxMJD: 

731 self.log.debug('MJD %f later than %.2f; rejecting', mjd, self.config.maxMJD) 

732 continue 

733 if self.config.doConfirmOverlap and not self.doesIntersectPolygon(visitSummary, patchPolygon): 

734 continue 

735 if (self.config.visitSummaryMinValues 

736 and not _meetsVisitSummaryMinValues(visit, visitSummary, self.config.visitSummaryMinValues, 

737 self.log)): 

738 continue 

739 

740 fwhmSizes.append(fwhm) 

741 visits.append(visit) 

742 

743 sortedVisits = [ind for (_, ind) in sorted(zip(fwhmSizes, visits))] 

744 if self.config.nVisitsMax < 0: 

745 output = sortedVisits 

746 else: 

747 output = sortedVisits[:self.config.nVisitsMax] 

748 

749 if len(output) == 0: 

750 self.log.info("All images rejected in BestSeeingSelectVisitsTask.") 

751 raise pipeBase.NoWorkFound(f"No good images found for {dataId}") 

752 else: 

753 self.log.info( 

754 "%d images selected with FWHM range of %f--%f arcseconds", 

755 len(output), 

756 fwhmSizes[visits.index(output[0])], 

757 fwhmSizes[visits.index(output[-1])], 

758 ) 

759 

760 # In order to store as a StructuredDataDict, convert list to dict 

761 goodVisits = {key: True for key in output} 

762 return pipeBase.Struct(goodVisits=goodVisits) 

763 

764 def makePatchPolygon(self, skyMap, dataId): 

765 """Return True if sky polygon overlaps visit. 

766 

767 Parameters 

768 ---------- 

769 skyMap : `lsst.afw.table.ExposureCatalog` 

770 Exposure catalog with per-detector geometry. 

771 dataId : `dict` of dataId keys 

772 For retrieving patch info. 

773 

774 Returns 

775 ------- 

776 result : `lsst.sphgeom.ConvexPolygon.convexHull` 

777 Polygon of patch's outer bbox. 

778 """ 

779 wcs = skyMap[dataId['tract']].getWcs() 

780 bbox = skyMap[dataId['tract']][dataId['patch']].getOuterBBox() 

781 sphCorners = wcs.pixelToSky(lsst.geom.Box2D(bbox).getCorners()) 

782 result = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in sphCorners]) 

783 return result 

784 

785 def doesIntersectPolygon(self, visitSummary, polygon): 

786 """Return True if sky polygon overlaps visit. 

787 

788 Parameters 

789 ---------- 

790 visitSummary : `lsst.afw.table.ExposureCatalog` 

791 Exposure catalog with per-detector geometry. 

792 polygon :` lsst.sphgeom.ConvexPolygon.convexHull` 

793 Polygon to check overlap. 

794 

795 Returns 

796 ------- 

797 doesIntersect : `bool` 

798 True if the visit overlaps the polygon. 

799 """ 

800 doesIntersect = False 

801 for detectorSummary in visitSummary: 

802 if (np.all(np.isfinite(detectorSummary['raCorners'])) 

803 and np.all(np.isfinite(detectorSummary['decCorners']))): 

804 corners = [lsst.geom.SpherePoint(ra, dec, units=lsst.geom.degrees).getVector() 

805 for (ra, dec) in 

806 zip(detectorSummary['raCorners'], detectorSummary['decCorners'])] 

807 detectorPolygon = lsst.sphgeom.ConvexPolygon.convexHull(corners) 

808 if detectorPolygon.intersects(polygon): 

809 doesIntersect = True 

810 break 

811 return doesIntersect 

812 

813 

814class BestSeeingQuantileSelectVisitsConfig(pipeBase.PipelineTaskConfig, 

815 pipelineConnections=BestSeeingSelectVisitsConnections): 

816 qMin = pexConfig.RangeField( 

817 doc="Lower bound of quantile range to select. Sorts visits by seeing from narrow to wide, " 

818 "and select those in the interquantile range (qMin, qMax). Set qMin to 0 for Best Seeing. " 

819 "This config should be changed from zero only for exploratory diffIm testing.", 

820 dtype=float, 

821 default=0, 

822 min=0, 

823 max=1, 

824 ) 

825 qMax = pexConfig.RangeField( 

826 doc="Upper bound of quantile range to select. Sorts visits by seeing from narrow to wide, " 

827 "and select those in the interquantile range (qMin, qMax). Set qMax to 1 for Worst Seeing.", 

828 dtype=float, 

829 default=0.33, 

830 min=0, 

831 max=1, 

832 ) 

833 maxPsfFwhm = pexConfig.Field( 

834 dtype=float, 

835 doc="Maximum PSF FWHM (in arcseconds) to select. Any visits with larger FWHM will be excluded " 

836 "before quantile selection.", 

837 default=None, 

838 optional=True 

839 ) 

840 nVisitsMin = pexConfig.Field( 

841 doc="The minimum number of visits to select, if qMin and qMax alone would have " 

842 "selected fewer. In regimes with many visits, at least this number of visits will be " 

843 "selected, superceding quantile when necessary. " 

844 "For example, if 10 visits cover this patch, qMin=0, qMax=0.33, and nVisitsMin=5, " 

845 "the best 5 visits will be selected, even though 5 > 0.33*10. " 

846 "In regimes with few visits, all available visits will be selected. " 

847 "For example, if 2 visits cover this patch and nVisitsMin=12, " 

848 "both visits will be selected regardless of qMin and qMax.", 

849 dtype=int, 

850 default=6, 

851 ) 

852 doConfirmOverlap = pexConfig.Field( 

853 dtype=bool, 

854 doc="Do remove visits that do not actually overlap the patch?", 

855 default=True, 

856 ) 

857 minMJD = pexConfig.Field( 

858 dtype=float, 

859 doc="Minimum visit MJD to select", 

860 default=None, 

861 optional=True 

862 ) 

863 maxMJD = pexConfig.Field( 

864 dtype=float, 

865 doc="Maximum visit MJD to select", 

866 default=None, 

867 optional=True 

868 ) 

869 visitSummaryMinValues = pexConfig.DictField( 

870 keytype=str, 

871 itemtype=float, 

872 doc=("Optional thresholds for visit selection. Keys are visit_summary column names" 

873 "(e.g. 'effTimeZeroPointScale'), and values are minimum allowed values."), 

874 default=None, 

875 optional=True 

876 ) 

877 

878 

879class BestSeeingQuantileSelectVisitsTask(BestSeeingSelectVisitsTask): 

880 """Select a quantile of the best-seeing visits. 

881 

882 Selects the best (for example, third) full visits based on the average 

883 PSF width in the entire visit. It can also be used for difference imaging 

884 experiments that require templates with the worst seeing visits. 

885 For example, selecting the worst third can be acheived by 

886 changing the config parameters qMin to 0.66 and qMax to 1.