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
matchBackgrounds.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 # MERCHANTABILIY 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 import numpy
22 import lsst.afw.image as afwImage
23 import lsst.afw.math as afwMath
24 import lsst.geom as geom
25 import lsst.pex.config as pexConfig
26 import lsst.pipe.base as pipeBase
27 import lsstDebug
28 from lsst.utils.timer import timeMethod
29 
30 
31 class MatchBackgroundsConfig(pexConfig.Config):
32 
33  usePolynomial = pexConfig.Field(
34  dtype=bool,
35  doc="Fit background difference with Chebychev polynomial interpolation "
36  "(using afw.math.Approximate)? If False, fit with spline interpolation using afw.math.Background",
37  default=False
38  )
39  order = pexConfig.Field(
40  dtype=int,
41  doc="Order of Chebyshev polynomial background model. Ignored if usePolynomial False",
42  default=8
43  )
44  badMaskPlanes = pexConfig.ListField(
45  doc="Names of mask planes to ignore while estimating the background",
46  dtype=str, default=["NO_DATA", "DETECTED", "DETECTED_NEGATIVE", "SAT", "BAD", "INTRP", "CR"],
47  itemCheck=lambda x: x in afwImage.Mask().getMaskPlaneDict(),
48  )
49  gridStatistic = pexConfig.ChoiceField(
50  dtype=str,
51  doc="Type of statistic to estimate pixel value for the grid points",
52  default="MEAN",
53  allowed={
54  "MEAN": "mean",
55  "MEDIAN": "median",
56  "MEANCLIP": "clipped mean"
57  }
58  )
59  undersampleStyle = pexConfig.ChoiceField(
60  doc="Behaviour if there are too few points in grid for requested interpolation style. "
61  "Note: INCREASE_NXNYSAMPLE only allowed for usePolynomial=True.",
62  dtype=str,
63  default="REDUCE_INTERP_ORDER",
64  allowed={
65  "THROW_EXCEPTION": "throw an exception if there are too few points",
66  "REDUCE_INTERP_ORDER": "use an interpolation style with a lower order.",
67  "INCREASE_NXNYSAMPLE": "Increase the number of samples used to make the interpolation grid.",
68  }
69  )
70  binSize = pexConfig.Field(
71  doc="Bin size for gridding the difference image and fitting a spatial model",
72  dtype=int,
73  default=256
74  )
75  interpStyle = pexConfig.ChoiceField(
76  dtype=str,
77  doc="Algorithm to interpolate the background values; ignored if usePolynomial is True"
78  "Maps to an enum; see afw.math.Background",
79  default="AKIMA_SPLINE",
80  allowed={
81  "CONSTANT": "Use a single constant value",
82  "LINEAR": "Use linear interpolation",
83  "NATURAL_SPLINE": "cubic spline with zero second derivative at endpoints",
84  "AKIMA_SPLINE": "higher-level nonlinear spline that is more robust to outliers",
85  "NONE": "No background estimation is to be attempted",
86  }
87  )
88  numSigmaClip = pexConfig.Field(
89  dtype=int,
90  doc="Sigma for outlier rejection; ignored if gridStatistic != 'MEANCLIP'.",
91  default=3
92  )
93  numIter = pexConfig.Field(
94  dtype=int,
95  doc="Number of iterations of outlier rejection; ignored if gridStatistic != 'MEANCLIP'.",
96  default=2
97  )
98  bestRefWeightCoverage = pexConfig.RangeField(
99  dtype=float,
100  doc="Weight given to coverage (number of pixels that overlap with patch), "
101  "when calculating best reference exposure. Higher weight prefers exposures with high coverage."
102  "Ignored when reference visit is supplied",
103  default=0.4,
104  min=0., max=1.
105  )
106  bestRefWeightVariance = pexConfig.RangeField(
107  dtype=float,
108  doc="Weight given to image variance when calculating best reference exposure. "
109  "Higher weight prefers exposures with low image variance. Ignored when reference visit is supplied",
110  default=0.4,
111  min=0., max=1.
112  )
113  bestRefWeightLevel = pexConfig.RangeField(
114  dtype=float,
115  doc="Weight given to mean background level when calculating best reference exposure. "
116  "Higher weight prefers exposures with low mean background level. "
117  "Ignored when reference visit is supplied.",
118  default=0.2,
119  min=0., max=1.
120  )
121  approxWeighting = pexConfig.Field(
122  dtype=bool,
123  doc=("Use inverse-variance weighting when approximating background offset model? "
124  "This will fail when the background offset is constant "
125  "(this is usually only the case in testing with artificial images)."
126  "(usePolynomial=True)"),
127  default=True,
128  )
129  gridStdevEpsilon = pexConfig.RangeField(
130  dtype=float,
131  doc="Tolerance on almost zero standard deviation in a background-offset grid bin. "
132  "If all bins have a standard deviation below this value, the background offset model "
133  "is approximated without inverse-variance weighting. (usePolynomial=True)",
134  default=1e-8,
135  min=0.
136  )
137 
138 
139 class MatchBackgroundsTask(pipeBase.Task):
140  ConfigClass = MatchBackgroundsConfig
141  _DefaultName = "matchBackgrounds"
142 
143  def __init__(self, *args, **kwargs):
144  pipeBase.Task.__init__(self, *args, **kwargs)
145 
147  self.sctrlsctrl.setAndMask(afwImage.Mask.getPlaneBitMask(self.config.badMaskPlanes))
148  self.sctrlsctrl.setNanSafe(True)
149 
150  @timeMethod
151  def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=None, refImageScaler=None):
152  """Match the backgrounds of a list of coadd temp exposures to a reference coadd temp exposure.
153 
154  Choose a refExpDataRef automatically if none supplied.
155 
156  @param[in] expRefList: list of data references to science exposures to be background-matched;
157  all exposures must exist.
158  @param[in] expDatasetType: dataset type of exposures, e.g. 'goodSeeingCoadd_tempExp'
159  @param[in] imageScalerList: list of image scalers (coaddUtils.ImageScaler);
160  if None then the images are not scaled
161  @param[in] refExpDataRef: data reference for the reference exposure.
162  If None, then this task selects the best exposures from expRefList.
163  if not None then must be one of the exposures in expRefList.
164  @param[in] refImageScaler: image scaler for reference image;
165  ignored if refExpDataRef is None, else scaling is not performed if None
166 
167  @return: a pipBase.Struct containing these fields:
168  - backgroundInfoList: a list of pipeBase.Struct, one per exposure in expRefList,
169  each of which contains these fields:
170  - isReference: this is the reference exposure (only one returned Struct will
171  contain True for this value, unless the ref exposure is listed multiple times)
172  - backgroundModel: differential background model (afw.Math.Background or afw.Math.Approximate).
173  Add this to the science exposure to match the reference exposure.
174  - fitRMS: rms of the fit. This is the sqrt(mean(residuals**2)).
175  - matchedMSE: the MSE of the reference and matched images: mean((refImage - matchedSciImage)**2);
176  should be comparable to difference image's mean variance.
177  - diffImVar: the mean variance of the difference image.
178  All fields except isReference will be None if isReference True or the fit failed.
179 
180  @warning: all exposures must exist on disk
181  """
182 
183  numExp = len(expRefList)
184  if numExp < 1:
185  raise pipeBase.TaskError("No exposures to match")
186 
187  if expDatasetType is None:
188  raise pipeBase.TaskError("Must specify expDatasetType")
189 
190  if imageScalerList is None:
191  self.log.info("imageScalerList is None; no scaling will be performed")
192  imageScalerList = [None] * numExp
193 
194  if len(expRefList) != len(imageScalerList):
195  raise RuntimeError("len(expRefList) = %s != %s = len(imageScalerList)" %
196  (len(expRefList), len(imageScalerList)))
197 
198  refInd = None
199  if refExpDataRef is None:
200  # select the best reference exposure from expRefList
201  refInd = self.selectRefExposureselectRefExposure(
202  expRefList=expRefList,
203  imageScalerList=imageScalerList,
204  expDatasetType=expDatasetType,
205  )
206  refExpDataRef = expRefList[refInd]
207  refImageScaler = imageScalerList[refInd]
208 
209  # refIndSet is the index of all exposures in expDataList that match the reference.
210  # It is used to avoid background-matching an exposure to itself. It is a list
211  # because it is possible (though unlikely) that expDataList will contain duplicates.
212  expKeyList = refExpDataRef.butlerSubset.butler.getKeys(expDatasetType)
213  refMatcher = DataRefMatcher(refExpDataRef.butlerSubset.butler, expDatasetType)
214  refIndSet = set(refMatcher.matchList(ref0=refExpDataRef, refList=expRefList))
215 
216  if refInd is not None and refInd not in refIndSet:
217  raise RuntimeError("Internal error: selected reference %s not found in expRefList")
218 
219  refExposure = refExpDataRef.get(expDatasetType, immediate=True)
220  if refImageScaler is not None:
221  refMI = refExposure.getMaskedImage()
222  refImageScaler.scaleMaskedImage(refMI)
223 
224  debugIdKeyList = tuple(set(expKeyList) - set(['tract', 'patch']))
225 
226  self.log.info("Matching %d Exposures", numExp)
227 
228  backgroundInfoList = []
229  for ind, (toMatchRef, imageScaler) in enumerate(zip(expRefList, imageScalerList)):
230  if ind in refIndSet:
231  backgroundInfoStruct = pipeBase.Struct(
232  isReference=True,
233  backgroundModel=None,
234  fitRMS=0.0,
235  matchedMSE=None,
236  diffImVar=None,
237  )
238  else:
239  self.log.info("Matching background of %s to %s", toMatchRef.dataId, refExpDataRef.dataId)
240  try:
241  toMatchExposure = toMatchRef.get(expDatasetType, immediate=True)
242  if imageScaler is not None:
243  toMatchMI = toMatchExposure.getMaskedImage()
244  imageScaler.scaleMaskedImage(toMatchMI)
245  # store a string specifying the visit to label debug plot
246  self.debugDataIdStringdebugDataIdString = ''.join([str(toMatchRef.dataId[vk]) for vk in debugIdKeyList])
247  backgroundInfoStruct = self.matchBackgroundsmatchBackgrounds(
248  refExposure=refExposure,
249  sciExposure=toMatchExposure,
250  )
251  backgroundInfoStruct.isReference = False
252  except Exception as e:
253  self.log.warning("Failed to fit background %s: %s", toMatchRef.dataId, e)
254  backgroundInfoStruct = pipeBase.Struct(
255  isReference=False,
256  backgroundModel=None,
257  fitRMS=None,
258  matchedMSE=None,
259  diffImVar=None,
260  )
261 
262  backgroundInfoList.append(backgroundInfoStruct)
263 
264  return pipeBase.Struct(
265  backgroundInfoList=backgroundInfoList)
266 
267  @timeMethod
268  def selectRefExposure(self, expRefList, imageScalerList, expDatasetType):
269  """Find best exposure to use as the reference exposure
270 
271  Calculate an appropriate reference exposure by minimizing a cost function that penalizes
272  high variance, high background level, and low coverage. Use the following config parameters:
273  - bestRefWeightCoverage
274  - bestRefWeightVariance
275  - bestRefWeightLevel
276 
277  @param[in] expRefList: list of data references to exposures.
278  Retrieves dataset type specified by expDatasetType.
279  If an exposure is not found, it is skipped with a warning.
280  @param[in] imageScalerList: list of image scalers (coaddUtils.ImageScaler);
281  must be the same length as expRefList
282  @param[in] expDatasetType: dataset type of exposure: e.g. 'goodSeeingCoadd_tempExp'
283 
284  @return: index of best exposure
285 
286  @raise pipeBase.TaskError if none of the exposures in expRefList are found.
287  """
288  self.log.info("Calculating best reference visit")
289  varList = []
290  meanBkgdLevelList = []
291  coverageList = []
292 
293  if len(expRefList) != len(imageScalerList):
294  raise RuntimeError("len(expRefList) = %s != %s = len(imageScalerList)" %
295  (len(expRefList), len(imageScalerList)))
296 
297  for expRef, imageScaler in zip(expRefList, imageScalerList):
298  exposure = expRef.get(expDatasetType, immediate=True)
299  maskedImage = exposure.getMaskedImage()
300  if imageScaler is not None:
301  try:
302  imageScaler.scaleMaskedImage(maskedImage)
303  except Exception:
304  # need to put a place holder in Arr
305  varList.append(numpy.nan)
306  meanBkgdLevelList.append(numpy.nan)
307  coverageList.append(numpy.nan)
308  continue
309  statObjIm = afwMath.makeStatistics(maskedImage.getImage(), maskedImage.getMask(),
310  afwMath.MEAN | afwMath.NPOINT | afwMath.VARIANCE, self.sctrlsctrl)
311  meanVar, meanVarErr = statObjIm.getResult(afwMath.VARIANCE)
312  meanBkgdLevel, meanBkgdLevelErr = statObjIm.getResult(afwMath.MEAN)
313  npoints, npointsErr = statObjIm.getResult(afwMath.NPOINT)
314  varList.append(meanVar)
315  meanBkgdLevelList.append(meanBkgdLevel)
316  coverageList.append(npoints)
317  if not coverageList:
318  raise pipeBase.TaskError(
319  "None of the candidate %s exist; cannot select best reference exposure" % (expDatasetType,))
320 
321  # Normalize metrics to range from 0 to 1
322  varArr = numpy.array(varList)/numpy.nanmax(varList)
323  meanBkgdLevelArr = numpy.array(meanBkgdLevelList)/numpy.nanmax(meanBkgdLevelList)
324  coverageArr = numpy.nanmin(coverageList)/numpy.array(coverageList)
325 
326  costFunctionArr = self.config.bestRefWeightVariance * varArr
327  costFunctionArr += self.config.bestRefWeightLevel * meanBkgdLevelArr
328  costFunctionArr += self.config.bestRefWeightCoverage * coverageArr
329  return numpy.nanargmin(costFunctionArr)
330 
331  @timeMethod
332  def matchBackgrounds(self, refExposure, sciExposure):
333  """
334  Match science exposure's background level to that of reference exposure.
335 
336  Process creates a difference image of the reference exposure minus the science exposure, and then
337  generates an afw.math.Background object. It assumes (but does not require/check) that the mask plane
338  already has detections set. If detections have not been set/masked, sources will bias the
339  background estimation.
340  The 'background' of the difference image is smoothed by spline interpolation (by the Background class)
341  or by polynomial interpolation by the Approximate class. This model of difference image
342  is added to the science exposure in memory.
343  Fit diagnostics are also calculated and returned.
344 
345  @param[in] refExposure: reference exposure
346  @param[in,out] sciExposure: science exposure; modified by changing the background level
347  to match that of the reference exposure
348  @returns a pipBase.Struct with fields:
349  - backgroundModel: an afw.math.Approximate or an afw.math.Background.
350  - fitRMS: rms of the fit. This is the sqrt(mean(residuals**2)).
351  - matchedMSE: the MSE of the reference and matched images: mean((refImage - matchedSciImage)**2);
352  should be comparable to difference image's mean variance.
353  - diffImVar: the mean variance of the difference image.
354  """
355 
356  if lsstDebug.Info(__name__).savefits:
357  refExposure.writeFits(lsstDebug.Info(__name__).figpath + 'refExposure.fits')
358  sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciExposure.fits')
359 
360  # Check Configs for polynomials:
361  if self.config.usePolynomial:
362  x, y = sciExposure.getDimensions()
363  shortSideLength = min(x, y)
364  if shortSideLength < self.config.binSize:
365  raise ValueError("%d = config.binSize > shorter dimension = %d" % (self.config.binSize,
366  shortSideLength))
367  npoints = shortSideLength // self.config.binSize
368  if shortSideLength % self.config.binSize != 0:
369  npoints += 1
370 
371  if self.config.order > npoints - 1:
372  raise ValueError("%d = config.order > npoints - 1 = %d" % (self.config.order, npoints - 1))
373 
374  # Check that exposures are same shape
375  if (sciExposure.getDimensions() != refExposure.getDimensions()):
376  wSci, hSci = sciExposure.getDimensions()
377  wRef, hRef = refExposure.getDimensions()
378  raise RuntimeError(
379  "Exposures are different dimensions. sci:(%i, %i) vs. ref:(%i, %i)" %
380  (wSci, hSci, wRef, hRef))
381 
382  statsFlag = getattr(afwMath, self.config.gridStatistic)
383  self.sctrlsctrl.setNumSigmaClip(self.config.numSigmaClip)
384  self.sctrlsctrl.setNumIter(self.config.numIter)
385 
386  im = refExposure.getMaskedImage()
387  diffMI = im.Factory(im, True)
388  diffMI -= sciExposure.getMaskedImage()
389 
390  width = diffMI.getWidth()
391  height = diffMI.getHeight()
392  nx = width // self.config.binSize
393  if width % self.config.binSize != 0:
394  nx += 1
395  ny = height // self.config.binSize
396  if height % self.config.binSize != 0:
397  ny += 1
398 
399  bctrl = afwMath.BackgroundControl(nx, ny, self.sctrlsctrl, statsFlag)
400  bctrl.setUndersampleStyle(self.config.undersampleStyle)
401 
402  bkgd = afwMath.makeBackground(diffMI, bctrl)
403 
404  # Some config and input checks if config.usePolynomial:
405  # 1) Check that order/bin size make sense:
406  # 2) Change binsize or order if underconstrained.
407  if self.config.usePolynomial:
408  order = self.config.order
409  bgX, bgY, bgZ, bgdZ = self._gridImage_gridImage(diffMI, self.config.binSize, statsFlag)
410  minNumberGridPoints = min(len(set(bgX)), len(set(bgY)))
411  if len(bgZ) == 0:
412  raise ValueError("No overlap with reference. Nothing to match")
413  elif minNumberGridPoints <= self.config.order:
414  # must either lower order or raise number of bins or throw exception
415  if self.config.undersampleStyle == "THROW_EXCEPTION":
416  raise ValueError("Image does not cover enough of ref image for order and binsize")
417  elif self.config.undersampleStyle == "REDUCE_INTERP_ORDER":
418  self.log.warning("Reducing order to %d", (minNumberGridPoints - 1))
419  order = minNumberGridPoints - 1
420  elif self.config.undersampleStyle == "INCREASE_NXNYSAMPLE":
421  newBinSize = (minNumberGridPoints*self.config.binSize) // (self.config.order + 1)
422  bctrl.setNxSample(newBinSize)
423  bctrl.setNySample(newBinSize)
424  bkgd = afwMath.makeBackground(diffMI, bctrl) # do over
425  self.log.warning("Decreasing binsize to %d", newBinSize)
426 
427  # If there is no variance in any image pixels, do not weight bins by inverse variance
428  isUniformImageDiff = not numpy.any(bgdZ > self.config.gridStdevEpsilon)
429  weightByInverseVariance = False if isUniformImageDiff else self.config.approxWeighting
430 
431  # Add offset to sciExposure
432  try:
433  if self.config.usePolynomial:
434  actrl = afwMath.ApproximateControl(afwMath.ApproximateControl.CHEBYSHEV,
435  order, order, weightByInverseVariance)
436  undersampleStyle = getattr(afwMath, self.config.undersampleStyle)
437  approx = bkgd.getApproximate(actrl, undersampleStyle)
438  bkgdImage = approx.getImage()
439  else:
440  bkgdImage = bkgd.getImageF(self.config.interpStyle, self.config.undersampleStyle)
441  except Exception as e:
442  raise RuntimeError("Background/Approximation failed to interp image %s: %s" % (
443  self.debugDataIdStringdebugDataIdString, e))
444 
445  sciMI = sciExposure.getMaskedImage()
446  sciMI += bkgdImage
447  del sciMI
448 
449  # Need RMS from fit: 2895 will replace this:
450  rms = 0.0
451  X, Y, Z, dZ = self._gridImage_gridImage(diffMI, self.config.binSize, statsFlag)
452  x0, y0 = diffMI.getXY0()
453  modelValueArr = numpy.empty(len(Z))
454  for i in range(len(X)):
455  modelValueArr[i] = bkgdImage[int(X[i]-x0), int(Y[i]-y0), afwImage.LOCAL]
456  resids = Z - modelValueArr
457  rms = numpy.sqrt(numpy.mean(resids[~numpy.isnan(resids)]**2))
458 
459  if lsstDebug.Info(__name__).savefits:
460  sciExposure.writeFits(lsstDebug.Info(__name__).figpath + 'sciMatchedExposure.fits')
461 
462  if lsstDebug.Info(__name__).savefig:
463  bbox = geom.Box2D(refExposure.getMaskedImage().getBBox())
464  try:
465  self._debugPlot_debugPlot(X, Y, Z, dZ, bkgdImage, bbox, modelValueArr, resids)
466  except Exception as e:
467  self.log.warning('Debug plot not generated: %s', e)
468 
469  meanVar = afwMath.makeStatistics(diffMI.getVariance(), diffMI.getMask(),
470  afwMath.MEANCLIP, self.sctrlsctrl).getValue()
471 
472  diffIm = diffMI.getImage()
473  diffIm -= bkgdImage # diffMI should now have a mean ~ 0
474  del diffIm
475  mse = afwMath.makeStatistics(diffMI, afwMath.MEANSQUARE, self.sctrlsctrl).getValue()
476 
477  outBkgd = approx if self.config.usePolynomial else bkgd
478 
479  return pipeBase.Struct(
480  backgroundModel=outBkgd,
481  fitRMS=rms,
482  matchedMSE=mse,
483  diffImVar=meanVar)
484 
485  def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids):
486  """Generate a plot showing the background fit and residuals.
487 
488  It is called when lsstDebug.Info(__name__).savefig = True
489  Saves the fig to lsstDebug.Info(__name__).figpath
490  Displays on screen if lsstDebug.Info(__name__).display = True
491 
492  @param X: array of x positions
493  @param Y: array of y positions
494  @param Z: array of the grid values that were interpolated
495  @param dZ: array of the error on the grid values
496  @param modelImage: image ofthe model of the fit
497  @param model: array of len(Z) containing the grid values predicted by the model
498  @param resids: Z - model
499  """
500  import matplotlib.pyplot as plt
501  import matplotlib.colors
502  from mpl_toolkits.axes_grid1 import ImageGrid
503  zeroIm = afwImage.MaskedImageF(geom.Box2I(bbox))
504  zeroIm += modelImage
505  x0, y0 = zeroIm.getXY0()
506  dx, dy = zeroIm.getDimensions()
507  if len(Z) == 0:
508  self.log.warning("No grid. Skipping plot generation.")
509  else:
510  max, min = numpy.max(Z), numpy.min(Z)
511  norm = matplotlib.colors.normalize(vmax=max, vmin=min)
512  maxdiff = numpy.max(numpy.abs(resids))
513  diffnorm = matplotlib.colors.normalize(vmax=maxdiff, vmin=-maxdiff)
514  rms = numpy.sqrt(numpy.mean(resids**2))
515  fig = plt.figure(1, (8, 6))
516  meanDz = numpy.mean(dZ)
517  grid = ImageGrid(fig, 111, nrows_ncols=(1, 2), axes_pad=0.1,
518  share_all=True, label_mode="L", cbar_mode="each",
519  cbar_size="7%", cbar_pad="2%", cbar_location="top")
520  im = grid[0].imshow(zeroIm.getImage().getArray(),
521  extent=(x0, x0+dx, y0+dy, y0), norm=norm,
522  cmap='Spectral')
523  im = grid[0].scatter(X, Y, c=Z, s=15.*meanDz/dZ, edgecolor='none', norm=norm,
524  marker='o', cmap='Spectral')
525  im2 = grid[1].scatter(X, Y, c=resids, edgecolor='none', norm=diffnorm,
526  marker='s', cmap='seismic')
527  grid.cbar_axes[0].colorbar(im)
528  grid.cbar_axes[1].colorbar(im2)
529  grid[0].axis([x0, x0+dx, y0+dy, y0])
530  grid[1].axis([x0, x0+dx, y0+dy, y0])
531  grid[0].set_xlabel("model and grid")
532  grid[1].set_xlabel("residuals. rms = %0.3f"%(rms))
533  if lsstDebug.Info(__name__).savefig:
534  fig.savefig(lsstDebug.Info(__name__).figpath + self.debugDataIdStringdebugDataIdString + '.png')
535  if lsstDebug.Info(__name__).display:
536  plt.show()
537  plt.clf()
538 
539  def _gridImage(self, maskedImage, binsize, statsFlag):
540  """Private method to grid an image for debugging"""
541  width, height = maskedImage.getDimensions()
542  x0, y0 = maskedImage.getXY0()
543  xedges = numpy.arange(0, width, binsize)
544  yedges = numpy.arange(0, height, binsize)
545  xedges = numpy.hstack((xedges, width)) # add final edge
546  yedges = numpy.hstack((yedges, height)) # add final edge
547 
548  # Use lists/append to protect against the case where
549  # a bin has no valid pixels and should not be included in the fit
550  bgX = []
551  bgY = []
552  bgZ = []
553  bgdZ = []
554 
555  for ymin, ymax in zip(yedges[0:-1], yedges[1:]):
556  for xmin, xmax in zip(xedges[0:-1], xedges[1:]):
557  subBBox = geom.Box2I(geom.PointI(int(x0 + xmin), int(y0 + ymin)),
558  geom.PointI(int(x0 + xmax-1), int(y0 + ymax-1)))
559  subIm = afwImage.MaskedImageF(maskedImage, subBBox, afwImage.PARENT, False)
560  stats = afwMath.makeStatistics(subIm,
561  afwMath.MEAN | afwMath.MEANCLIP | afwMath.MEDIAN
562  | afwMath.NPOINT | afwMath.STDEV,
563  self.sctrlsctrl)
564  npoints, _ = stats.getResult(afwMath.NPOINT)
565  if npoints >= 2:
566  stdev, _ = stats.getResult(afwMath.STDEV)
567  if stdev < self.config.gridStdevEpsilon:
568  stdev = self.config.gridStdevEpsilon
569  bgX.append(0.5 * (x0 + xmin + x0 + xmax))
570  bgY.append(0.5 * (y0 + ymin + y0 + ymax))
571  bgdZ.append(stdev/numpy.sqrt(npoints))
572  est, _ = stats.getResult(statsFlag)
573  bgZ.append(est)
574 
575  return numpy.array(bgX), numpy.array(bgY), numpy.array(bgZ), numpy.array(bgdZ)
576 
577 
579  """Match data references for a specified dataset type
580 
581  Note that this is not exact, but should suffice for this task
582  until there is better support for this kind of thing in the butler.
583  """
584 
585  def __init__(self, butler, datasetType):
586  """Construct a DataRefMatcher
587 
588  @param[in] butler
589  @param[in] datasetType: dataset type to match
590  """
591  self._datasetType_datasetType = datasetType # for diagnostics
592  self._keyNames_keyNames = butler.getKeys(datasetType)
593 
594  def _makeKey(self, ref):
595  """Return a tuple of values for the specified keyNames
596 
597  @param[in] ref: data reference
598 
599  @raise KeyError if ref.dataId is missing a key in keyNames
600  """
601  return tuple(ref.dataId[key] for key in self._keyNames_keyNames)
602 
603  def isMatch(self, ref0, ref1):
604  """Return True if ref0 == ref1
605 
606  @param[in] ref0: data ref 0
607  @param[in] ref1: data ref 1
608 
609  @raise KeyError if either ID is missing a key in keyNames
610  """
611  return self._makeKey_makeKey(ref0) == self._makeKey_makeKey(ref1)
612 
613  def matchList(self, ref0, refList):
614  """Return a list of indices of matches
615 
616  @return tuple of indices of matches
617 
618  @raise KeyError if any ID is missing a key in keyNames
619  """
620  key0 = self._makeKey_makeKey(ref0)
621  return tuple(ind for ind, ref in enumerate(refList) if self._makeKey_makeKey(ref) == key0)
int min
Represent a 2-dimensional array of bitmask pixels.
Definition: Mask.h:77
Control how to make an approximation.
Definition: Approximate.h:48
Pass parameters to a Background object.
Definition: Background.h:56
Pass parameters to a Statistics object.
Definition: Statistics.h:92
A floating-point coordinate rectangle geometry.
Definition: Box.h:413
An integer coordinate rectangle.
Definition: Box.h:55
def selectRefExposure(self, expRefList, imageScalerList, expDatasetType)
def run(self, expRefList, expDatasetType, imageScalerList=None, refExpDataRef=None, refImageScaler=None)
def _gridImage(self, maskedImage, binsize, statsFlag)
def _debugPlot(self, X, Y, Z, dZ, modelImage, bbox, model, resids)
def matchBackgrounds(self, refExposure, sciExposure)
daf::base::PropertySet * set
Definition: fits.cc:912
Backwards-compatibility support for depersisting the old Calib (FluxMag0/FluxMag0Err) objects.
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:359
std::shared_ptr< Background > makeBackground(ImageT const &img, BackgroundControl const &bgCtrl)
A convenience function that uses function overloading to make the correct type of Background.
Definition: Background.h:526