Coverage for python/lsst/pipe/tasks/selectImages.py: 23%
258 statements
« prev ^ index » next coverage.py v7.10.6, created at 2025-10-11 02:06 -0700
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/>.
22__all__ = ["BaseSelectImagesTask", "BaseExposureInfo", "WcsSelectImagesTask", "PsfWcsSelectImagesTask",
23 "DatabaseSelectImagesConfig", "BestSeeingSelectVisitsTask",
24 "BestSeeingQuantileSelectVisitsTask"]
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
37def _meetsVisitSummaryMinValues(visit, visitSummary, visitSummaryMinValues, logger=None):
38 """Check if the visitSummary meets minimum values in visitSummaryMinValues.
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.
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
69class DatabaseSelectImagesConfig(pexConfig.Config):
70 """Base configuration for subclasses of BaseSelectImagesTask that use a
71 database.
72 """
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 )
93class BaseExposureInfo(pipeBase.Struct):
94 """Data about a selected exposure.
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 """
105 def __init__(self, dataId, coordList):
106 super(BaseExposureInfo, self).__init__(dataId=dataId, coordList=coordList)
109class BaseSelectImagesTask(pipeBase.Task):
110 """Base task for selecting images suitable for coaddition.
111 """
113 ConfigClass = pexConfig.Config
114 _DefaultName = "selectImages"
116 @timeMethod
117 def run(self, coordList):
118 """Select images suitable for coaddition in a particular region.
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.
127 Returns
128 -------
129 result : `pipeBase.Struct`
130 Results as a struct with attributes:
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()
142def _extractKeyValue(dataList, keys=None):
143 """Extract the keys and values from a list of dataIds.
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.
149 Parameters
150 ----------
151 dataList : `Unknown`
152 keys : `Unknown`
154 Returns
155 -------
156 keys : `Unknown`
157 values : `Unknown`
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
177class SelectStruct(pipeBase.Struct):
178 """A container for data to be passed to the WcsSelectImagesTask.
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 """
190 def __init__(self, dataRef, wcs, bbox):
191 super(SelectStruct, self).__init__(dataRef=dataRef, wcs=wcs, bbox=bbox)
194class WcsSelectImagesConfig(pexConfig.Config):
195 excludeDetectors = pexConfig.ListField(
196 dtype=int,
197 default=[],
198 doc="Detectors to exclude from selection.",
199 optional=True,
200 )
203class WcsSelectImagesTask(BaseSelectImagesTask):
204 """Select images using their Wcs.
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.
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 """
215 ConfigClass = WcsSelectImagesConfig
216 _DefaultName = "WcsSelectImages"
218 def run(self, wcsList, bboxList, coordList, dataIds=None, **kwargs):
219 """Return indices of provided lists that meet the selection criteria.
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.
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
256 def getValidImageCorners(self, imageWcs, imageBox, patchPoly, dataId=None):
257 """Return corners or `None` if bad.
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
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
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
283 return None
286class PsfWcsSelectImagesConnections(pipeBase.PipelineTaskConnections,
287 dimensions=("tract", "patch", "skymap", "instrument", "visit"),
288 defaultTemplates={"coaddName": "deep"}):
289 pass
292class PsfWcsSelectImagesConfig(pipeBase.PipelineTaskConfig,
293 pipelineConnections=PsfWcsSelectImagesConnections):
294 maxEllipResidual = pexConfig.Field(
295 doc="Maximum median ellipticity residual",
296 dtype=float,
297 default=0.007,
298 optional=True,
299 )
300 maxSizeScatter = pexConfig.Field(
301 doc="Maximum scatter in the size residuals",
302 dtype=float,
303 optional=True,
304 )
305 maxScaledSizeScatter = pexConfig.Field(
306 doc="Maximum scatter in the size residuals, scaled by the median size",
307 dtype=float,
308 default=0.019,
309 optional=True,
310 )
311 maxPsfTraceRadiusDelta = pexConfig.Field(
312 doc="Maximum delta (max - min) of model PSF trace radius values evaluated on a grid on "
313 "the unmasked detector pixels (pixel).",
314 dtype=float,
315 default=0.7,
316 optional=True,
317 )
318 maxPsfApFluxDelta = pexConfig.Field(
319 doc="Maximum delta (max - min) of model PSF aperture flux (with aperture radius of "
320 "max(2, 3*psfSigma)) values evaluated on a grid on the unmasked detector pixels (based "
321 "on a normalized-to-one flux).",
322 dtype=float,
323 default=0.24,
324 optional=True,
325 )
326 maxPsfApCorrSigmaScaledDelta = pexConfig.Field(
327 doc="Maximum delta (max - min) of model PSF aperture correction values evaluated on a grid "
328 "on the unmasked detector pixels scaled (divided) by the measured model psfSigma.",
329 dtype=float,
330 default=0.22,
331 optional=True,
332 )
333 minNPsfStarPerBand = pexConfig.DictField(
334 keytype=str,
335 itemtype=int,
336 default={
337 "u": 6,
338 "g": 6,
339 "r": 6,
340 "i": 6,
341 "z": 6,
342 "y": 6,
343 "fallback": 6,
344 },
345 doc="Minimum number of PSF stars for the final PSF model to be considered "
346 "well-constrained and suitible for inclusion in the coadd. This number should "
347 "take into consideration the spatial order used for the PSF model. If the current "
348 "band for the exposure is not included as a key in this dict, the value associated "
349 "with the \"fallback\" key will be used.",
350 )
351 excludeDetectors = pexConfig.ListField(
352 dtype=int,
353 default=[],
354 doc="Detectors to exclude from selection.",
355 optional=True,
356 )
358 def validate(self):
359 super().validate()
360 if "fallback" not in self.minNPsfStarPerBand:
361 msg = ("Must include a \"fallback\" key in the config.minNPsfStarPerBand config dict. "
362 f"It is currenly: {self.minNPsfStarPerBand}.")
363 raise pexConfig.FieldValidationError(self.__class__.minNPsfStarPerBand, self, msg)
366class PsfWcsSelectImagesTask(WcsSelectImagesTask):
367 """Select images using their Wcs and cuts on the PSF properties.
369 The PSF quality criteria are based on the size and ellipticity
370 residuals from the adaptive second moments of the star and the PSF.
372 The criteria are:
373 - the median of the ellipticty residuals.
374 - the robust scatter of the size residuals (using the median absolute
375 deviation).
376 - the robust scatter of the size residuals scaled by the square of
377 the median size.
378 """
380 ConfigClass = PsfWcsSelectImagesConfig
381 _DefaultName = "PsfWcsSelectImages"
383 def run(self, wcsList, bboxList, coordList, visitSummary, dataIds=None, **kwargs):
384 """Return indices of provided lists that meet the selection criteria.
386 Parameters
387 ----------
388 wcsList : `list` [`lsst.afw.geom.SkyWcs`]
389 Specifying the WCS's of the input ccds to be selected.
390 bboxList : `list` [`lsst.geom.Box2I`]
391 Specifying the bounding boxes of the input ccds to be selected.
392 coordList : `list` [`lsst.geom.SpherePoint`]
393 ICRS coordinates specifying boundary of the patch.
394 visitSummary : `list` [`lsst.afw.table.ExposureCatalog`]
395 containing the PSF shape information for the input ccds to be
396 selected.
397 dataIds : iterable [`lsst.daf.butler.dataId`] or `None`, optional
398 An iterable object of dataIds which point to reference catalogs.
399 **kwargs
400 Additional keyword arguments.
402 Returns
403 -------
404 goodPsf : `list` [`int`]
405 The indices of selected ccds.
406 """
407 goodWcs = super(PsfWcsSelectImagesTask, self).run(wcsList=wcsList, bboxList=bboxList,
408 coordList=coordList, dataIds=dataIds)
410 goodPsf = []
412 for i, dataId in enumerate(dataIds):
413 if i not in goodWcs:
414 continue
415 if self.isValid(visitSummary, dataId["detector"]):
416 goodPsf.append(i)
418 return goodPsf
420 def isValid(self, visitSummary, detectorId):
421 """Should this ccd be selected based on its PSF shape information.
423 Parameters
424 ----------
425 visitSummary : `lsst.afw.table.ExposureCatalog`
426 Exposure catalog with per-detector summary information.
427 detectorId : `int`
428 Detector identifier.
430 Returns
431 -------
432 valid : `bool`
433 True if selected.
434 """
435 row = visitSummary.find(detectorId)
436 if row is None:
437 # This is not listed, so it must be bad.
438 self.log.warning("Removing detector %d because summary stats not available.", detectorId)
439 return False
441 band = row["band"]
442 if band in self.config.minNPsfStarPerBand:
443 minNPsfStar = self.config.minNPsfStarPerBand[band]
444 else:
445 minNPsfStar = self.config.minNPsfStarPerBand["fallback"]
446 nPsfStar = row["nPsfStar"]
447 medianE = np.sqrt(row["psfStarDeltaE1Median"]**2. + row["psfStarDeltaE2Median"]**2.)
448 scatterSize = row["psfStarDeltaSizeScatter"]
449 scaledScatterSize = row["psfStarScaledDeltaSizeScatter"]
450 psfTraceRadiusDelta = row["psfTraceRadiusDelta"]
451 psfApFluxDelta = row["psfApFluxDelta"]
452 psfApCorrSigmaScaledDelta = row["psfApCorrSigmaScaledDelta"]
454 valid = True
455 if self.config.maxEllipResidual and not (medianE <= self.config.maxEllipResidual):
456 self.log.info("Removing visit %d detector %d because median e residual too large: %f vs %f",
457 row["visit"], detectorId, medianE, self.config.maxEllipResidual)
458 valid = False
459 elif self.config.maxSizeScatter and not (scatterSize <= self.config.maxSizeScatter):
460 self.log.info("Removing visit %d detector %d because size scatter too large: %f vs %f",
461 row["visit"], detectorId, scatterSize, self.config.maxSizeScatter)
462 valid = False
463 elif self.config.maxScaledSizeScatter and not (scaledScatterSize <= self.config.maxScaledSizeScatter):
464 self.log.info("Removing visit %d detector %d because scaled size scatter too large: %f vs %f",
465 row["visit"], detectorId, scaledScatterSize, self.config.maxScaledSizeScatter)
466 valid = False
467 elif (
468 self.config.maxPsfTraceRadiusDelta is not None
469 and not (psfTraceRadiusDelta <= self.config.maxPsfTraceRadiusDelta)
470 ):
471 self.log.info(
472 "Removing visit %d detector %d because max-min delta of model PSF trace radius values "
473 "across the unmasked detector pixels is not finite or too large: %.3f vs %.3f (pixels)",
474 row["visit"], detectorId, psfTraceRadiusDelta, self.config.maxPsfTraceRadiusDelta
475 )
476 valid = False
477 elif (
478 self.config.maxPsfApFluxDelta is not None
479 and not (psfApFluxDelta <= self.config.maxPsfApFluxDelta)
480 ):
481 self.log.info(
482 "Removing visit %d detector %d because max-min delta of model PSF aperture flux values "
483 "across the unmasked detector pixels is not finite or too large: %.3f vs %.3f (pixels)",
484 row["visit"], detectorId, psfApFluxDelta, self.config.maxPsfApFluxDelta
485 )
486 valid = False
487 elif (
488 self.config.maxPsfApCorrSigmaScaledDelta is not None
489 and not (psfApCorrSigmaScaledDelta <= self.config.maxPsfApCorrSigmaScaledDelta)
490 ):
491 self.log.info(
492 "Removing visit %d detector %d because max-min delta of the model PSF apterture correction "
493 "values scaled (divided) by the psfSigma across the unmasked detector pixels is not "
494 "finite or too large: %.3f vs %.3f (fatcor)",
495 row["visit"], detectorId, psfApCorrSigmaScaledDelta, self.config.maxPsfApCorrSigmaScaledDelta
496 )
497 valid = False
498 elif minNPsfStar and (nPsfStar < minNPsfStar):
499 self.log.info(
500 "Removing visit %d detector %d because the PSF model had too few stars: %d vs %d",
501 row["visit"], detectorId, nPsfStar, minNPsfStar
502 )
503 valid = False
505 return valid
508class BestSeeingSelectVisitsConnections(pipeBase.PipelineTaskConnections,
509 dimensions=("tract", "patch", "skymap", "band", "instrument"),
510 defaultTemplates={"coaddName": "goodSeeing",
511 "calexpType": ""}):
512 skyMap = pipeBase.connectionTypes.Input(
513 doc="Input definition of geometry/bbox and projection/wcs for coadded exposures",
514 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
515 storageClass="SkyMap",
516 dimensions=("skymap",),
517 )
518 visitSummaries = pipeBase.connectionTypes.Input(
519 doc="Per-visit consolidated exposure metadata",
520 name="finalVisitSummary",
521 storageClass="ExposureCatalog",
522 dimensions=("instrument", "visit",),
523 multiple=True,
524 deferLoad=True
525 )
526 goodVisits = pipeBase.connectionTypes.Output(
527 doc="Selected visits to be coadded.",
528 name="{coaddName}Visits",
529 storageClass="StructuredDataDict",
530 dimensions=("instrument", "tract", "patch", "skymap", "band"),
531 )
534class BestSeeingSelectVisitsConfig(pipeBase.PipelineTaskConfig,
535 pipelineConnections=BestSeeingSelectVisitsConnections):
536 nVisitsMax = pexConfig.RangeField(
537 dtype=int,
538 doc="Maximum number of visits to select; use -1 to select all.",
539 default=12,
540 min=-1,
541 )
542 maxPsfFwhm = pexConfig.Field(
543 dtype=float,
544 doc="Maximum PSF FWHM (in arcseconds) to select",
545 default=1.5,
546 optional=True
547 )
548 minPsfFwhm = pexConfig.Field(
549 dtype=float,
550 doc="Minimum PSF FWHM (in arcseconds) to select",
551 default=0.,
552 optional=True
553 )
554 doConfirmOverlap = pexConfig.Field(
555 dtype=bool,
556 doc="Do remove visits that do not actually overlap the patch?",
557 default=True,
558 )
559 minMJD = pexConfig.Field(
560 dtype=float,
561 doc="Minimum visit MJD to select",
562 default=None,
563 optional=True
564 )
565 maxMJD = pexConfig.Field(
566 dtype=float,
567 doc="Maximum visit MJD to select",
568 default=None,
569 optional=True
570 )
571 visitSummaryMinValues = pexConfig.DictField(
572 keytype=str,
573 itemtype=float,
574 doc=("Optional thresholds for visit selection. Keys are visit_summary column names"
575 "(e.g. 'effTimeZeroPointScale'), and values are minimum allowed values."),
576 default=None,
577 optional=True
578 )
581class BestSeeingSelectVisitsTask(pipeBase.PipelineTask):
582 """Select up to a maximum number of the best-seeing visits.
584 Don't exceed the FWHM range specified by configs min(max)PsfFwhm.
585 This Task is a port of the Gen2 image-selector used in the AP pipeline:
586 BestSeeingSelectImagesTask. This Task selects full visits based on the
587 average PSF of the entire visit.
588 """
590 ConfigClass = BestSeeingSelectVisitsConfig
591 _DefaultName = 'bestSeeingSelectVisits'
593 def runQuantum(self, butlerQC, inputRefs, outputRefs):
594 inputs = butlerQC.get(inputRefs)
595 quantumDataId = butlerQC.quantum.dataId
596 outputs = self.run(**inputs, dataId=quantumDataId)
597 butlerQC.put(outputs, outputRefs)
599 def run(self, visitSummaries, skyMap, dataId):
600 """Run task.
602 Parameters
603 ----------
604 visitSummary : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
605 List of `lsst.pipe.base.connections.DeferredDatasetRef` of
606 visitSummary tables of type `lsst.afw.table.ExposureCatalog`.
607 skyMap : `lsst.skyMap.SkyMap`
608 SkyMap for checking visits overlap patch.
609 dataId : `dict` of dataId keys
610 For retrieving patch info for checking visits overlap patch.
612 Returns
613 -------
614 result : `lsst.pipe.base.Struct`
615 Results as a struct with attributes:
617 ``goodVisits``
618 A `dict` with selected visit ids as keys,
619 so that it can be be saved as a StructuredDataDict.
620 StructuredDataList's are currently limited.
621 """
622 if self.config.doConfirmOverlap:
623 patchPolygon = self.makePatchPolygon(skyMap, dataId)
625 inputVisits = [visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries]
626 fwhmSizes = []
627 visits = []
628 for visit, visitSummary in zip(inputVisits, visitSummaries):
629 # read in one-by-one and only once. There may be hundreds
630 visitSummary = visitSummary.get()
632 # mjd is guaranteed to be the same for every detector in the
633 # visitSummary.
634 mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD)
636 pixelScales = np.array([vs['pixelScale'] for vs in visitSummary if vs.getWcs()])
637 # psfSigma is PSF model determinant radius at chip center in pixels
638 psfSigmas = np.array([vs['psfSigma'] for vs in visitSummary if vs.getWcs()])
639 fwhm = np.nanmean(psfSigmas * pixelScales) * np.sqrt(8.*np.log(2.))
641 if self.config.maxPsfFwhm and fwhm > self.config.maxPsfFwhm:
642 continue
643 if self.config.minPsfFwhm and fwhm < self.config.minPsfFwhm:
644 continue
645 if self.config.minMJD and mjd < self.config.minMJD:
646 self.log.debug('MJD %f earlier than %.2f; rejecting', mjd, self.config.minMJD)
647 continue
648 if self.config.maxMJD and mjd > self.config.maxMJD:
649 self.log.debug('MJD %f later than %.2f; rejecting', mjd, self.config.maxMJD)
650 continue
651 if self.config.doConfirmOverlap and not self.doesIntersectPolygon(visitSummary, patchPolygon):
652 continue
653 if (self.config.visitSummaryMinValues
654 and not _meetsVisitSummaryMinValues(visit, visitSummary, self.config.visitSummaryMinValues,
655 self.log)):
656 continue
658 fwhmSizes.append(fwhm)
659 visits.append(visit)
661 sortedVisits = [ind for (_, ind) in sorted(zip(fwhmSizes, visits))]
662 if self.config.nVisitsMax < 0:
663 output = sortedVisits
664 else:
665 output = sortedVisits[:self.config.nVisitsMax]
667 if len(output) == 0:
668 self.log.info("All images rejected in BestSeeingSelectVisitsTask.")
669 raise pipeBase.NoWorkFound(f"No good images found for {dataId}")
670 else:
671 self.log.info(
672 "%d images selected with FWHM range of %f--%f arcseconds",
673 len(output),
674 fwhmSizes[visits.index(output[0])],
675 fwhmSizes[visits.index(output[-1])],
676 )
678 # In order to store as a StructuredDataDict, convert list to dict
679 goodVisits = {key: True for key in output}
680 return pipeBase.Struct(goodVisits=goodVisits)
682 def makePatchPolygon(self, skyMap, dataId):
683 """Return True if sky polygon overlaps visit.
685 Parameters
686 ----------
687 skyMap : `lsst.afw.table.ExposureCatalog`
688 Exposure catalog with per-detector geometry.
689 dataId : `dict` of dataId keys
690 For retrieving patch info.
692 Returns
693 -------
694 result : `lsst.sphgeom.ConvexPolygon.convexHull`
695 Polygon of patch's outer bbox.
696 """
697 wcs = skyMap[dataId['tract']].getWcs()
698 bbox = skyMap[dataId['tract']][dataId['patch']].getOuterBBox()
699 sphCorners = wcs.pixelToSky(lsst.geom.Box2D(bbox).getCorners())
700 result = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in sphCorners])
701 return result
703 def doesIntersectPolygon(self, visitSummary, polygon):
704 """Return True if sky polygon overlaps visit.
706 Parameters
707 ----------
708 visitSummary : `lsst.afw.table.ExposureCatalog`
709 Exposure catalog with per-detector geometry.
710 polygon :` lsst.sphgeom.ConvexPolygon.convexHull`
711 Polygon to check overlap.
713 Returns
714 -------
715 doesIntersect : `bool`
716 True if the visit overlaps the polygon.
717 """
718 doesIntersect = False
719 for detectorSummary in visitSummary:
720 if (np.all(np.isfinite(detectorSummary['raCorners']))
721 and np.all(np.isfinite(detectorSummary['decCorners']))):
722 corners = [lsst.geom.SpherePoint(ra, dec, units=lsst.geom.degrees).getVector()
723 for (ra, dec) in
724 zip(detectorSummary['raCorners'], detectorSummary['decCorners'])]
725 detectorPolygon = lsst.sphgeom.ConvexPolygon.convexHull(corners)
726 if detectorPolygon.intersects(polygon):
727 doesIntersect = True
728 break
729 return doesIntersect
732class BestSeeingQuantileSelectVisitsConfig(pipeBase.PipelineTaskConfig,
733 pipelineConnections=BestSeeingSelectVisitsConnections):
734 qMin = pexConfig.RangeField(
735 doc="Lower bound of quantile range to select. Sorts visits by seeing from narrow to wide, "
736 "and select those in the interquantile range (qMin, qMax). Set qMin to 0 for Best Seeing. "
737 "This config should be changed from zero only for exploratory diffIm testing.",
738 dtype=float,
739 default=0,
740 min=0,
741 max=1,
742 )
743 qMax = pexConfig.RangeField(
744 doc="Upper bound of quantile range to select. Sorts visits by seeing from narrow to wide, "
745 "and select those in the interquantile range (qMin, qMax). Set qMax to 1 for Worst Seeing.",
746 dtype=float,
747 default=0.33,
748 min=0,
749 max=1,
750 )
751 nVisitsMin = pexConfig.Field(
752 doc="The minimum number of visits to select, if qMin and qMax alone would have "
753 "selected fewer. In regimes with many visits, at least this number of visits will be "
754 "selected, superceding quantile when necessary. "
755 "For example, if 10 visits cover this patch, qMin=0, qMax=0.33, and nVisitsMin=5, "
756 "the best 5 visits will be selected, even though 5 > 0.33*10. "
757 "In regimes with few visits, all available visits will be selected. "
758 "For example, if 2 visits cover this patch and nVisitsMin=12, "
759 "both visits will be selected regardless of qMin and qMax.",
760 dtype=int,
761 default=6,
762 )
763 doConfirmOverlap = pexConfig.Field(
764 dtype=bool,
765 doc="Do remove visits that do not actually overlap the patch?",
766 default=True,
767 )
768 minMJD = pexConfig.Field(
769 dtype=float,
770 doc="Minimum visit MJD to select",
771 default=None,
772 optional=True
773 )
774 maxMJD = pexConfig.Field(
775 dtype=float,
776 doc="Maximum visit MJD to select",
777 default=None,
778 optional=True
779 )
780 visitSummaryMinValues = pexConfig.DictField(
781 keytype=str,
782 itemtype=float,
783 doc=("Optional thresholds for visit selection. Keys are visit_summary column names"
784 "(e.g. 'effTimeZeroPointScale'), and values are minimum allowed values."),
785 default=None,
786 optional=True
787 )
790class BestSeeingQuantileSelectVisitsTask(BestSeeingSelectVisitsTask):
791 """Select a quantile of the best-seeing visits.
793 Selects the best (for example, third) full visits based on the average
794 PSF width in the entire visit. It can also be used for difference imaging
795 experiments that require templates with the worst seeing visits.
796 For example, selecting the worst third can be acheived by
797 changing the config parameters qMin to 0.66 and qMax to 1.
798 """
799 ConfigClass = BestSeeingQuantileSelectVisitsConfig
800 _DefaultName = 'bestSeeingQuantileSelectVisits'
802 def run(self, visitSummaries, skyMap, dataId):
803 if self.config.doConfirmOverlap:
804 patchPolygon = self.makePatchPolygon(skyMap, dataId)
805 visits = np.array([visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries])
806 radius = np.empty(len(visits))
807 intersects = np.full(len(visits), True)
808 for i, visitSummary in enumerate(visitSummaries):
809 # read in one-by-one and only once. There may be hundreds
810 visitSummary = visitSummary.get()
811 # psfSigma is PSF model determinant radius at chip center in pixels
812 psfSigma = np.nanmedian([vs['psfSigma'] for vs in visitSummary])
813 radius[i] = psfSigma
814 if self.config.doConfirmOverlap:
815 intersects[i] = self.doesIntersectPolygon(visitSummary, patchPolygon)
816 if self.config.minMJD or self.config.maxMJD:
817 # mjd is guaranteed to be the same for every detector in the
818 # visitSummary.
819 mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD)
820 aboveMin = mjd > self.config.minMJD if self.config.minMJD else True
821 belowMax = mjd < self.config.maxMJD if self.config.maxMJD else True
822 intersects[i] = intersects[i] and aboveMin and belowMax
823 if (self.config.visitSummaryMinValues
824 and not _meetsVisitSummaryMinValues(visits[i], visitSummary,
825 self.config.visitSummaryMinValues, self.log)):
826 intersects[i] = False
828 sortedVisits = [v for rad, v in sorted(zip(radius[intersects], visits[intersects]))]
829 lowerBound = min(int(np.round(self.config.qMin*len(visits[intersects]))),
830 max(0, len(visits[intersects]) - self.config.nVisitsMin))
831 upperBound = max(int(np.round(self.config.qMax*len(visits[intersects]))), self.config.nVisitsMin)
833 # In order to store as a StructuredDataDict, convert list to dict
834 goodVisits = {int(visit): True for visit in sortedVisits[lowerBound:upperBound]}
835 return pipeBase.Struct(goodVisits=goodVisits)