Coverage for python/lsst/meas/algorithms/measureApCorr.py: 79%
172 statements
« prev ^ index » next coverage.py v7.14.3, created at 2026-06-29 22:25 +0000
1#
2# LSST Data Management System
3#
4# Copyright 2008-2017 AURA/LSST.
5#
6# This product includes software developed by the
7# LSST Project (http://www.lsst.org/).
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 LSST License Statement and
20# the GNU General Public License along with this program. If not,
21# see <https://www.lsstcorp.org/LegalNotices/>.
22#
24__all__ = ("MeasureApCorrConfig", "MeasureApCorrTask", "MeasureApCorrError")
26import numpy as np
27from scipy.stats import median_abs_deviation
29import lsst.pex.config
30from lsst.afw.image import ApCorrMap
31from lsst.afw.math import ChebyshevBoundedField, ChebyshevBoundedFieldConfig
32from lsst.pipe.base import Task, Struct, AlgorithmError
33from lsst.meas.base.apCorrRegistry import getApCorrNameSet
35from .sourceSelector import sourceSelectorRegistry
38class MeasureApCorrError(AlgorithmError):
39 """Raised if Aperture Correction fails in a non-recoverable way.
41 Parameters
42 ----------
43 name : `str`
44 Name of the kind of aperture correction that failed; typically an
45 instFlux catalog field.
46 nSources : `int`
47 Number of sources available to the fitter at the point of failure.
48 ndof : `int`
49 Number of degrees of freedom required at the point of failure.
50 iteration : `int`, optional
51 Which fit iteration the failure occurred in.
52 """
53 def __init__(self, *, name, nSources, ndof, iteration=None):
54 msg = f"Unable to measure aperture correction for '{name}'"
55 if iteration is not None:
56 msg += f" after {iteration} steps:"
57 else:
58 msg += ":"
59 msg += f" only {nSources} sources, but require at least {ndof}."
60 super().__init__(msg)
61 self.name = name
62 self.nSources = nSources
63 self.ndof = ndof
64 self.iteration = iteration
66 @property
67 def metadata(self):
68 metadata = {"name": self.name,
69 "nSources": self.nSources,
70 "ndof": self.ndof,
71 }
72 # NOTE: have to do this because task metadata doesn't allow None.
73 if self.iteration is not None:
74 metadata["iteration"] = self.iteration
75 return metadata
78class _FluxNames:
79 """A collection of flux-related names for a given flux measurement algorithm.
81 Parameters
82 ----------
83 name : `str`
84 Name of flux measurement algorithm, e.g. ``base_PsfFlux``.
85 schema : `lsst.afw.table.Schema`
86 Catalog schema containing the flux field. The ``{name}_instFlux``,
87 ``{name}_instFluxErr``, ``{name}_flag`` fields are checked for
88 existence, and the ``apcorr_{name}_used`` field is added.
90 Raises
91 ------
92 KeyError if any of instFlux, instFluxErr, or flag fields is missing.
93 """
94 def __init__(self, name, schema):
95 self.fluxName = name + "_instFlux"
96 if self.fluxName not in schema:
97 raise KeyError("Could not find " + self.fluxName)
98 self.errName = name + "_instFluxErr"
99 if self.errName not in schema: 99 ↛ 100line 99 didn't jump to line 100 because the condition on line 99 was never true
100 raise KeyError("Could not find " + self.errName)
101 self.flagName = name + "_flag"
102 if self.flagName not in schema: 102 ↛ 103line 102 didn't jump to line 103 because the condition on line 102 was never true
103 raise KeyError("Cound not find " + self.flagName)
104 self.usedName = "apcorr_" + name + "_used"
105 schema.addField(self.usedName, type="Flag",
106 doc="Set if source was used in measuring aperture correction.")
109class MeasureApCorrConfig(lsst.pex.config.Config):
110 """Configuration for MeasureApCorrTask.
111 """
112 refFluxName = lsst.pex.config.Field(
113 doc="Field name prefix for the flux other measurements should be aperture corrected to match",
114 dtype=str,
115 default="slot_CalibFlux",
116 )
117 sourceSelector = sourceSelectorRegistry.makeField(
118 doc="Selector that sets the stars that aperture corrections will be measured from.",
119 default="science",
120 )
121 minDegreesOfFreedom = lsst.pex.config.RangeField(
122 doc="Minimum number of degrees of freedom (# of valid data points - # of parameters);"
123 " if this is exceeded, the order of the fit is decreased (in both dimensions), and"
124 " if we can't decrease it enough, we'll raise ValueError.",
125 dtype=int,
126 default=1,
127 min=1,
128 )
129 fitConfig = lsst.pex.config.ConfigField(
130 doc="Configuration used in fitting the aperture correction fields.",
131 dtype=ChebyshevBoundedFieldConfig,
132 )
133 numIter = lsst.pex.config.Field(
134 doc="Number of iterations for robust MAD sigma clipping.",
135 dtype=int,
136 default=4,
137 )
138 numSigmaClip = lsst.pex.config.Field(
139 doc="Number of robust MAD sigma to do clipping.",
140 dtype=float,
141 default=4.0,
142 )
143 doFinalMedianShift = lsst.pex.config.Field(
144 doc="Do final shift to ensure medians match.",
145 dtype=bool,
146 default=True,
147 )
148 allowFailure = lsst.pex.config.ListField(
149 doc="Allow these measurement algorithms to fail without an exception.",
150 dtype=str,
151 default=[],
152 )
154 def setDefaults(self):
155 selector = self.sourceSelector["science"]
157 selector.doFlags = True
158 selector.doUnresolved = True
159 selector.doSignalToNoise = True
160 selector.doIsolated = False
161 selector.flags.good = []
162 selector.flags.bad = [
163 "base_PixelFlags_flag_edge",
164 "base_PixelFlags_flag_nodata",
165 "base_PixelFlags_flag_interpolatedCenter",
166 "base_PixelFlags_flag_saturatedCenter",
167 "base_PixelFlags_flag_crCenter",
168 "base_PixelFlags_flag_bad",
169 "base_PixelFlags_flag_interpolated",
170 "base_PixelFlags_flag_saturated",
171 ]
172 selector.signalToNoise.minimum = 200.0
173 selector.signalToNoise.maximum = None
174 selector.signalToNoise.fluxField = "base_PsfFlux_instFlux"
175 selector.signalToNoise.errField = "base_PsfFlux_instFluxErr"
177 def validate(self):
178 lsst.pex.config.Config.validate(self)
179 if self.sourceSelector.target.usesMatches: 179 ↛ 180line 179 didn't jump to line 180 because the condition on line 179 was never true
180 raise lsst.pex.config.FieldValidationError(
181 MeasureApCorrConfig.sourceSelector,
182 self,
183 "Star selectors that require matches are not permitted.")
186class MeasureApCorrTask(Task):
187 """Task to measure aperture correction.
189 For every name to correct, a new field apcorr_{name}_used will
190 be added, and will be logged in self.toCorrect.
192 Parameters
193 ----------
194 schema : `lsst.afw.table.Schema`
195 Schema for the input table; will be modified in place to
196 add ``apcorr_{name}_used`` fields.
197 namesToCorrect : `list` [`str`], optional
198 List of names to correct. If `None` then the set of
199 registered fields in lsst.meas.base.getApCorrNameSet()
200 will be used.
201 **kwargs : `dict`
202 Additional kwargs to pass to lsst.pipe.base.Task.__init__()
204 Raises
205 ------
206 MeasureApCorrError if any of the names to correct fails and is
207 not in the config.allowFailure list.
208 """
209 ConfigClass = MeasureApCorrConfig
210 _DefaultName = "measureApCorr"
212 def __init__(self, schema, namesToCorrect=None, **kwargs):
213 Task.__init__(self, **kwargs)
214 self.refFluxNames = _FluxNames(self.config.refFluxName, schema)
215 self.toCorrect = {} # dict of flux field name prefix: FluxKeys instance
216 names = namesToCorrect if namesToCorrect else getApCorrNameSet()
217 for name in sorted(names):
218 try:
219 self.toCorrect[name] = _FluxNames(name, schema)
220 except KeyError:
221 # if a field in the registry is missing, just ignore it.
222 pass
223 self.makeSubtask("sourceSelector")
225 def run(self, exposure, catalog):
226 """Measure aperture correction
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.
239 Returns
240 -------
241 Struct : `lsst.pipe.base.Struct`
242 Contains the following:
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
252 display = lsstDebug.Info(__name__).display
253 doPause = lsstDebug.Info(__name__).doPause
255 self.log.info("Measuring aperture corrections for %d flux fields", len(self.toCorrect))
257 # First, create a subset of the catalog that contains only selected stars
258 # with non-flagged reference fluxes.
259 selected = self.sourceSelector.run(catalog, exposure=exposure)
261 use = (
262 ~selected.sourceCat[self.refFluxNames.flagName]
263 & (np.isfinite(selected.sourceCat[self.refFluxNames.fluxName]))
264 )
265 goodRefCat = selected.sourceCat[use].copy()
267 apCorrMap = ApCorrMap()
269 # Outer loop over the fields we want to correct
270 for name, fluxNames in self.toCorrect.items():
271 # Create a more restricted subset with only the objects where the to-be-correct flux
272 # is not flagged.
273 fluxes = goodRefCat[fluxNames.fluxName]
274 with np.errstate(invalid="ignore"): # suppress NaN warnings.
275 isGood = (
276 (~goodRefCat[fluxNames.flagName])
277 & (np.isfinite(fluxes))
278 & (fluxes > 0.0)
279 )
281 # The 1 is the minimum number of ctrl.computeSize() when the order
282 # drops to 0 in both x and y.
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)
293 goodCat = goodRefCat[isGood].copy()
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']
300 # We start with an initial fit that is the median offset; this
301 # works well in practice.
302 fitValues = np.median(z)
304 ctrl = self.config.fitConfig.makeControl()
306 allBad = False
307 for iteration in range(self.config.numIter):
308 resid = z - fitValues
309 # We add a small (epsilon) amount of floating-point slop because
310 # the median_abs_deviation may give a value that is just larger than 0
311 # even if given a completely flat residual field (as in tests).
312 apCorrErr = median_abs_deviation(resid, scale="normal") + 1e-7
313 keep = np.abs(resid) <= self.config.numSigmaClip * apCorrErr
315 self.log.debug("Removing %d sources as outliers.", len(resid) - keep.sum())
317 x = x[keep]
318 y = y[keep]
319 z = z[keep]
320 ids = ids[keep]
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: 328 ↛ 331line 328 didn't jump to line 331 because the condition on line 328 was always true
329 ctrl.orderY -= 1
330 else:
331 allBad = True
332 break
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)
345 apCorrField = ChebyshevBoundedField.fit(bbox, x, y, z, ctrl)
346 fitValues = apCorrField.evaluate(x, y)
348 if allBad:
349 continue
351 if self.config.doFinalMedianShift: 351 ↛ 358line 351 didn't jump to line 358 because the condition on line 351 was always true
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)
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 )
367 if display: 367 ↛ 368line 367 didn't jump to line 368 because the condition on line 367 was never true
368 plotApCorr(bbox, x, y, z, apCorrField, "%s, final" % (name,), doPause)
370 # Record which sources were used.
371 used = np.zeros(len(catalog), dtype=bool)
372 used[np.searchsorted(catalog['id'], ids)] = True
373 catalog[fluxNames.usedName] = used
375 # Save the result in the output map
376 # The error is constant spatially (we could imagine being
377 # more clever, but we're not yet sure if it's worth the effort).
378 # We save the errors as a 0th-order ChebyshevBoundedField
379 apCorrMap[fluxNames.fluxName] = apCorrField
380 apCorrMap[fluxNames.errName] = ChebyshevBoundedField(
381 bbox,
382 np.array([[apCorrErr]]),
383 )
385 return Struct(
386 apCorrMap=apCorrMap,
387 )
390def plotApCorr(bbox, xx, yy, zzMeasure, field, title, doPause):
391 """Plot aperture correction fit residuals
393 There are two subplots: residuals against x and y.
395 Intended for debugging.
397 Parameters
398 ----------
399 bbox : `lsst.geom.Box2I`
400 Bounding box (for bounds)
401 xx, yy : `numpy.ndarray`, (N)
402 x and y coordinates
403 zzMeasure : `float`
404 Measured value of the aperture correction
405 field : `lsst.afw.math.ChebyshevBoundedField`
406 Fit aperture correction field
407 title : 'str'
408 Title for plot
409 doPause : `bool`
410 Pause to inspect the residuals plot? If
411 False, there will be a 4 second delay to
412 allow for inspection of the plot before
413 closing it and moving on.
414 """
415 import matplotlib.pyplot as plt
417 zzFit = field.evaluate(xx, yy)
418 residuals = zzMeasure - zzFit
420 fig, axes = plt.subplots(2, 1)
422 axes[0].scatter(xx, residuals, s=3, marker='o', lw=0, alpha=0.7)
423 axes[1].scatter(yy, residuals, s=3, marker='o', lw=0, alpha=0.7)
424 for ax in axes:
425 ax.set_ylabel("ApCorr Fit Residual")
426 ax.set_ylim(0.9*residuals.min(), 1.1*residuals.max())
427 axes[0].set_xlabel("x")
428 axes[0].set_xlim(bbox.getMinX(), bbox.getMaxX())
429 axes[1].set_xlabel("y")
430 axes[1].set_xlim(bbox.getMinY(), bbox.getMaxY())
431 plt.suptitle(title)
433 if not doPause:
434 try:
435 plt.pause(4)
436 plt.close()
437 except Exception:
438 print("%s: plt.pause() failed. Please close plots when done." % __name__)
439 plt.show()
440 else:
441 print("%s: Please close plots when done." % __name__)
442 plt.show()