LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
selectImages.py
Go to the documentation of this file.
1 #
2 # LSST Data Management System
3 # Copyright 2008, 2009, 2010 LSST Corporation.
4 #
5 # This product includes software developed by the
6 # LSST Project (http://www.lsst.org/).
7 #
8 # This program is free software: you can redistribute it and/or modify
9 # it under the terms of the GNU General Public License as published by
10 # the Free Software Foundation, either version 3 of the License, or
11 # (at your option) any later version.
12 #
13 # This program is distributed in the hope that it will be useful,
14 # but WITHOUT ANY WARRANTY; without even the implied warranty of
15 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 # GNU General Public License for more details.
17 #
18 # You should have received a copy of the LSST License Statement and
19 # the GNU General Public License along with this program. If not,
20 # see <http://www.lsstcorp.org/LegalNotices/>.
21 #
22 import numpy as np
23 import lsst.sphgeom
24 import lsst.utils as utils
25 import lsst.pex.config as pexConfig
26 import lsst.pex.exceptions as pexExceptions
27 import lsst.geom as geom
28 import lsst.pipe.base as pipeBase
29 from lsst.skymap import BaseSkyMap
30 from lsst.daf.base import DateTime
31 from lsst.utils.timer import timeMethod
32 
33 __all__ = ["BaseSelectImagesTask", "BaseExposureInfo", "WcsSelectImagesTask", "PsfWcsSelectImagesTask",
34  "DatabaseSelectImagesConfig", "BestSeeingWcsSelectImagesTask", "BestSeeingSelectVisitsTask",
35  "BestSeeingQuantileSelectVisitsTask"]
36 
37 
38 class DatabaseSelectImagesConfig(pexConfig.Config):
39  """Base configuration for subclasses of BaseSelectImagesTask that use a database"""
40  host = pexConfig.Field(
41  doc="Database server host name",
42  dtype=str,
43  )
44  port = pexConfig.Field(
45  doc="Database server port",
46  dtype=int,
47  )
48  database = pexConfig.Field(
49  doc="Name of database",
50  dtype=str,
51  )
52  maxExposures = pexConfig.Field(
53  doc="maximum exposures to select; intended for debugging; ignored if None",
54  dtype=int,
55  optional=True,
56  )
57 
58 
59 class BaseExposureInfo(pipeBase.Struct):
60  """Data about a selected exposure
61  """
62 
63  def __init__(self, dataId, coordList):
64  """Create exposure information that can be used to generate data references
65 
66  The object has the following fields:
67  - dataId: data ID of exposure (a dict)
68  - coordList: ICRS coordinates of the corners of the exposure (list of lsst.geom.SpherePoint)
69  plus any others items that are desired
70  """
71  super(BaseExposureInfo, self).__init__(dataId=dataId, coordList=coordList)
72 
73 
74 class BaseSelectImagesTask(pipeBase.Task):
75  """Base task for selecting images suitable for coaddition
76  """
77  ConfigClass = pexConfig.Config
78  _DefaultName = "selectImages"
79 
80  @timeMethod
81  def run(self, coordList):
82  """Select images suitable for coaddition in a particular region
83 
84  @param[in] coordList: list of coordinates defining region of interest; if None then select all images
85  subclasses may add additional keyword arguments, as required
86 
87  @return a pipeBase Struct containing:
88  - exposureInfoList: a list of exposure information objects (subclasses of BaseExposureInfo),
89  which have at least the following fields:
90  - dataId: data ID dictionary
91  - coordList: ICRS coordinates of the corners of the exposure (list of lsst.geom.SpherePoint)
92  """
93  raise NotImplementedError()
94 
95  def _runArgDictFromDataId(self, dataId):
96  """Extract keyword arguments for run (other than coordList) from a data ID
97 
98  @return keyword arguments for run (other than coordList), as a dict
99  """
100  raise NotImplementedError()
101 
102  def runDataRef(self, dataRef, coordList, makeDataRefList=True, selectDataList=[]):
103  """Run based on a data reference
104 
105  This delegates to run() and _runArgDictFromDataId() to do the actual
106  selection. In the event that the selectDataList is non-empty, this will
107  be used to further restrict the selection, providing the user with
108  additional control over the selection.
109 
110  @param[in] dataRef: data reference; must contain any extra keys needed by the subclass
111  @param[in] coordList: list of coordinates defining region of interest; if None, search the whole sky
112  @param[in] makeDataRefList: if True, return dataRefList
113  @param[in] selectDataList: List of SelectStruct with dataRefs to consider for selection
114  @return a pipeBase Struct containing:
115  - exposureInfoList: a list of objects derived from ExposureInfo
116  - dataRefList: a list of data references (None if makeDataRefList False)
117  """
118  runArgDict = self._runArgDictFromDataId_runArgDictFromDataId(dataRef.dataId)
119  exposureInfoList = self.runrun(coordList, **runArgDict).exposureInfoList
120 
121  if len(selectDataList) > 0 and len(exposureInfoList) > 0:
122  # Restrict the exposure selection further
123  ccdKeys, ccdValues = _extractKeyValue(exposureInfoList)
124  inKeys, inValues = _extractKeyValue([s.dataRef for s in selectDataList], keys=ccdKeys)
125  inValues = set(inValues)
126  newExposureInfoList = []
127  for info, ccdVal in zip(exposureInfoList, ccdValues):
128  if ccdVal in inValues:
129  newExposureInfoList.append(info)
130  else:
131  self.log.info("De-selecting exposure %s: not in selectDataList", info.dataId)
132  exposureInfoList = newExposureInfoList
133 
134  if makeDataRefList:
135  butler = dataRef.butlerSubset.butler
136  dataRefList = [butler.dataRef(datasetType="calexp",
137  dataId=expInfo.dataId,
138  ) for expInfo in exposureInfoList]
139  else:
140  dataRefList = None
141 
142  return pipeBase.Struct(
143  dataRefList=dataRefList,
144  exposureInfoList=exposureInfoList,
145  )
146 
147 
148 def _extractKeyValue(dataList, keys=None):
149  """Extract the keys and values from a list of dataIds
150 
151  The input dataList is a list of objects that have 'dataId' members.
152  This allows it to be used for both a list of data references and a
153  list of ExposureInfo
154  """
155  assert len(dataList) > 0
156  if keys is None:
157  keys = sorted(dataList[0].dataId.keys())
158  keySet = set(keys)
159  values = list()
160  for data in dataList:
161  thisKeys = set(data.dataId.keys())
162  if thisKeys != keySet:
163  raise RuntimeError("DataId keys inconsistent: %s vs %s" % (keySet, thisKeys))
164  values.append(tuple(data.dataId[k] for k in keys))
165  return keys, values
166 
167 
168 class SelectStruct(pipeBase.Struct):
169  """A container for data to be passed to the WcsSelectImagesTask"""
170 
171  def __init__(self, dataRef, wcs, bbox):
172  super(SelectStruct, self).__init__(dataRef=dataRef, wcs=wcs, bbox=bbox)
173 
174 
176  """Select images using their Wcs
177 
178  We use the "convexHull" method of lsst.sphgeom.ConvexPolygon to define
179  polygons on the celestial sphere, and test the polygon of the
180  patch for overlap with the polygon of the image.
181 
182  We use "convexHull" instead of generating a ConvexPolygon
183  directly because the standard for the inputs to ConvexPolygon
184  are pretty high and we don't want to be responsible for reaching them.
185  """
186 
187  def runDataRef(self, dataRef, coordList, makeDataRefList=True, selectDataList=[]):
188  """Select images in the selectDataList that overlap the patch
189 
190  This method is the old entry point for the Gen2 commandline tasks and drivers
191  Will be deprecated in v22.
192 
193  @param dataRef: Data reference for coadd/tempExp (with tract, patch)
194  @param coordList: List of ICRS coordinates (lsst.geom.SpherePoint) specifying boundary of patch
195  @param makeDataRefList: Construct a list of data references?
196  @param selectDataList: List of SelectStruct, to consider for selection
197  """
198  dataRefList = []
199  exposureInfoList = []
200 
201  patchVertices = [coord.getVector() for coord in coordList]
202  patchPoly = lsst.sphgeom.ConvexPolygon.convexHull(patchVertices)
203 
204  for data in selectDataList:
205  dataRef = data.dataRef
206  imageWcs = data.wcs
207  imageBox = data.bbox
208 
209  imageCorners = self.getValidImageCornersgetValidImageCorners(imageWcs, imageBox, patchPoly, dataId=None)
210  if imageCorners:
211  dataRefList.append(dataRef)
212  exposureInfoList.append(BaseExposureInfo(dataRef.dataId, imageCorners))
213 
214  return pipeBase.Struct(
215  dataRefList=dataRefList if makeDataRefList else None,
216  exposureInfoList=exposureInfoList,
217  )
218 
219  def run(self, wcsList, bboxList, coordList, dataIds=None, **kwargs):
220  """Return indices of provided lists that meet the selection criteria
221 
222  Parameters:
223  -----------
224  wcsList : `list` of `lsst.afw.geom.SkyWcs`
225  specifying the WCS's of the input ccds to be selected
226  bboxList : `list` of `lsst.geom.Box2I`
227  specifying the bounding boxes of the input ccds to be selected
228  coordList : `list` of `lsst.geom.SpherePoint`
229  ICRS coordinates specifying boundary of the patch.
230 
231  Returns:
232  --------
233  result: `list` of `int`
234  of indices of selected ccds
235  """
236  if dataIds is None:
237  dataIds = [None] * len(wcsList)
238  patchVertices = [coord.getVector() for coord in coordList]
239  patchPoly = lsst.sphgeom.ConvexPolygon.convexHull(patchVertices)
240  result = []
241  for i, (imageWcs, imageBox, dataId) in enumerate(zip(wcsList, bboxList, dataIds)):
242  imageCorners = self.getValidImageCornersgetValidImageCorners(imageWcs, imageBox, patchPoly, dataId)
243  if imageCorners:
244  result.append(i)
245  return result
246 
247  def getValidImageCorners(self, imageWcs, imageBox, patchPoly, dataId=None):
248  "Return corners or None if bad"
249  try:
250  imageCorners = [imageWcs.pixelToSky(pix) for pix in geom.Box2D(imageBox).getCorners()]
252  # Protecting ourselves from awful Wcs solutions in input images
253  self.log.debug("WCS error in testing calexp %s (%s): deselecting", dataId, e)
254  return
255 
256  imagePoly = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in imageCorners])
257  if imagePoly is None:
258  self.log.debug("Unable to create polygon from image %s: deselecting", dataId)
259  return
260 
261  if patchPoly.intersects(imagePoly):
262  # "intersects" also covers "contains" or "is contained by"
263  self.log.info("Selecting calexp %s", dataId)
264  return imageCorners
265 
266 
267 def sigmaMad(array):
268  "Return median absolute deviation scaled to normally distributed data"
269  return 1.4826*np.median(np.abs(array - np.median(array)))
270 
271 
272 class PsfWcsSelectImagesConnections(pipeBase.PipelineTaskConnections,
273  dimensions=("tract", "patch", "skymap", "instrument", "visit"),
274  defaultTemplates={"coaddName": "deep"}):
275  pass
276 
277 
278 class PsfWcsSelectImagesConfig(pipeBase.PipelineTaskConfig,
279  pipelineConnections=PsfWcsSelectImagesConnections):
280  maxEllipResidual = pexConfig.Field(
281  doc="Maximum median ellipticity residual",
282  dtype=float,
283  default=0.007,
284  optional=True,
285  )
286  maxSizeScatter = pexConfig.Field(
287  doc="Maximum scatter in the size residuals",
288  dtype=float,
289  optional=True,
290  )
291  maxScaledSizeScatter = pexConfig.Field(
292  doc="Maximum scatter in the size residuals, scaled by the median size",
293  dtype=float,
294  default=0.009,
295  optional=True,
296  )
297  starSelection = pexConfig.Field(
298  doc="select star with this field",
299  dtype=str,
300  default='calib_psf_used',
301  deprecated=('This field has been moved to ComputeExposureSummaryStatsTask and '
302  'will be removed after v24.')
303  )
304  starShape = pexConfig.Field(
305  doc="name of star shape",
306  dtype=str,
307  default='base_SdssShape',
308  deprecated=('This field has been moved to ComputeExposureSummaryStatsTask and '
309  'will be removed after v24.')
310  )
311  psfShape = pexConfig.Field(
312  doc="name of psf shape",
313  dtype=str,
314  default='base_SdssShape_psf',
315  deprecated=('This field has been moved to ComputeExposureSummaryStatsTask and '
316  'will be removed after v24.')
317  )
318 
319 
320 class PsfWcsSelectImagesTask(WcsSelectImagesTask):
321  """Select images using their Wcs and cuts on the PSF properties
322 
323  The PSF quality criteria are based on the size and ellipticity residuals from the
324  adaptive second moments of the star and the PSF.
325 
326  The criteria are:
327  - the median of the ellipticty residuals
328  - the robust scatter of the size residuals (using the median absolute deviation)
329  - the robust scatter of the size residuals scaled by the square of
330  the median size
331  """
332 
333  ConfigClass = PsfWcsSelectImagesConfig
334  _DefaultName = "PsfWcsSelectImages"
335 
336  def runDataRef(self, dataRef, coordList, makeDataRefList=True, selectDataList=[]):
337  """Select images in the selectDataList that overlap the patch and satisfy PSF quality critera.
338 
339  This method is the old entry point for the Gen2 commandline tasks and drivers
340  Will be deprecated in v22.
341 
342  @param dataRef: Data reference for coadd/tempExp (with tract, patch)
343  @param coordList: List of ICRS coordinates (lsst.geom.SpherePoint) specifying boundary of patch
344  @param makeDataRefList: Construct a list of data references?
345  @param selectDataList: List of SelectStruct, to consider for selection
346  """
347  result = super(PsfWcsSelectImagesTask, self).runDataRef(dataRef, coordList, makeDataRefList,
348  selectDataList)
349 
350  dataRefList = []
351  exposureInfoList = []
352  for dataRef, exposureInfo in zip(result.dataRefList, result.exposureInfoList):
353  butler = dataRef.butlerSubset.butler
354  srcCatalog = butler.get('src', dataRef.dataId)
355  valid = self.isValidGen2(srcCatalog, dataRef.dataId)
356  if valid is False:
357  continue
358 
359  dataRefList.append(dataRef)
360  exposureInfoList.append(exposureInfo)
361 
362  return pipeBase.Struct(
363  dataRefList=dataRefList,
364  exposureInfoList=exposureInfoList,
365  )
366 
367  def run(self, wcsList, bboxList, coordList, visitSummary, dataIds=None, **kwargs):
368  """Return indices of provided lists that meet the selection criteria
369 
370  Parameters:
371  -----------
372  wcsList : `list` of `lsst.afw.geom.SkyWcs`
373  specifying the WCS's of the input ccds to be selected
374  bboxList : `list` of `lsst.geom.Box2I`
375  specifying the bounding boxes of the input ccds to be selected
376  coordList : `list` of `lsst.geom.SpherePoint`
377  ICRS coordinates specifying boundary of the patch.
378  visitSummary : `list` of `lsst.afw.table.ExposureCatalog`
379  containing the PSF shape information for the input ccds to be selected
380 
381  Returns:
382  --------
383  goodPsf: `list` of `int`
384  of indices of selected ccds
385  """
386  goodWcs = super(PsfWcsSelectImagesTask, self).run(wcsList=wcsList, bboxList=bboxList,
387  coordList=coordList, dataIds=dataIds)
388 
389  goodPsf = []
390 
391  for i, dataId in enumerate(dataIds):
392  if i not in goodWcs:
393  continue
394  if self.isValid(visitSummary, dataId['detector']):
395  goodPsf.append(i)
396 
397  return goodPsf
398 
399  def isValid(self, visitSummary, detectorId):
400  """Should this ccd be selected based on its PSF shape information.
401 
402  Parameters
403  ----------
404  visitSummary : `lsst.afw.table.ExposureCatalog`
405  detectorId : `int`
406  Detector identifier.
407 
408  Returns
409  -------
410  valid : `bool`
411  True if selected.
412  """
413  row = visitSummary.find(detectorId)
414  if row is None:
415  # This is not listed, so it must be bad.
416  self.log.warning("Removing visit %d detector %d because summary stats not available.",
417  row["visit"], detectorId)
418  return False
419 
420  medianE = np.sqrt(row["psfStarDeltaE1Median"]**2. + row["psfStarDeltaE2Median"]**2.)
421  scatterSize = row["psfStarDeltaSizeScatter"]
422  scaledScatterSize = row["psfStarScaledDeltaSizeScatter"]
423 
424  valid = True
425  if self.config.maxEllipResidual and medianE > self.config.maxEllipResidual:
426  self.log.info("Removing visit %d detector %d because median e residual too large: %f vs %f",
427  row["visit"], detectorId, medianE, self.config.maxEllipResidual)
428  valid = False
429  elif self.config.maxSizeScatter and scatterSize > self.config.maxSizeScatter:
430  self.log.info("Removing visit %d detector %d because size scatter too large: %f vs %f",
431  row["visit"], detectorId, scatterSize, self.config.maxSizeScatter)
432  valid = False
433  elif self.config.maxScaledSizeScatter and scaledScatterSize > self.config.maxScaledSizeScatter:
434  self.log.info("Removing visit %d detector %d because scaled size scatter too large: %f vs %f",
435  row["visit"], detectorId, scaledScatterSize, self.config.maxScaledSizeScatter)
436  valid = False
437 
438  return valid
439 
440  def isValidGen2(self, srcCatalog, dataId=None):
441  """Should this ccd be selected based on its PSF shape information.
442 
443  This routine is only used in Gen2 processing, and can be
444  removed when Gen2 is retired.
445 
446  Parameters
447  ----------
448  srcCatalog : `lsst.afw.table.SourceCatalog`
449  dataId : `dict` of dataId keys, optional.
450  Used only for logging. Defaults to None.
451 
452  Returns
453  -------
454  valid : `bool`
455  True if selected.
456  """
457  mask = srcCatalog[self.config.starSelection]
458 
459  starXX = srcCatalog[self.config.starShape+'_xx'][mask]
460  starYY = srcCatalog[self.config.starShape+'_yy'][mask]
461  starXY = srcCatalog[self.config.starShape+'_xy'][mask]
462  psfXX = srcCatalog[self.config.psfShape+'_xx'][mask]
463  psfYY = srcCatalog[self.config.psfShape+'_yy'][mask]
464  psfXY = srcCatalog[self.config.psfShape+'_xy'][mask]
465 
466  starSize = np.power(starXX*starYY - starXY**2, 0.25)
467  starE1 = (starXX - starYY)/(starXX + starYY)
468  starE2 = 2*starXY/(starXX + starYY)
469  medianSize = np.median(starSize)
470 
471  psfSize = np.power(psfXX*psfYY - psfXY**2, 0.25)
472  psfE1 = (psfXX - psfYY)/(psfXX + psfYY)
473  psfE2 = 2*psfXY/(psfXX + psfYY)
474 
475  medianE1 = np.abs(np.median(starE1 - psfE1))
476  medianE2 = np.abs(np.median(starE2 - psfE2))
477  medianE = np.sqrt(medianE1**2 + medianE2**2)
478 
479  scatterSize = sigmaMad(starSize - psfSize)
480  scaledScatterSize = scatterSize/medianSize**2
481 
482  valid = True
483  if self.config.maxEllipResidual and medianE > self.config.maxEllipResidual:
484  self.log.info("Removing visit %s because median e residual too large: %f vs %f",
485  dataId, medianE, self.config.maxEllipResidual)
486  valid = False
487  elif self.config.maxSizeScatter and scatterSize > self.config.maxSizeScatter:
488  self.log.info("Removing visit %s because size scatter is too large: %f vs %f",
489  dataId, scatterSize, self.config.maxSizeScatter)
490  valid = False
491  elif self.config.maxScaledSizeScatter and scaledScatterSize > self.config.maxScaledSizeScatter:
492  self.log.info("Removing visit %s because scaled size scatter is too large: %f vs %f",
493  dataId, scaledScatterSize, self.config.maxScaledSizeScatter)
494  valid = False
495 
496  return valid
497 
498 
499 class BestSeeingWcsSelectImageConfig(WcsSelectImagesTask.ConfigClass):
500  """Base configuration for BestSeeingSelectImagesTask.
501  """
502  nImagesMax = pexConfig.RangeField(
503  dtype=int,
504  doc="Maximum number of images to select",
505  default=5,
506  min=0)
507  maxPsfFwhm = pexConfig.Field(
508  dtype=float,
509  doc="Maximum PSF FWHM (in arcseconds) to select",
510  default=1.5,
511  optional=True)
512  minPsfFwhm = pexConfig.Field(
513  dtype=float,
514  doc="Minimum PSF FWHM (in arcseconds) to select",
515  default=0.,
516  optional=True)
517 
518 
519 class BestSeeingWcsSelectImagesTask(WcsSelectImagesTask):
520  """Select up to a maximum number of the best-seeing images using their Wcs.
521  """
522  ConfigClass = BestSeeingWcsSelectImageConfig
523 
524  def runDataRef(self, dataRef, coordList, makeDataRefList=True,
525  selectDataList=None):
526  """Select the best-seeing images in the selectDataList that overlap the patch.
527 
528  This method is the old entry point for the Gen2 commandline tasks and drivers
529  Will be deprecated in v22.
530 
531  Parameters
532  ----------
533  dataRef : `lsst.daf.persistence.ButlerDataRef`
534  Data reference for coadd/tempExp (with tract, patch)
535  coordList : `list` of `lsst.geom.SpherePoint`
536  List of ICRS sky coordinates specifying boundary of patch
537  makeDataRefList : `boolean`, optional
538  Construct a list of data references?
539  selectDataList : `list` of `SelectStruct`
540  List of SelectStruct, to consider for selection
541 
542  Returns
543  -------
544  result : `lsst.pipe.base.Struct`
545  Result struct with components:
546  - ``exposureList``: the selected exposures
547  (`list` of `lsst.pipe.tasks.selectImages.BaseExposureInfo`).
548  - ``dataRefList``: the optional data references corresponding to
549  each element of ``exposureList``
550  (`list` of `lsst.daf.persistence.ButlerDataRef`, or `None`).
551  """
552  psfSizes = []
553  dataRefList = []
554  exposureInfoList = []
555 
556  if selectDataList is None:
557  selectDataList = []
558 
559  result = super().runDataRef(dataRef, coordList, makeDataRefList=True, selectDataList=selectDataList)
560 
561  for dataRef, exposureInfo in zip(result.dataRefList, result.exposureInfoList):
562  cal = dataRef.get("calexp", immediate=True)
563 
564  # if min/max PSF values are defined, remove images out of bounds
565  pixToArcseconds = cal.getWcs().getPixelScale().asArcseconds()
566  psfSize = cal.getPsf().computeShape().getDeterminantRadius()*pixToArcseconds
567  sizeFwhm = psfSize * np.sqrt(8.*np.log(2.))
568  if self.config.maxPsfFwhm and sizeFwhm > self.config.maxPsfFwhm:
569  continue
570  if self.config.minPsfFwhm and sizeFwhm < self.config.minPsfFwhm:
571  continue
572  psfSizes.append(sizeFwhm)
573  dataRefList.append(dataRef)
574  exposureInfoList.append(exposureInfo)
575 
576  if len(psfSizes) > self.config.nImagesMax:
577  sortedIndices = np.argsort(psfSizes)[:self.config.nImagesMax]
578  filteredDataRefList = [dataRefList[i] for i in sortedIndices]
579  filteredExposureInfoList = [exposureInfoList[i] for i in sortedIndices]
580  self.log.info("%d images selected with FWHM range of %f--%f arcseconds",
581  len(sortedIndices), psfSizes[sortedIndices[0]], psfSizes[sortedIndices[-1]])
582 
583  else:
584  if len(psfSizes) == 0:
585  self.log.warning("0 images selected.")
586  else:
587  self.log.debug("%d images selected with FWHM range of %d--%d arcseconds",
588  len(psfSizes), psfSizes[0], psfSizes[-1])
589  filteredDataRefList = dataRefList
590  filteredExposureInfoList = exposureInfoList
591 
592  return pipeBase.Struct(
593  dataRefList=filteredDataRefList if makeDataRefList else None,
594  exposureInfoList=filteredExposureInfoList,
595  )
596 
597 
598 class BestSeeingSelectVisitsConnections(pipeBase.PipelineTaskConnections,
599  dimensions=("tract", "patch", "skymap", "band", "instrument"),
600  defaultTemplates={"coaddName": "goodSeeing"}):
601  skyMap = pipeBase.connectionTypes.Input(
602  doc="Input definition of geometry/bbox and projection/wcs for coadded exposures",
603  name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
604  storageClass="SkyMap",
605  dimensions=("skymap",),
606  )
607  visitSummaries = pipeBase.connectionTypes.Input(
608  doc="Per-visit consolidated exposure metadata from ConsolidateVisitSummaryTask",
609  name="visitSummary",
610  storageClass="ExposureCatalog",
611  dimensions=("instrument", "visit",),
612  multiple=True,
613  deferLoad=True
614  )
615  goodVisits = pipeBase.connectionTypes.Output(
616  doc="Selected visits to be coadded.",
617  name="{coaddName}Visits",
618  storageClass="StructuredDataDict",
619  dimensions=("instrument", "tract", "patch", "skymap", "band"),
620  )
621 
622 
623 class BestSeeingSelectVisitsConfig(pipeBase.PipelineTaskConfig,
624  pipelineConnections=BestSeeingSelectVisitsConnections):
625  nVisitsMax = pexConfig.RangeField(
626  dtype=int,
627  doc="Maximum number of visits to select",
628  default=12,
629  min=0
630  )
631  maxPsfFwhm = pexConfig.Field(
632  dtype=float,
633  doc="Maximum PSF FWHM (in arcseconds) to select",
634  default=1.5,
635  optional=True
636  )
637  minPsfFwhm = pexConfig.Field(
638  dtype=float,
639  doc="Minimum PSF FWHM (in arcseconds) to select",
640  default=0.,
641  optional=True
642  )
643  doConfirmOverlap = pexConfig.Field(
644  dtype=bool,
645  doc="Do remove visits that do not actually overlap the patch?",
646  default=True,
647  )
648  minMJD = pexConfig.Field(
649  dtype=float,
650  doc="Minimum visit MJD to select",
651  default=None,
652  optional=True
653  )
654  maxMJD = pexConfig.Field(
655  dtype=float,
656  doc="Maximum visit MJD to select",
657  default=None,
658  optional=True
659  )
660 
661 
662 class BestSeeingSelectVisitsTask(pipeBase.PipelineTask):
663  """Select up to a maximum number of the best-seeing visits
664 
665  Don't exceed the FWHM range specified by configs min(max)PsfFwhm.
666  This Task is a port of the Gen2 image-selector used in the AP pipeline:
667  BestSeeingSelectImagesTask. This Task selects full visits based on the
668  average PSF of the entire visit.
669  """
670  ConfigClass = BestSeeingSelectVisitsConfig
671  _DefaultName = 'bestSeeingSelectVisits'
672 
673  def runQuantum(self, butlerQC, inputRefs, outputRefs):
674  inputs = butlerQC.get(inputRefs)
675  quantumDataId = butlerQC.quantum.dataId
676  outputs = self.run(**inputs, dataId=quantumDataId)
677  butlerQC.put(outputs, outputRefs)
678 
679  def run(self, visitSummaries, skyMap, dataId):
680  """Run task
681 
682  Parameters:
683  -----------
684  visitSummary : `list`
685  List of `lsst.pipe.base.connections.DeferredDatasetRef` of
686  visitSummary tables of type `lsst.afw.table.ExposureCatalog`
687  skyMap : `lsst.skyMap.SkyMap`
688  SkyMap for checking visits overlap patch
689  dataId : `dict` of dataId keys
690  For retrieving patch info for checking visits overlap patch
691 
692  Returns
693  -------
694  result : `lsst.pipe.base.Struct`
695  Result struct with components:
696 
697  - `goodVisits`: `dict` with selected visit ids as keys,
698  so that it can be be saved as a StructuredDataDict.
699  StructuredDataList's are currently limited.
700  """
701 
702  if self.config.doConfirmOverlap:
703  patchPolygon = self.makePatchPolygon(skyMap, dataId)
704 
705  inputVisits = [visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries]
706  fwhmSizes = []
707  visits = []
708  for visit, visitSummary in zip(inputVisits, visitSummaries):
709  # read in one-by-one and only once. There may be hundreds
710  visitSummary = visitSummary.get()
711 
712  # mjd is guaranteed to be the same for every detector in the visitSummary.
713  mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD)
714 
715  pixToArcseconds = [vs.getWcs().getPixelScale(vs.getBBox().getCenter()).asArcseconds()
716  for vs in visitSummary]
717  # psfSigma is PSF model determinant radius at chip center in pixels
718  psfSigmas = np.array([vs['psfSigma'] for vs in visitSummary])
719  fwhm = np.nanmean(psfSigmas * pixToArcseconds) * np.sqrt(8.*np.log(2.))
720 
721  if self.config.maxPsfFwhm and fwhm > self.config.maxPsfFwhm:
722  continue
723  if self.config.minPsfFwhm and fwhm < self.config.minPsfFwhm:
724  continue
725  if self.config.minMJD and mjd < self.config.minMJD:
726  self.log.debug('MJD %f earlier than %.2f; rejecting', mjd, self.config.minMJD)
727  continue
728  if self.config.maxMJD and mjd > self.config.maxMJD:
729  self.log.debug('MJD %f later than %.2f; rejecting', mjd, self.config.maxMJD)
730  continue
731  if self.config.doConfirmOverlap and not self.doesIntersectPolygon(visitSummary, patchPolygon):
732  continue
733 
734  fwhmSizes.append(fwhm)
735  visits.append(visit)
736 
737  sortedVisits = [ind for (_, ind) in sorted(zip(fwhmSizes, visits))]
738  output = sortedVisits[:self.config.nVisitsMax]
739  self.log.info("%d images selected with FWHM range of %d--%d arcseconds",
740  len(output), fwhmSizes[visits.index(output[0])], fwhmSizes[visits.index(output[-1])])
741 
742  # In order to store as a StructuredDataDict, convert list to dict
743  goodVisits = {key: True for key in output}
744  return pipeBase.Struct(goodVisits=goodVisits)
745 
746  def makePatchPolygon(self, skyMap, dataId):
747  """Return True if sky polygon overlaps visit
748 
749  Parameters:
750  -----------
751  skyMap : `lsst.afw.table.ExposureCatalog`
752  Exposure catalog with per-detector geometry
753  dataId : `dict` of dataId keys
754  For retrieving patch info
755 
756  Returns:
757  --------
758  result :` lsst.sphgeom.ConvexPolygon.convexHull`
759  Polygon of patch's outer bbox
760  """
761  wcs = skyMap[dataId['tract']].getWcs()
762  bbox = skyMap[dataId['tract']][dataId['patch']].getOuterBBox()
763  sphCorners = wcs.pixelToSky(lsst.geom.Box2D(bbox).getCorners())
764  result = lsst.sphgeom.ConvexPolygon.convexHull([coord.getVector() for coord in sphCorners])
765  return result
766 
767  def doesIntersectPolygon(self, visitSummary, polygon):
768  """Return True if sky polygon overlaps visit
769 
770  Parameters:
771  -----------
772  visitSummary : `lsst.afw.table.ExposureCatalog`
773  Exposure catalog with per-detector geometry
774  polygon :` lsst.sphgeom.ConvexPolygon.convexHull`
775  Polygon to check overlap
776 
777  Returns:
778  --------
779  doesIntersect: `bool`
780  Does the visit overlap the polygon
781  """
782  doesIntersect = False
783  for detectorSummary in visitSummary:
784  corners = [lsst.geom.SpherePoint(ra, decl, units=lsst.geom.degrees).getVector() for (ra, decl) in
785  zip(detectorSummary['raCorners'], detectorSummary['decCorners'])]
786  detectorPolygon = lsst.sphgeom.ConvexPolygon.convexHull(corners)
787  if detectorPolygon.intersects(polygon):
788  doesIntersect = True
789  break
790  return doesIntersect
791 
792 
793 class BestSeeingQuantileSelectVisitsConfig(pipeBase.PipelineTaskConfig,
794  pipelineConnections=BestSeeingSelectVisitsConnections):
795  qMin = pexConfig.RangeField(
796  doc="Lower bound of quantile range to select. Sorts visits by seeing from narrow to wide, "
797  "and select those in the interquantile range (qMin, qMax). Set qMin to 0 for Best Seeing. "
798  "This config should be changed from zero only for exploratory diffIm testing.",
799  dtype=float,
800  default=0,
801  min=0,
802  max=1,
803  )
804  qMax = pexConfig.RangeField(
805  doc="Upper bound of quantile range to select. Sorts visits by seeing from narrow to wide, "
806  "and select those in the interquantile range (qMin, qMax). Set qMax to 1 for Worst Seeing.",
807  dtype=float,
808  default=0.33,
809  min=0,
810  max=1,
811  )
812  nVisitsMin = pexConfig.Field(
813  doc="At least this number of visits selected and supercedes quantile. For example, if 10 visits "
814  "cover this patch, qMin=0.33, and nVisitsMin=5, the best 5 visits will be selected.",
815  dtype=int,
816  default=6,
817  )
818  doConfirmOverlap = pexConfig.Field(
819  dtype=bool,
820  doc="Do remove visits that do not actually overlap the patch?",
821  default=True,
822  )
823  minMJD = pexConfig.Field(
824  dtype=float,
825  doc="Minimum visit MJD to select",
826  default=None,
827  optional=True
828  )
829  maxMJD = pexConfig.Field(
830  dtype=float,
831  doc="Maximum visit MJD to select",
832  default=None,
833  optional=True
834  )
835 
836 
837 class BestSeeingQuantileSelectVisitsTask(BestSeeingSelectVisitsTask):
838  """Select a quantile of the best-seeing visits
839 
840  Selects the best (for example, third) full visits based on the average
841  PSF width in the entire visit. It can also be used for difference imaging
842  experiments that require templates with the worst seeing visits.
843  For example, selecting the worst third can be acheived by
844  changing the config parameters qMin to 0.66 and qMax to 1.
845  """
846  ConfigClass = BestSeeingQuantileSelectVisitsConfig
847  _DefaultName = 'bestSeeingQuantileSelectVisits'
848 
849  @utils.inheritDoc(BestSeeingSelectVisitsTask)
850  def run(self, visitSummaries, skyMap, dataId):
851  if self.config.doConfirmOverlap:
852  patchPolygon = self.makePatchPolygon(skyMap, dataId)
853  visits = np.array([visitSummary.ref.dataId['visit'] for visitSummary in visitSummaries])
854  radius = np.empty(len(visits))
855  intersects = np.full(len(visits), True)
856  for i, visitSummary in enumerate(visitSummaries):
857  # read in one-by-one and only once. There may be hundreds
858  visitSummary = visitSummary.get()
859  # psfSigma is PSF model determinant radius at chip center in pixels
860  psfSigma = np.nanmedian([vs['psfSigma'] for vs in visitSummary])
861  radius[i] = psfSigma
862  if self.config.doConfirmOverlap:
863  intersects[i] = self.doesIntersectPolygon(visitSummary, patchPolygon)
864  if self.config.minMJD or self.config.maxMJD:
865  # mjd is guaranteed to be the same for every detector in the visitSummary.
866  mjd = visitSummary[0].getVisitInfo().getDate().get(system=DateTime.MJD)
867  aboveMin = mjd > self.config.minMJD if self.config.minMJD else True
868  belowMax = mjd < self.config.maxMJD if self.config.maxMJD else True
869  intersects[i] = intersects[i] and aboveMin and belowMax
870 
871  sortedVisits = [v for rad, v in sorted(zip(radius[intersects], visits[intersects]))]
872  lowerBound = min(int(np.round(self.config.qMin*len(visits[intersects]))),
873  max(0, len(visits[intersects]) - self.config.nVisitsMin))
874  upperBound = max(int(np.round(self.config.qMax*len(visits[intersects]))), self.config.nVisitsMin)
875 
876  # In order to store as a StructuredDataDict, convert list to dict
877  goodVisits = {int(visit): True for visit in sortedVisits[lowerBound:upperBound]}
878  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
def __init__(self, dataId, coordList)
Definition: selectImages.py:63
def runDataRef(self, dataRef, coordList, makeDataRefList=True, selectDataList=[])
def __init__(self, dataRef, wcs, bbox)
def getValidImageCorners(self, imageWcs, imageBox, patchPoly, dataId=None)
def runDataRef(self, dataRef, coordList, makeDataRefList=True, selectDataList=[])
def run(self, wcsList, bboxList, coordList, dataIds=None, **kwargs)
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...
Definition: ConvexPolygon.h:65
daf::base::PropertyList * list
Definition: fits.cc:913
daf::base::PropertySet * set
Definition: fits.cc:912
bool isValid
Definition: fits.cc:399
def run(self, coaddExposures, bbox, wcs)
Definition: getTemplate.py:603