222 def run(self, catalog, apCorrMap):
223 """Apply aperture corrections to a catalog of sources.
224
225 Parameters
226 ----------
227 catalog : `lsst.afw.table.SourceCatalog`
228 Catalog of sources. Will be updated in place.
229 apCorrMap : `lsst.afw.image.ApCorrMap`
230 Aperture correction map
231
232 Notes
233 -----
234 If you show debug-level log messages then you will see statistics for
235 the effects of aperture correction.
236 """
237 self.log.info("Applying aperture corrections to %d instFlux fields", len(self.apCorrInfoDict))
238 if UseNaiveFluxErr:
239 self.log.debug("Use naive instFlux sigma computation")
240 else:
241 self.log.debug("Use complex instFlux sigma computation that double-counts photon noise "
242 "and thus over-estimates instFlux uncertainty")
243
244
245 periodicLog = PeriodicLogger(self.log)
246
247 for apCorrIndex, apCorrInfo in enumerate(self.apCorrInfoDict.values()):
248 apCorrModel = apCorrMap.get(apCorrInfo.modelName)
249 apCorrErrModel = apCorrMap.get(apCorrInfo.modelSigmaName)
250 if None in (apCorrModel, apCorrErrModel):
251 missingNames = [(apCorrInfo.modelName, apCorrInfo.modelSigmaName)[i]
252 for i, model in enumerate((apCorrModel, apCorrErrModel)) if model is None]
253 self.log.warning("Cannot aperture correct %s because could not find %s in apCorrMap",
254 apCorrInfo.name, " or ".join(missingNames))
255 for source in catalog:
256 source.set(apCorrInfo.apCorrFlagKey, True)
257 continue
258
259
260 catalog[apCorrInfo.apCorrFlagKey] = True
261 oldFluxFlagState = np.zeros(len(catalog), dtype=np.bool_)
262 if self.config.doFlagApCorrFailures:
263 oldFluxFlagState = catalog[apCorrInfo.fluxFlagKey]
264 catalog[apCorrInfo.fluxFlagKey] = True
265
266 apCorr = apCorrModel.evaluate(catalog["slot_Centroid_x"], catalog["slot_Centroid_y"])
267 if not UseNaiveFluxErr:
268 apCorrErr = apCorrErrModel.evaluate(
269 catalog["slot_Centroid_x"],
270 catalog["slot_Centroid_y"],
271 )
272 else:
273 apCorrErr = np.zeros(len(catalog))
274
275 if apCorrInfo.doApCorrColumn:
276 catalog[apCorrInfo.apCorrKey] = apCorr
277 catalog[apCorrInfo.apCorrErrKey] = apCorrErr
278
279
280 good = np.isfinite(apCorr) & (apCorr > 0.0) & (apCorrErr >= 0.0)
281
282 if good.sum() == 0:
283 continue
284
285 instFlux = catalog[apCorrInfo.instFluxKey]
286 instFluxErr = catalog[apCorrInfo.instFluxErrKey]
287 catalog[apCorrInfo.instFluxKey][good] = instFlux[good]*apCorr[good]
288 if UseNaiveFluxErr:
289 catalog[apCorrInfo.instFluxErrKey][good] = instFluxErr[good]*apCorr[good]
290 else:
291 a = instFluxErr[good]/instFlux[good]
292 b = apCorrErr[good]/apCorr[good]
293 err = np.abs(instFlux[good]*apCorr[good])*np.sqrt(a*a + b*b)
294 catalog[apCorrInfo.instFluxErrKey][good] = err
295
296 flags = catalog[apCorrInfo.apCorrFlagKey].copy()
297 flags[good] = False
298 catalog[apCorrInfo.apCorrFlagKey] = flags
299 if self.config.doFlagApCorrFailures:
300 fluxFlags = catalog[apCorrInfo.fluxFlagKey].copy()
301 fluxFlags[good] = oldFluxFlagState[good]
302 catalog[apCorrInfo.fluxFlagKey] = fluxFlags
303
304
305 periodicLog.log("Aperture corrections applied to %d fields out of %d",
306 apCorrIndex + 1, len(self.apCorrInfoDict))
307
308 if self.log.isEnabledFor(self.log.DEBUG):
309
310 apCorrArr = np.array([s.get(apCorrInfo.apCorrKey) for s in catalog])
311 apCorrErrArr = np.array([s.get(apCorrInfo.apCorrErrKey) for s in catalog])
312 self.log.debug("For instFlux field %r: mean apCorr=%s, stdDev apCorr=%s, "
313 "mean apCorrErr=%s, stdDev apCorrErr=%s for %s sources",
314 apCorrInfo.name, apCorrArr.mean(), apCorrArr.std(),
315 apCorrErrArr.mean(), apCorrErrArr.std(), len(catalog))