225 def run(self, exposure, catalog):
226 """Measure aperture correction
227
228 Parameters
229 ----------
230 exposure : `lsst.afw.image.Exposure`
231 Exposure aperture corrections are being measured on. The
232 bounding box is retrieved from it, and it is passed to the
233 sourceSelector. The output aperture correction map is *not*
234 added to the exposure; this is left to the caller.
235 catalog : `lsst.afw.table.SourceCatalog`
236 SourceCatalog containing measurements to be used to
237 compute aperture corrections.
238
239 Returns
240 -------
241 Struct : `lsst.pipe.base.Struct`
242 Contains the following:
243
244 ``apCorrMap``
245 aperture correction map (`lsst.afw.image.ApCorrMap`)
246 that contains two entries for each flux field:
247 - flux field (e.g. base_PsfFlux_instFlux): 2d model
248 - flux sigma field (e.g. base_PsfFlux_instFluxErr): 2d model of error
249 """
250 bbox = exposure.getBBox()
251 import lsstDebug
254
255 self.log.info("Measuring aperture corrections for %d flux fields", len(self.toCorrect))
256
257
258
259 selected = self.sourceSelector.run(catalog, exposure=exposure)
260
261 use = (
262 ~selected.sourceCat[self.refFluxNames.flagName]
263 & (np.isfinite(selected.sourceCat[self.refFluxNames.fluxName]))
264 )
265 goodRefCat = selected.sourceCat[use].copy()
266
267 apCorrMap = ApCorrMap()
268
269
270 for name, fluxNames in self.toCorrect.items():
271
272
273 fluxes = goodRefCat[fluxNames.fluxName]
274 with np.errstate(invalid="ignore"):
275 isGood = (
276 (~goodRefCat[fluxNames.flagName])
277 & (np.isfinite(fluxes))
278 & (fluxes > 0.0)
279 )
280
281
282
283 if (isGood.sum() - 1) < self.config.minDegreesOfFreedom:
284 if name in self.config.allowFailure:
285 self.log.warning("Unable to measure aperture correction for '%s': "
286 "only %d sources, but require at least %d.",
287 name, isGood.sum(), self.config.minDegreesOfFreedom + 1)
288 continue
289 else:
290 raise MeasureApCorrError(name=name, nSources=isGood.sum(),
291 ndof=self.config.minDegreesOfFreedom + 1)
292
293 goodCat = goodRefCat[isGood].copy()
294
295 x = goodCat['slot_Centroid_x']
296 y = goodCat['slot_Centroid_y']
297 z = goodCat[self.refFluxNames.fluxName]/goodCat[fluxNames.fluxName]
298 ids = goodCat['id']
299
300
301
302 fitValues = np.median(z)
303
304 ctrl = self.config.fitConfig.makeControl()
305
306 allBad = False
307 for iteration in range(self.config.numIter):
308 resid = z - fitValues
309
310
311
312 apCorrErr = median_abs_deviation(resid, scale="normal") + 1e-7
313 keep = np.abs(resid) <= self.config.numSigmaClip * apCorrErr
314
315 self.log.debug("Removing %d sources as outliers.", len(resid) - keep.sum())
316
317 x = x[keep]
318 y = y[keep]
319 z = z[keep]
320 ids = ids[keep]
321
322 while (len(x) - ctrl.computeSize()) < self.config.minDegreesOfFreedom:
323 if ctrl.orderX > 0:
324 ctrl.orderX -= 1
325 else:
326 allBad = True
327 break
328 if ctrl.orderY > 0:
329 ctrl.orderY -= 1
330 else:
331 allBad = True
332 break
333
334 if allBad:
335 if name in self.config.allowFailure:
336 self.log.warning("Unable to measure aperture correction for '%s': "
337 "only %d sources remain, but require at least %d." %
338 (name, keep.sum(), self.config.minDegreesOfFreedom + 1))
339 break
340 else:
341 raise MeasureApCorrError(name=name, nSources=keep.sum(),
342 ndof=self.config.minDegreesOfFreedom + 1,
343 iteration=iteration+1)
344
345 apCorrField = ChebyshevBoundedField.fit(bbox, x, y, z, ctrl)
346 fitValues = apCorrField.evaluate(x, y)
347
348 if allBad:
349 continue
350
351 if self.config.doFinalMedianShift:
352 med = np.median(fitValues - z)
353 coeffs = apCorrField.getCoefficients().copy()
354 coeffs[0, 0] -= med
355 apCorrField = ChebyshevBoundedField(bbox, coeffs)
356 fitValues = apCorrField.evaluate(x, y)
357
358 self.log.info(
359 "Aperture correction for %s from %d stars: med %f, MAD %f, RMS %f",
360 name,
361 len(x),
362 np.median(fitValues - z),
363 median_abs_deviation(fitValues - z, scale="normal"),
364 np.mean((fitValues - z)**2.)**0.5,
365 )
366
367 if display:
368 plotApCorr(bbox, x, y, z, apCorrField, "%s, final" % (name,), doPause)
369
370
371 used = np.zeros(len(catalog), dtype=bool)
372 used[np.searchsorted(catalog['id'], ids)] = True
373 catalog[fluxNames.usedName] = used
374
375
376
377
378
379 apCorrMap[fluxNames.fluxName] = apCorrField
380 apCorrMap[fluxNames.errName] = ChebyshevBoundedField(
381 bbox,
382 np.array([[apCorrErr]]),
383 )
384
385 return Struct(
386 apCorrMap=apCorrMap,
387 )
388
389