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