LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
selectImages.py
Go to the documentation of this file.
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
37class DatabaseSelectImagesConfig(pexConfig.Config):
38 """Base configuration for subclasses of BaseSelectImagesTask that use a
39 database.
40 """
41
42 host = pexConfig.Field(
43 doc="Database server host name",
44 dtype=str,
45 )
46 port = pexConfig.Field(
47 doc="Database server port",
48 dtype=int,
49 )
50 database = pexConfig.Field(
51 doc="Name of database",
52 dtype=str,
53 )
54 maxExposures = pexConfig.Field(
55 doc="maximum exposures to select; intended for debugging; ignored if None",
56 dtype=int,
57 optional=True,
58 )
59
60
61class BaseExposureInfo(pipeBase.Struct):
62 """Data about a selected exposure.
63
64 Parameters
65 ----------
66 dataId : `dict`
67 Data ID keys of exposure.
68 coordList : `list` [`lsst.afw.geom.SpherePoint`]
69 ICRS coordinates of the corners of the exposure
70 plus any others items that are desired.
71 """
72
73 def __init__(self, dataId, coordList):
74 super(BaseExposureInfo, self).__init__(dataId=dataId, coordList=coordList)
75
76
77class BaseSelectImagesTask(pipeBase.Task):
78 """Base task for selecting images suitable for coaddition.
79 """
80
81 ConfigClass = pexConfig.Config
82 _DefaultName = "selectImages"
83
84 @timeMethod
85 def run(self, coordList):
86 """Select images suitable for coaddition in a particular region.
87
88 Parameters
89 ----------
90 coordList : `list` [`lsst.geom.SpherePoint`] or `None`
91 List of coordinates defining region of interest; if `None`, then
92 select all images subclasses may add additional keyword arguments,
93 as required.
94
95 Returns
96 -------
97 result : `pipeBase.Struct`
98 Results as a struct with attributes:
99
100 ``exposureInfoList``
101 A list of exposure information objects (subclasses of
102 BaseExposureInfo), which have at least the following fields:
103 - dataId: Data ID dictionary (`dict`).
104 - coordList: ICRS coordinates of the corners of the exposure.
105 (`list` [`lsst.geom.SpherePoint`])
106 """
107 raise NotImplementedError()
108
109
110def _extractKeyValue(dataList, keys=None):
111 """Extract the keys and values from a list of dataIds.
112
113 The input dataList is a list of objects that have 'dataId' members.
114 This allows it to be used for both a list of data references and a
115 list of ExposureInfo.
116
117 Parameters
118 ----------
119 dataList : `Unknown`
120 keys : `Unknown`
121
122 Returns
123 -------
124 keys : `Unknown`
125 values : `Unknown`
126
127 Raises
128 ------
129 RuntimeError
130 Raised if DataId keys are inconsistent.
131 """
132 assert len(dataList) > 0
133 if keys is None:
134 keys = sorted(dataList[0].dataId.keys())
135 keySet = set(keys)
136 values = list()
137 for data in dataList:
138 thisKeys = set(data.dataId.keys())
139 if thisKeys != keySet:
140 raise RuntimeError("DataId keys inconsistent: %s vs %s" % (keySet, thisKeys))
141 values.append(tuple(data.dataId[k] for k in keys))
142 return keys, values
143
144
145class SelectStruct(pipeBase.Struct):
146 """A container for data to be passed to the WcsSelectImagesTask.
147
148 Parameters
149 ----------
150 dataRef : `Unknown`
151 Data reference.
152 wcs : `lsst.afw.geom.SkyWcs`
153 Coordinate system definition (wcs).
154 bbox : `lsst.geom.box.Box2I`
155 Integer bounding box for image.
156 """
157
158 def __init__(self, dataRef, wcs, bbox):
159 super(SelectStruct, self).__init__(dataRef=dataRef, wcs=wcs, bbox=bbox)
160
161
163 """Select images using their Wcs.
164
165 We use the "convexHull" method of lsst.sphgeom.ConvexPolygon to define
166 polygons on the celestial sphere, and test the polygon of the
167 patch for overlap with the polygon of the image.
168
169 We use "convexHull" instead of generating a ConvexPolygon
170 directly because the standard for the inputs to ConvexPolygon
171 are pretty high and we don't want to be responsible for reaching them.
172 """
173
174 def run(self, wcsList, bboxList, coordList, dataIds=None, **kwargs):
175 """Return indices of provided lists that meet the selection criteria.
176
177 Parameters
178 ----------
179 wcsList : `list` [`lsst.afw.geom.SkyWcs`]
180 Specifying the WCS's of the input ccds to be selected.
181 bboxList : `list` [`lsst.geom.Box2I`]
182 Specifying the bounding boxes of the input ccds to be selected.
183 coordList : `list` [`lsst.geom.SpherePoint`]
184 ICRS coordinates specifying boundary of the patch.
185 dataIds : iterable [`lsst.daf.butler.dataId`] or `None`, optional
186 An iterable object of dataIds which point to reference catalogs.
187 **kwargs
188 Additional keyword arguments.
189
190 Returns
191 -------
192 result : `list` [`int`]
193 The indices of selected ccds.
194 """
195 if dataIds is None:
196 dataIds = [None] * len(wcsList)
197 patchVertices = [coord.getVector() for coord in coordList]
198 patchPoly = lsst.sphgeom.ConvexPolygon.convexHull(patchVertices)
199 result = []
200 for i, (imageWcs, imageBox, dataId) in enumerate(zip(wcsList, bboxList, dataIds)):
201 if imageWcs is None:
202 self.log.info("De-selecting exposure %s: Exposure has no WCS.", dataId)
203 else:
204 imageCorners = self.getValidImageCorners(imageWcs, imageBox, patchPoly, dataId)
205 if imageCorners:
206 result.append(i)
207 return result
208
209 def getValidImageCorners(self, imageWcs, imageBox, patchPoly, dataId=None):
210 """Return corners or `None` if bad.
211
212 Parameters
213 ----------
214 imageWcs : `Unknown`
215 imageBox : `Unknown`
216 patchPoly : `Unknown`
217 dataId : `Unknown`
218 """
219 try:
220 imageCorners = [imageWcs.pixelToSky(pix) for pix in geom.Box2D(imageBox).getCorners()]
222 # Protecting ourselves from awful Wcs solutions in input images
223 self.log.debug("WCS error in testing calexp %s (%s): deselecting", dataId, e)
224 return None
225
226 imagePoly = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in imageCorners])
227 if imagePoly is None:
228 self.log.debug("Unable to create polygon from image %s: deselecting", dataId)
229 return None
230
231 if patchPoly.intersects(imagePoly):
232 # "intersects" also covers "contains" or "is contained by"
233 self.log.info("Selecting calexp %s", dataId)
234 return imageCorners
235
236 return None
237
238
239class PsfWcsSelectImagesConnections(pipeBase.PipelineTaskConnections,
240 dimensions=("tract", "patch", "skymap", "instrument", "visit"),
241 defaultTemplates={"coaddName": "deep"}):
242 pass
243
244
245class PsfWcsSelectImagesConfig(pipeBase.PipelineTaskConfig,
246 pipelineConnections=PsfWcsSelectImagesConnections):
247 maxEllipResidual = pexConfig.Field(
248 doc="Maximum median ellipticity residual",
249 dtype=float,
250 default=0.007,
251 optional=True,
252 )
253 maxSizeScatter = pexConfig.Field(
254 doc="Maximum scatter in the size residuals",
255 dtype=float,
256 optional=True,
257 )
258 maxScaledSizeScatter = pexConfig.Field(
259 doc="Maximum scatter in the size residuals, scaled by the median size",
260 dtype=float,
261 default=0.019,
262 optional=True,
263 )
264 maxPsfTraceRadiusDelta = pexConfig.Field(
265 doc="Maximum delta (max - min) of model PSF trace radius values evaluated on a grid on "
266 "the unmasked detector pixels (pixel).",
267 dtype=float,
268 default=0.7,
269 optional=True,
270 )
271
272
273class PsfWcsSelectImagesTask(WcsSelectImagesTask):
274 """Select images using their Wcs and cuts on the PSF properties.
275
276 The PSF quality criteria are based on the size and ellipticity
277 residuals from the adaptive second moments of the star and the PSF.
278
279 The criteria are:
280 - the median of the ellipticty residuals.
281 - the robust scatter of the size residuals (using the median absolute
282 deviation).
283 - the robust scatter of the size residuals scaled by the square of
284 the median size.
285 """
286
287 ConfigClass = PsfWcsSelectImagesConfig
288 _DefaultName = "PsfWcsSelectImages"
289
290 def run(self, wcsList, bboxList, coordList, visitSummary, dataIds=None, **kwargs):
291 """Return indices of provided lists that meet the selection criteria.
292
293 Parameters
294 ----------
295 wcsList : `list` [`lsst.afw.geom.SkyWcs`]
296 Specifying the WCS's of the input ccds to be selected.
297 bboxList : `list` [`lsst.geom.Box2I`]
298 Specifying the bounding boxes of the input ccds to be selected.
299 coordList : `list` [`lsst.geom.SpherePoint`]
300 ICRS coordinates specifying boundary of the patch.
301 visitSummary : `list` [`lsst.afw.table.ExposureCatalog`]
302 containing the PSF shape information for the input ccds to be
303 selected.
304 dataIds : iterable [`lsst.daf.butler.dataId`] or `None`, optional
305 An iterable object of dataIds which point to reference catalogs.
306 **kwargs
307 Additional keyword arguments.
308
309 Returns
310 -------
311 goodPsf : `list` [`int`]
312 The indices of selected ccds.
313 """
314 goodWcs = super(PsfWcsSelectImagesTask, self).run(wcsList=wcsList, bboxList=bboxList,
315 coordList=coordList, dataIds=dataIds)
316
317 goodPsf = []
318
319 for i, dataId in enumerate(dataIds):
320 if i not in goodWcs:
321 continue
322 if self.isValid(visitSummary, dataId["detector"]):
323 goodPsf.append(i)
324
325 return goodPsf
326
327 def isValid(self, visitSummary, detectorId):
328 """Should this ccd be selected based on its PSF shape information.
329
330 Parameters
331 ----------
332 visitSummary : `lsst.afw.table.ExposureCatalog`
333 Exposure catalog with per-detector summary information.
334 detectorId : `int`
335 Detector identifier.
336
337 Returns
338 -------
339 valid : `bool`
340 True if selected.
341 """
342 row = visitSummary.find(detectorId)
343 if row is None:
344 # This is not listed, so it must be bad.
345 self.log.warning("Removing detector %d because summary stats not available.", detectorId)
346 return False
347
348 medianE = np.sqrt(row["psfStarDeltaE1Median"]**2. + row["psfStarDeltaE2Median"]**2.)
349 scatterSize = row["psfStarDeltaSizeScatter"]
350 scaledScatterSize = row["psfStarScaledDeltaSizeScatter"]
351 psfTraceRadiusDelta = row["psfTraceRadiusDelta"]
352
353 valid = True
354 if self.config.maxEllipResidual and not (medianE <= self.config.maxEllipResidual):
355 self.log.info("Removing visit %d detector %d because median e residual too large: %f vs %f",
356 row["visit"], detectorId, medianE, self.config.maxEllipResidual)
357 valid = False
358 elif self.config.maxSizeScatter and not (scatterSize <= self.config.maxSizeScatter):
359 self.log.info("Removing visit %d detector %d because size scatter too large: %f vs %f",
360 row["visit"], detectorId, scatterSize, self.config.maxSizeScatter)
361 valid = False
362 elif self.config.maxScaledSizeScatter and not (scaledScatterSize <= self.config.maxScaledSizeScatter):
363 self.log.info("Removing visit %d detector %d because scaled size scatter too large: %f vs %f",
364 row["visit"], detectorId, scaledScatterSize, self.config.maxScaledSizeScatter)
365 valid = False
366 elif (
367 self.config.maxPsfTraceRadiusDelta is not None
368 and not (psfTraceRadiusDelta <= self.config.maxPsfTraceRadiusDelta)
369 ):
370 self.log.info(
371 "Removing visit %d detector %d because max-min delta of model PSF trace radius values "
372 "across the unmasked detector pixels is not finite or too large: %.3f vs %.3f (pixels)",
373 row["visit"], detectorId, psfTraceRadiusDelta, self.config.maxPsfTraceRadiusDelta
374 )
375 valid = False
376
377 return valid
378
379
380class BestSeeingSelectVisitsConnections(pipeBase.PipelineTaskConnections,
381 dimensions=("tract", "patch", "skymap", "band", "instrument"),
382 defaultTemplates={"coaddName": "goodSeeing",
383 "calexpType": ""}):
384 skyMap = pipeBase.connectionTypes.Input(
385 doc="Input definition of geometry/bbox and projection/wcs for coadded exposures",
386 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
387 storageClass="SkyMap",
388 dimensions=("skymap",),
389 )
390 visitSummaries = pipeBase.connectionTypes.Input(
391 doc="Per-visit consolidated exposure metadata",
392 name="finalVisitSummary",
393 storageClass="ExposureCatalog",
394 dimensions=("instrument", "visit",),
395 multiple=True,
396 deferLoad=True
397 )
398 goodVisits = pipeBase.connectionTypes.Output(
399 doc="Selected visits to be coadded.",
400 name="{coaddName}Visits",
401 storageClass="StructuredDataDict",
402 dimensions=("instrument", "tract", "patch", "skymap", "band"),
403 )
404
405
406class BestSeeingSelectVisitsConfig(pipeBase.PipelineTaskConfig,
407 pipelineConnections=BestSeeingSelectVisitsConnections):
408 nVisitsMax = pexConfig.RangeField(
409 dtype=int,
410 doc="Maximum number of visits to select; use -1 to select all.",
411 default=12,
412 min=-1,
413 )
414 maxPsfFwhm = pexConfig.Field(
415 dtype=float,
416 doc="Maximum PSF FWHM (in arcseconds) to select",
417 default=1.5,
418 optional=True
419 )
420 minPsfFwhm = pexConfig.Field(
421 dtype=float,
422 doc="Minimum PSF FWHM (in arcseconds) to select",
423 default=0.,
424 optional=True
425 )
426 doConfirmOverlap = pexConfig.Field(
427 dtype=bool,
428 doc="Do remove visits that do not actually overlap the patch?",
429 default=True,
430 )
431 minMJD = pexConfig.Field(
432 dtype=float,
433 doc="Minimum visit MJD to select",
434 default=None,
435 optional=True
436 )
437 maxMJD = pexConfig.Field(
438 dtype=float,
439 doc="Maximum visit MJD to select",
440 default=None,
441 optional=True
442 )
443
444
445class BestSeeingSelectVisitsTask(pipeBase.PipelineTask):
446 """Select up to a maximum number of the best-seeing visits.
447
448 Don't exceed the FWHM range specified by configs min(max)PsfFwhm.
449 This Task is a port of the Gen2 image-selector used in the AP pipeline:
450 BestSeeingSelectImagesTask. This Task selects full visits based on the
451 average PSF of the entire visit.
452 """
453
454 ConfigClass = BestSeeingSelectVisitsConfig
455 _DefaultName = 'bestSeeingSelectVisits'
456
457 def runQuantum(self, butlerQC, inputRefs, outputRefs):
458 inputs = butlerQC.get(inputRefs)
459 quantumDataId = butlerQC.quantum.dataId
460 outputs = self.run(**inputs, dataId=quantumDataId)
461 butlerQC.put(outputs, outputRefs)
462
463 def run(self, visitSummaries, skyMap, dataId):
464 """Run task.
465
466 Parameters
467 ----------
468 visitSummary : `list` [`lsst.pipe.base.connections.DeferredDatasetRef`]
469 List of `lsst.pipe.base.connections.DeferredDatasetRef` of
470 visitSummary tables of type `lsst.afw.table.ExposureCatalog`.
471 skyMap : `lsst.skyMap.SkyMap`
472 SkyMap for checking visits overlap patch.
473 dataId : `dict` of dataId keys
474 For retrieving patch info for checking visits overlap patch.
475
476 Returns
477 -------
478 result : `lsst.pipe.base.Struct`
479 Results as a struct with attributes:
480
481 ``goodVisits``
482 A `dict` with selected visit ids as keys,
483 so that it can be be saved as a StructuredDataDict.
484 StructuredDataList's are currently limited.
485 """
486 if self.config.doConfirmOverlap:
487 patchPolygon = self.makePatchPolygon(skyMap, dataId)
488
489 inputVisits = [visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries]
490 fwhmSizes = []
491 visits = []
492 for visit, visitSummary in zip(inputVisits, visitSummaries):
493 # read in one-by-one and only once. There may be hundreds
494 visitSummary = visitSummary.get()
495
496 # mjd is guaranteed to be the same for every detector in the
497 # visitSummary.
498 mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD)
499
500 pixToArcseconds = [vs.getWcs().getPixelScale(vs.getBBox().getCenter()).asArcseconds()
501 for vs in visitSummary if vs.getWcs()]
502 # psfSigma is PSF model determinant radius at chip center in pixels
503 psfSigmas = np.array([vs['psfSigma'] for vs in visitSummary if vs.getWcs()])
504 fwhm = np.nanmean(psfSigmas * pixToArcseconds) * np.sqrt(8.*np.log(2.))
505
506 if self.config.maxPsfFwhm and fwhm > self.config.maxPsfFwhm:
507 continue
508 if self.config.minPsfFwhm and fwhm < self.config.minPsfFwhm:
509 continue
510 if self.config.minMJD and mjd < self.config.minMJD:
511 self.log.debug('MJD %f earlier than %.2f; rejecting', mjd, self.config.minMJD)
512 continue
513 if self.config.maxMJD and mjd > self.config.maxMJD:
514 self.log.debug('MJD %f later than %.2f; rejecting', mjd, self.config.maxMJD)
515 continue
516 if self.config.doConfirmOverlap and not self.doesIntersectPolygon(visitSummary, patchPolygon):
517 continue
518
519 fwhmSizes.append(fwhm)
520 visits.append(visit)
521
522 sortedVisits = [ind for (_, ind) in sorted(zip(fwhmSizes, visits))]
523 if self.config.nVisitsMax < 0:
524 output = sortedVisits
525 else:
526 output = sortedVisits[:self.config.nVisitsMax]
527
528 if len(output) == 0:
529 self.log.info("All images rejected in BestSeeingSelectVisitsTask.")
530 raise pipeBase.NoWorkFound(f"No good images found for {dataId}")
531 else:
532 self.log.info(
533 "%d images selected with FWHM range of %f--%f arcseconds",
534 len(output),
535 fwhmSizes[visits.index(output[0])],
536 fwhmSizes[visits.index(output[-1])],
537 )
538
539 # In order to store as a StructuredDataDict, convert list to dict
540 goodVisits = {key: True for key in output}
541 return pipeBase.Struct(goodVisits=goodVisits)
542
543 def makePatchPolygon(self, skyMap, dataId):
544 """Return True if sky polygon overlaps visit.
545
546 Parameters
547 ----------
548 skyMap : `lsst.afw.table.ExposureCatalog`
549 Exposure catalog with per-detector geometry.
550 dataId : `dict` of dataId keys
551 For retrieving patch info.
552
553 Returns
554 -------
555 result : `lsst.sphgeom.ConvexPolygon.convexHull`
556 Polygon of patch's outer bbox.
557 """
558 wcs = skyMap[dataId['tract']].getWcs()
559 bbox = skyMap[dataId['tract']][dataId['patch']].getOuterBBox()
560 sphCorners = wcs.pixelToSky(lsst.geom.Box2D(bbox).getCorners())
561 result = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in sphCorners])
562 return result
563
564 def doesIntersectPolygon(self, visitSummary, polygon):
565 """Return True if sky polygon overlaps visit.
566
567 Parameters
568 ----------
569 visitSummary : `lsst.afw.table.ExposureCatalog`
570 Exposure catalog with per-detector geometry.
571 polygon :` lsst.sphgeom.ConvexPolygon.convexHull`
572 Polygon to check overlap.
573
574 Returns
575 -------
576 doesIntersect : `bool`
577 True if the visit overlaps the polygon.
578 """
579 doesIntersect = False
580 for detectorSummary in visitSummary:
581 if (np.all(np.isfinite(detectorSummary['raCorners']))
582 and np.all(np.isfinite(detectorSummary['decCorners']))):
583 corners = [lsst.geom.SpherePoint(ra, dec, units=lsst.geom.degrees).getVector()
584 for (ra, dec) in
585 zip(detectorSummary['raCorners'], detectorSummary['decCorners'])]
586 detectorPolygon = lsst.sphgeom.ConvexPolygon.convexHull(corners)
587 if detectorPolygon.intersects(polygon):
588 doesIntersect = True
589 break
590 return doesIntersect
591
592
593class BestSeeingQuantileSelectVisitsConfig(pipeBase.PipelineTaskConfig,
594 pipelineConnections=BestSeeingSelectVisitsConnections):
595 qMin = pexConfig.RangeField(
596 doc="Lower bound of quantile range to select. Sorts visits by seeing from narrow to wide, "
597 "and select those in the interquantile range (qMin, qMax). Set qMin to 0 for Best Seeing. "
598 "This config should be changed from zero only for exploratory diffIm testing.",
599 dtype=float,
600 default=0,
601 min=0,
602 max=1,
603 )
604 qMax = pexConfig.RangeField(
605 doc="Upper bound of quantile range to select. Sorts visits by seeing from narrow to wide, "
606 "and select those in the interquantile range (qMin, qMax). Set qMax to 1 for Worst Seeing.",
607 dtype=float,
608 default=0.33,
609 min=0,
610 max=1,
611 )
612 nVisitsMin = pexConfig.Field(
613 doc="At least this number of visits selected and supercedes quantile. For example, if 10 visits "
614 "cover this patch, qMin=0.33, and nVisitsMin=5, the best 5 visits will be selected.",
615 dtype=int,
616 default=6,
617 )
618 doConfirmOverlap = pexConfig.Field(
619 dtype=bool,
620 doc="Do remove visits that do not actually overlap the patch?",
621 default=True,
622 )
623 minMJD = pexConfig.Field(
624 dtype=float,
625 doc="Minimum visit MJD to select",
626 default=None,
627 optional=True
628 )
629 maxMJD = pexConfig.Field(
630 dtype=float,
631 doc="Maximum visit MJD to select",
632 default=None,
633 optional=True
634 )
635
636
637class BestSeeingQuantileSelectVisitsTask(BestSeeingSelectVisitsTask):
638 """Select a quantile of the best-seeing visits.
639
640 Selects the best (for example, third) full visits based on the average
641 PSF width in the entire visit. It can also be used for difference imaging
642 experiments that require templates with the worst seeing visits.
643 For example, selecting the worst third can be acheived by
644 changing the config parameters qMin to 0.66 and qMax to 1.
645 """
646 ConfigClass = BestSeeingQuantileSelectVisitsConfig
647 _DefaultName = 'bestSeeingQuantileSelectVisits'
648
649 def run(self, visitSummaries, skyMap, dataId):
650 if self.config.doConfirmOverlap:
651 patchPolygon = self.makePatchPolygon(skyMap, dataId)
652 visits = np.array([visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries])
653 radius = np.empty(len(visits))
654 intersects = np.full(len(visits), True)
655 for i, visitSummary in enumerate(visitSummaries):
656 # read in one-by-one and only once. There may be hundreds
657 visitSummary = visitSummary.get()
658 # psfSigma is PSF model determinant radius at chip center in pixels
659 psfSigma = np.nanmedian([vs['psfSigma'] for vs in visitSummary])
660 radius[i] = psfSigma
661 if self.config.doConfirmOverlap:
662 intersects[i] = self.doesIntersectPolygon(visitSummary, patchPolygon)
663 if self.config.minMJD or self.config.maxMJD:
664 # mjd is guaranteed to be the same for every detector in the
665 # visitSummary.
666 mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD)
667 aboveMin = mjd > self.config.minMJD if self.config.minMJD else True
668 belowMax = mjd < self.config.maxMJD if self.config.maxMJD else True
669 intersects[i] = intersects[i] and aboveMin and belowMax
670
671 sortedVisits = [v for rad, v in sorted(zip(radius[intersects], visits[intersects]))]
672 lowerBound = min(int(np.round(self.config.qMin*len(visits[intersects]))),
673 max(0, len(visits[intersects]) - self.config.nVisitsMin))
674 upperBound = max(int(np.round(self.config.qMax*len(visits[intersects]))), self.config.nVisitsMin)
675
676 # In order to store as a StructuredDataDict, convert list to dict
677 goodVisits = {int(visit): True for visit in sortedVisits[lowerBound:upperBound]}
678 return pipeBase.Struct(goodVisits=goodVisits)
int min
int max
A floating-point coordinate rectangle geometry.
Definition Box.h:413
Point in an unspecified spherical coordinate system.
Definition SpherePoint.h:57
Reports arguments outside the domain of an operation.
Definition Runtime.h:57
Reports errors that are due to events beyond the control of the program.
Definition Runtime.h:104
run(self, wcsList, bboxList, coordList, dataIds=None, **kwargs)
getValidImageCorners(self, imageWcs, imageBox, patchPoly, dataId=None)
static ConvexPolygon convexHull(std::vector< UnitVector3d > const &points)
convexHull returns the convex hull of the given set of points if it exists and throws an exception ot...
daf::base::PropertySet * set
Definition fits.cc:931
bool isValid
Definition fits.cc:404
_extractKeyValue(dataList, keys=None)