Coverage for python / lsst / meas / algorithms / measureApCorr.py: 15%
172 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-04-17 09:07 +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:
100 raise KeyError("Could not find " + self.errName)
101 self.flagName = name + "_flag"
102 if self.flagName not in schema:
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:
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:
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:
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:
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()