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