LSST Applications g063fba187b+cac8b7c890,g0f08755f38+6aee506743,g1653933729+a8ce1bb630,g168dd56ebc+a8ce1bb630,g1a2382251a+b4475c5878,g1dcb35cd9c+8f9bc1652e,g20f6ffc8e0+6aee506743,g217e2c1bcf+73dee94bd0,g28da252d5a+1f19c529b9,g2bbee38e9b+3f2625acfc,g2bc492864f+3f2625acfc,g3156d2b45e+6e55a43351,g32e5bea42b+1bb94961c2,g347aa1857d+3f2625acfc,g35bb328faa+a8ce1bb630,g3a166c0a6a+3f2625acfc,g3e281a1b8c+c5dd892a6c,g3e8969e208+a8ce1bb630,g414038480c+5927e1bc1e,g41af890bb2+8a9e676b2a,g7af13505b9+809c143d88,g80478fca09+6ef8b1810f,g82479be7b0+f568feb641,g858d7b2824+6aee506743,g89c8672015+f4add4ffd5,g9125e01d80+a8ce1bb630,ga5288a1d22+2903d499ea,gb58c049af0+d64f4d3760,gc28159a63d+3f2625acfc,gcab2d0539d+b12535109e,gcf0d15dbbd+46a3f46ba9,gda6a2b7d83+46a3f46ba9,gdaeeff99f8+1711a396fd,ge79ae78c31+3f2625acfc,gef2f8181fd+0a71e47438,gf0baf85859+c1f95f4921,gfa517265be+6aee506743,gfa999e8aa5+17cd334064,w.2024.51
LSST Data Management Base Package
Loading...
Searching...
No Matches
applyApCorr.py
Go to the documentation of this file.
1# This file is part of meas_base.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
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 GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
21
22import numpy as np
23
24import lsst.afw.image
25import lsst.pex.config
27import lsst.pipe.base
28from lsst.utils.logging import PeriodicLogger
29from .apCorrRegistry import getApCorrNameSet
30
31# If True then scale instFlux error by apCorr; if False then use a more complex computation
32# that over-estimates instFlux error (often grossly so) because it double-counts photon noise.
33# This flag is intended to be temporary until we figure out a better way to compute
34# the effects of aperture correction on instFlux uncertainty
35UseNaiveFluxErr = True
36
37__all__ = ("ApplyApCorrConfig", "ApplyApCorrTask")
38
39
41 """Catalog field names and keys needed to aperture correct a particular
42 instrument flux.
43
44 Parameters
45 ----------
46 schema : `lsst.afw.table`
47 Source catalog schema. Three fields are used to generate keys:
48 - ``{name}_instFlux``
49 - ``{name}_instFluxErr``
50 - ``{name}_flag``
51 Three fields are added:
52 - ``{name}_apCorr`` (only if not already added by proxy)
53 - ``{name}_apCorrErr`` (only if not already added by proxy)
54 - ``{name}_flag_apCorr``
55 model : `str`
56 Field name prefix for instFlux with aperture correction model, e.g.
57 "base_PsfFlux"
58 name : `str`
59 Field name prefix for instFlux needing aperture correction; may be
60 `None` if it is the same as ``model``
61
62 Notes
63 -----
64 The aperture correction can be derived from the meaasurements in the
65 column being aperture-corrected or from measurements in a different
66 column (a "proxy"). In the first case, we will add columns to contain
67 the aperture correction values; in the second case (using a proxy),
68 we will add an alias to the proxy's aperture correction values. In
69 all cases, we add a flag.
70 """
71
72 name = None
73 """Field name prefix for flux needing aperture correction (`str`).
74 """
75
76 modelName = None
77 """Field name for aperture correction model for flux (`str`).
78 """
79
80 modelSigmaName = None
81 """Field name for aperture correction model for fluxErr (`str`).
82 """
83
84 doApCorrColumn = None
85 """Should we write the aperture correction values (`bool`)?
86
87 They should not be written if they're already being written by a proxy.
88 """
89
90 instFluxName = None
91 """Name of ``instFlux`` field (`str`).
92 """
93
94 instFluxErrName = None
95 """Name of ``instFlux`` sigma field (`str`).
96 """
97
98 instFluxKey = None
99 """Key to ``instFlux`` field (`lsst.afw.table.schema.Key`).
100 """
101
102 instFluxErrKey = None
103 """Key to ``instFlux`` sigma field (`lsst.afw.table.schema.Key`).
104 """
105
106 fluxFlagKey = None
107 """Key to the flux flag field (`lsst.afw.table.schema.Key`).
108 """
109
110 apCorrKey = None
111 """Key to new aperture correction field (`lsst.afw.table.schema.Key`).
112 """
113
114 apCorrErrKey = None
115 """Key to new aperture correction sigma field (`lsst.afw.table.schema.Key`).
116 """
117
118 apCorrFlagKey = None
119 """Key to new aperture correction flag field (`lsst.afw.table.schema.Key`).
120 """
121
122 def __init__(self, schema, model, name=None):
123 if name is None:
124 name = model
125 self.name = name
126 self.modelName = model + "_instFlux"
127 self.modelSigmaName = model + "_instFluxErr"
128 self.instFluxName = name + "_instFlux"
129 self.instFluxErrName = name + "_instFluxErr"
130 self.instFluxKey = schema.find(self.instFluxName).key
131 self.instFluxErrKey = schema.find(self.instFluxErrName).key
132 self.fluxFlagKey = schema.find(name + "_flag").key
133
134 # No need to write the same aperture corrections multiple times
135 self.doApCorrColumn = (name == model or model + "_apCorr" not in schema)
136 if self.doApCorrColumn:
137 self.apCorrKey = schema.addField(
138 name + "_apCorr",
139 doc="aperture correction applied to %s" % (name,),
140 type=np.float64,
141 )
142 self.apCorrErrKey = schema.addField(
143 name + "_apCorrErr",
144 doc="standard deviation of aperture correction applied to %s" % (name,),
145 type=np.float64,
146 )
147 else:
148 aliases = schema.getAliasMap()
149 aliases.set(name + "_apCorr", model + "_apCorr")
150 aliases.set(name + "_apCorrErr", model + "_apCorrErr")
151 self.apCorrKey = schema.find(name + "_apCorr").key
152 self.apCorrErrKey = schema.find(name + "_apCorrErr").key
153
154 self.apCorrFlagKey = schema.addField(
155 name + "_flag_apCorr",
156 doc="set if unable to aperture correct %s" % (name,),
157 type="Flag",
158 )
159
160
162 """Aperture correction configuration.
163 """
164
166 doc="flux measurement algorithms in getApCorrNameSet() to ignore; "
167 "if a name is listed that does not appear in getApCorrNameSet() then a warning is logged",
168 dtype=str,
169 optional=False,
170 default=(),
171 )
172 doFlagApCorrFailures = lsst.pex.config.Field(
173 doc="set the general failure flag for a flux when it cannot be aperture-corrected?",
174 dtype=bool,
175 default=True,
176 )
178 doc="flux measurement algorithms to be aperture-corrected by reference to another algorithm; "
179 "this is a mapping alg1:alg2, where 'alg1' is the algorithm being corrected, and 'alg2' "
180 "is the algorithm supplying the corrections",
181 keytype=str,
182 itemtype=str,
183 default={},
184 )
185
186
187class ApplyApCorrTask(lsst.pipe.base.Task):
188 """Apply aperture corrections.
189
190 Parameters
191 ----------
192 schema : `lsst.afw.table.Schema`
193 """
194 ConfigClass = ApplyApCorrConfig
195 _DefaultName = "applyApCorr"
196
197 def __init__(self, schema, **kwds):
198 lsst.pipe.base.Task.__init__(self, **kwds)
199
200 self.apCorrInfoDict = dict()
201 apCorrNameSet = getApCorrNameSet()
202 ignoreSet = set(self.config.ignoreList)
203 missingNameSet = ignoreSet - set(apCorrNameSet)
204 if missingNameSet:
205 self.log.warning("Fields in ignoreList that are not in fluxCorrectList: %s",
206 sorted(missingNameSet))
207 for name in sorted(apCorrNameSet - ignoreSet):
208 if name + "_instFlux" not in schema:
209 # if a field in the registry is missing from the schema, silently ignore it.
210 continue
211 self.apCorrInfoDict[name] = ApCorrInfo(schema=schema, model=name)
212
213 for name, model in self.config.proxies.items():
214 if name in apCorrNameSet:
215 # Already done or ignored
216 continue
217 if name + "_instFlux" not in schema:
218 # Silently ignore
219 continue
220 self.apCorrInfoDict[name] = ApCorrInfo(schema=schema, model=model, name=name)
221
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 # Wrap the task logger to a periodic logger.
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 # Say we've failed before we start; we'll unset these flags on success.
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 # Check bad values that are not finite (possible for coadds).
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 # Log a message if it has been a while since the last log.
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 # log statistics on the effects of aperture correction
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))
__init__(self, schema, model, name=None)