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