190 def run(self, mapperResults, exposure, **kwargs):
191 """Reduce a list of items produced by `ImageMapper`.
193 Either stitch the passed `mapperResults` list
194 together into a new Exposure (default) or pass it through
195 (if `self.config.reduceOperation` is 'none').
197 If `self.config.reduceOperation` is not 'none', then expect
198 that the `pipeBase.Struct`s in the `mapperResults` list
199 contain sub-exposures named 'subExposure', to be stitched back
200 into a single Exposure with the same dimensions, PSF, and mask
201 as the input `exposure`. Otherwise, the `mapperResults` list
202 is simply returned directly.
206 mapperResults : `list`
207 list of `lsst.pipe.base.Struct` returned by `ImageMapper.run`.
208 exposure : `lsst.afw.image.Exposure`
209 the original exposure which is cloned to use as the
210 basis for the resulting exposure (if
211 ``self.config.mapper.reduceOperation`` is not 'None')
213 additional keyword arguments propagated from
214 `ImageMapReduceTask.run`.
218 A `lsst.pipe.base.Struct` containing either an `lsst.afw.image.Exposure`
219 (named 'exposure') or a list (named 'result'),
220 depending on `config.reduceOperation`.
224 1. This currently correctly handles overlapping sub-exposures.
225 For overlapping sub-exposures, use `config.reduceOperation='average'`.
226 2. This correctly handles varying PSFs, constructing the resulting
227 exposure's PSF via CoaddPsf (DM-9629).
231 1. To be done: correct handling of masks (nearly there)
232 2. This logic currently makes *two* copies of the original exposure
233 (one here and one in `mapper.run()`). Possibly of concern
234 for large images on memory-constrained systems.
237 if self.config.reduceOperation ==
'none':
238 return pipeBase.Struct(result=mapperResults)
240 if self.config.reduceOperation ==
'coaddPsf':
242 coaddPsf = self._constructPsf(mapperResults, exposure)
243 return pipeBase.Struct(result=coaddPsf)
245 newExp = exposure.clone()
246 newMI = newExp.getMaskedImage()
248 reduceOp = self.config.reduceOperation
249 if reduceOp ==
'copy':
251 newMI.getImage()[:, :] = np.nan
252 newMI.getVariance()[:, :] = np.nan
254 newMI.getImage()[:, :] = 0.
255 newMI.getVariance()[:, :] = 0.
256 if reduceOp ==
'average':
257 weights = afwImage.ImageI(newMI.getBBox())
259 for item
in mapperResults:
260 item = item.subExposure
261 if not (isinstance(item, afwImage.ExposureF)
or isinstance(item, afwImage.ExposureI)
262 or isinstance(item, afwImage.ExposureU)
or isinstance(item, afwImage.ExposureD)):
263 raise TypeError(
"""Expecting an Exposure type, got %s.
264 Consider using `reduceOperation="none".""" % str(
type(item)))
265 subExp = newExp.Factory(newExp, item.getBBox())
266 subMI = subExp.getMaskedImage()
267 patchMI = item.getMaskedImage()
268 isValid = ~np.isnan(patchMI.getImage().getArray() * patchMI.getVariance().getArray())
270 if reduceOp ==
'copy':
271 subMI.getImage().getArray()[isValid] = patchMI.getImage().getArray()[isValid]
272 subMI.getVariance().getArray()[isValid] = patchMI.getVariance().getArray()[isValid]
273 subMI.getMask().getArray()[:, :] |= patchMI.getMask().getArray()
275 if reduceOp ==
'sum' or reduceOp ==
'average':
276 subMI.getImage().getArray()[isValid] += patchMI.getImage().getArray()[isValid]
277 subMI.getVariance().getArray()[isValid] += patchMI.getVariance().getArray()[isValid]
278 subMI.getMask().getArray()[:, :] |= patchMI.getMask().getArray()
279 if reduceOp ==
'average':
281 wtsView = afwImage.ImageI(weights, item.getBBox())
282 wtsView.getArray()[isValid] += 1
285 mask = newMI.getMask()
286 for m
in self.config.badMaskPlanes:
288 bad = mask.getPlaneBitMask(self.config.badMaskPlanes)
290 isNan = np.where(np.isnan(newMI.getImage().getArray() * newMI.getVariance().getArray()))
291 if len(isNan[0]) > 0:
293 mask.getArray()[isNan[0], isNan[1]] |= bad
295 if reduceOp ==
'average':
296 wts = weights.getArray().astype(np.float)
297 self.log.
info(
'AVERAGE: Maximum overlap: %f', np.nanmax(wts))
298 self.log.
info(
'AVERAGE: Average overlap: %f', np.nanmean(wts))
299 self.log.
info(
'AVERAGE: Minimum overlap: %f', np.nanmin(wts))
300 wtsZero = np.equal(wts, 0.)
301 wtsZeroInds = np.where(wtsZero)
302 wtsZeroSum = len(wtsZeroInds[0])
303 self.log.
info(
'AVERAGE: Number of zero pixels: %f (%f%%)', wtsZeroSum,
304 wtsZeroSum * 100. / wtsZero.size)
305 notWtsZero = ~wtsZero
306 tmp = newMI.getImage().getArray()
307 np.divide(tmp, wts, out=tmp, where=notWtsZero)
308 tmp = newMI.getVariance().getArray()
309 np.divide(tmp, wts, out=tmp, where=notWtsZero)
310 if len(wtsZeroInds[0]) > 0:
311 newMI.getImage().getArray()[wtsZeroInds] = np.nan
312 newMI.getVariance().getArray()[wtsZeroInds] = np.nan
315 mask.getArray()[wtsZeroInds] |= bad
318 if reduceOp ==
'sum' or reduceOp ==
'average':
319 psf = self._constructPsf(mapperResults, exposure)
322 return pipeBase.Struct(exposure=newExp)
def run(self, coaddExposures, bbox, wcs)