Coverage for python / lsst / meas / base / applyApCorr.py: 17%

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 astropy.table 

23import numpy as np 

24 

25import lsst.afw.image 

26import lsst.pex.config 

27import lsst.pex.exceptions 

28import lsst.pipe.base 

29from lsst.utils.logging import PeriodicLogger 

30from .apCorrRegistry import getApCorrNameSet 

31 

32# If True then scale instFlux error by apCorr; if False then use a more complex computation 

33# that over-estimates instFlux error (often grossly so) because it double-counts photon noise. 

34# This flag is intended to be temporary until we figure out a better way to compute 

35# the effects of aperture correction on instFlux uncertainty 

36UseNaiveFluxErr = True 

37 

38__all__ = ("ApplyApCorrConfig", "ApplyApCorrTask") 

39 

40 

41class ApCorrInfo: 

42 """Catalog field names and keys needed to aperture correct a particular 

43 instrument flux. 

44 

45 Parameters 

46 ---------- 

47 schema : `lsst.afw.table` 

48 Source catalog schema. Three fields are used to generate keys: 

49 - ``{name}_instFlux`` 

50 - ``{name}_instFluxErr`` 

51 - ``{name}_flag`` 

52 Three fields are added: 

53 - ``{name}_apCorr`` (only if not already added by proxy) 

54 - ``{name}_apCorrErr`` (only if not already added by proxy) 

55 - ``{name}_flag_apCorr`` 

56 model : `str` 

57 Field name prefix for instFlux with aperture correction model, e.g. 

58 "base_PsfFlux" 

59 name : `str` 

60 Field name prefix for instFlux needing aperture correction; may be 

61 `None` if it is the same as ``model`` 

62 

63 Notes 

64 ----- 

65 The aperture correction can be derived from the meaasurements in the 

66 column being aperture-corrected or from measurements in a different 

67 column (a "proxy"). In the first case, we will add columns to contain 

68 the aperture correction values; in the second case (using a proxy), 

69 we will add an alias to the proxy's aperture correction values. In 

70 all cases, we add a flag. 

71 

72 Attributes 

73 ---------- 

74 name : str 

75 Field name prefix for flux needing aperture correction. 

76 modelName : str 

77 Field name for aperture correction model for flux. 

78 modelSigmaName : str 

79 Field name for aperture correction model for fluxErr. 

80 doApCorrColumn : bool 

81 Should we write the aperture correction values? 

82 They should not be written if they're already being written by a proxy. 

83 instFluxName : str 

84 Name of ``instFlux`` field. 

85 instFluxErrName : str 

86 Name of ``instFlux`` sigma field. 

87 """ 

88 

89 name = None 

90 modelName = None 

91 modelSigmaName = None 

92 doApCorrColumn = None 

93 instFluxName = None 

94 instFluxErrName = None 

95 apCorrName = None 

96 

97 def __init__(self, model, schema=None, name=None): 

98 if name is None: 

99 name = model 

100 self.name = name 

101 self.modelName = model + "_instFlux" 

102 self.modelSigmaName = model + "_instFluxErr" 

103 self.instFluxName = name + "_instFlux" 

104 self.instFluxErrName = name + "_instFluxErr" 

105 self.fluxFlagName = name + "_flag" 

106 self.apCorrName = name + "_apCorr" 

107 self.apCorrErrName = name + "_apCorrErr" 

108 self.apCorrFlagName = name + "_flag_apCorr" 

109 

110 # No need to write the same aperture corrections multiple times 

111 if schema is not None: 

112 if name == model or model + "_apCorr" not in schema: 

113 schema.addField( 

114 name + "_apCorr", 

115 doc="aperture correction applied to %s" % (name,), 

116 type=np.float64, 

117 ) 

118 schema.addField( 

119 name + "_apCorrErr", 

120 doc="standard deviation of aperture correction applied to %s" % (name,), 

121 type=np.float64, 

122 ) 

123 else: 

124 aliases = schema.getAliasMap() 

125 aliases.set(name + "_apCorr", model + "_apCorr") 

126 aliases.set(name + "_apCorrErr", model + "_apCorrErr") 

127 

128 if name + "_flag_apCorr" not in schema: 

129 schema.addField( 

130 name + "_flag_apCorr", 

131 doc="set if unable to aperture correct %s" % (name,), 

132 type="Flag", 

133 ) 

134 

135 

136class ApplyApCorrConfig(lsst.pex.config.Config): 

137 """Aperture correction configuration. 

138 """ 

139 

140 ignoreList = lsst.pex.config.ListField( 

141 doc="flux measurement algorithms in getApCorrNameSet() to ignore; " 

142 "if a name is listed that does not appear in getApCorrNameSet() then a warning is logged", 

143 dtype=str, 

144 optional=False, 

145 default=(), 

146 ) 

147 doFlagApCorrFailures = lsst.pex.config.Field( 

148 doc="set the general failure flag for a flux when it cannot be aperture-corrected?", 

149 dtype=bool, 

150 default=True, 

151 ) 

152 proxies = lsst.pex.config.DictField( 

153 doc="flux measurement algorithms to be aperture-corrected by reference to another algorithm; " 

154 "this is a mapping alg1:alg2, where 'alg1' is the algorithm being corrected, and 'alg2' " 

155 "is the algorithm supplying the corrections", 

156 keytype=str, 

157 itemtype=str, 

158 default={}, 

159 ) 

160 xColumn = lsst.pex.config.Field( 

161 doc="name of the x coordinate column in the catalog", 

162 dtype=str, 

163 default="slot_Centroid_x", 

164 ) 

165 yColumn = lsst.pex.config.Field( 

166 doc="name of the y coordinate column in the catalog", 

167 dtype=str, 

168 default="slot_Centroid_y", 

169 ) 

170 

171 

172class ApplyApCorrTask(lsst.pipe.base.Task): 

173 """Apply aperture corrections. 

174 

175 Parameters 

176 ---------- 

177 schema : `lsst.afw.table.Schema` 

178 """ 

179 ConfigClass = ApplyApCorrConfig 

180 _DefaultName = "applyApCorr" 

181 

182 def __init__(self, schema=None, **kwds): 

183 lsst.pipe.base.Task.__init__(self, **kwds) 

184 

185 self.apCorrInfoDict = dict() 

186 apCorrNameSet = getApCorrNameSet() 

187 ignoreSet = set(self.config.ignoreList) 

188 missingNameSet = ignoreSet - set(apCorrNameSet) 

189 if missingNameSet: 

190 self.log.warning("Fields in ignoreList that are not in fluxCorrectList: %s", 

191 sorted(missingNameSet)) 

192 for name in sorted(apCorrNameSet - ignoreSet): 

193 if schema is not None and name + "_instFlux" not in schema: 

194 # if a field in the registry is missing from the schema, silently ignore it. 

195 continue 

196 self.apCorrInfoDict[name] = ApCorrInfo(schema=schema, model=name) 

197 

198 for name, model in self.config.proxies.items(): 

199 if name in apCorrNameSet: 

200 # Already done or ignored 

201 continue 

202 if schema is not None and name + "_instFlux" not in schema: 

203 # Silently ignore 

204 continue 

205 self.apCorrInfoDict[name] = ApCorrInfo(schema=schema, model=model, name=name) 

206 

207 def run(self, catalog, apCorrMap): 

208 """Apply aperture corrections to a catalog of sources. 

209 

210 Parameters 

211 ---------- 

212 catalog : `lsst.afw.table.SourceCatalog` or `astropy.table.Table` 

213 Catalog of sources. Will be updated in place. 

214 apCorrMap : `lsst.afw.image.ApCorrMap` 

215 Aperture correction map 

216 

217 Notes 

218 ----- 

219 If you show debug-level log messages then you will see statistics for 

220 the effects of aperture correction. 

221 """ 

222 self.log.info("Applying aperture corrections to %d instFlux fields", len(self.apCorrInfoDict)) 

223 if UseNaiveFluxErr: 

224 self.log.debug("Use naive instFlux sigma computation") 

225 else: 

226 self.log.debug("Use complex instFlux sigma computation that double-counts photon noise " 

227 "and thus over-estimates instFlux uncertainty") 

228 

229 # Wrap the task logger to a periodic logger. 

230 periodicLog = PeriodicLogger(self.log) 

231 

232 if isinstance(catalog, astropy.table.Table): 

233 # Prepare the astropy table for use with this task. 

234 self._prepAstropyTable(catalog) 

235 xColumn = self.config.xColumn 

236 yColumn = self.config.yColumn 

237 

238 for apCorrIndex, apCorrInfo in enumerate(self.apCorrInfoDict.values()): 

239 apCorrModel = apCorrMap.get(apCorrInfo.modelName) 

240 apCorrErrModel = apCorrMap.get(apCorrInfo.modelSigmaName) 

241 if None in (apCorrModel, apCorrErrModel): 

242 missingNames = [(apCorrInfo.modelName, apCorrInfo.modelSigmaName)[i] 

243 for i, model in enumerate((apCorrModel, apCorrErrModel)) if model is None] 

244 self.log.warning("Cannot aperture correct %s because could not find %s in apCorrMap", 

245 apCorrInfo.name, " or ".join(missingNames)) 

246 catalog[apCorrInfo.apCorrFlagName] = np.ones(len(catalog), dtype=np.bool_) 

247 continue 

248 

249 # Say we've failed before we start; we'll unset these flags on success. 

250 catalog[apCorrInfo.apCorrFlagName] = True 

251 oldFluxFlagState = np.zeros(len(catalog), dtype=np.bool_) 

252 if self.config.doFlagApCorrFailures: 

253 oldFluxFlagState = catalog[apCorrInfo.fluxFlagName] 

254 catalog[apCorrInfo.fluxFlagName] = True 

255 

256 apCorr = apCorrModel.evaluate(catalog[xColumn], catalog[yColumn]) 

257 if not UseNaiveFluxErr: 

258 apCorrErr = apCorrErrModel.evaluate( 

259 catalog[xColumn], 

260 catalog[yColumn], 

261 ) 

262 else: 

263 apCorrErr = np.zeros(len(catalog)) 

264 

265 if apCorrInfo.doApCorrColumn: 

266 catalog[apCorrInfo.apCorrName] = apCorr 

267 catalog[apCorrInfo.apCorrErrName] = apCorrErr 

268 

269 # Check bad values that are not finite (possible for coadds). 

270 good = np.isfinite(apCorr) & (apCorr > 0.0) & (apCorrErr >= 0.0) 

271 

272 if good.sum() == 0: 

273 continue 

274 

275 instFlux = catalog[apCorrInfo.instFluxName] 

276 instFluxErr = catalog[apCorrInfo.instFluxErrName] 

277 catalog[apCorrInfo.instFluxName][good] = instFlux[good]*apCorr[good] 

278 if UseNaiveFluxErr: 

279 catalog[apCorrInfo.instFluxErrName][good] = instFluxErr[good]*apCorr[good] 

280 else: 

281 a = instFluxErr[good]/instFlux[good] 

282 b = apCorrErr[good]/apCorr[good] 

283 err = np.abs(instFlux[good]*apCorr[good])*np.sqrt(a*a + b*b) 

284 catalog[apCorrInfo.instFluxErrName][good] = err 

285 

286 flags = catalog[apCorrInfo.apCorrFlagName].copy() 

287 flags[good] = False 

288 catalog[apCorrInfo.apCorrFlagName] = flags 

289 if self.config.doFlagApCorrFailures: 

290 fluxFlags = catalog[apCorrInfo.fluxFlagName].copy() 

291 fluxFlags[good] = oldFluxFlagState[good] 

292 catalog[apCorrInfo.fluxFlagName] = fluxFlags 

293 

294 # Log a message if it has been a while since the last log. 

295 periodicLog.log("Aperture corrections applied to %d fields out of %d", 

296 apCorrIndex + 1, len(self.apCorrInfoDict)) 

297 

298 if self.log.isEnabledFor(self.log.DEBUG): 

299 # log statistics on the effects of aperture correction 

300 apCorrArr = np.asarray([s[apCorrInfo.instFluxName] for s in catalog]) 

301 apCorrErrArr = np.asarray([s[apCorrInfo.instFluxErrName] for s in catalog]) 

302 self.log.debug("For instFlux field %r: mean apCorr=%s, stdDev apCorr=%s, " 

303 "mean apCorrErr=%s, stdDev apCorrErr=%s for %s sources", 

304 apCorrInfo.name, apCorrArr.mean(), apCorrArr.std(), 

305 apCorrErrArr.mean(), apCorrErrArr.std(), len(catalog)) 

306 

307 def _prepAstropyTable(self, table): 

308 """Prepare an astropy table for use with this task. 

309 

310 Parameters 

311 ---------- 

312 table : `astropy.table.Table` 

313 Table to prepare. 

314 

315 Returns 

316 ------- 

317 `astropy.table.Table` 

318 Prepared table. 

319 """ 

320 for apCorrInfo in self.apCorrInfoDict.values(): 

321 if apCorrInfo.apCorrName not in table.colnames: 

322 table[apCorrInfo.apCorrName] = np.zeros(len(table), dtype=np.float64) 

323 if apCorrInfo.apCorrErrName not in table.colnames: 

324 table[apCorrInfo.apCorrErrName] = np.zeros(len(table), dtype=np.float64) 

325 if apCorrInfo.apCorrFlagName not in table.colnames: 

326 table[apCorrInfo.apCorrFlagName] = np.zeros(len(table), dtype=bool) 

327 if apCorrInfo.fluxFlagName not in table.colnames: 

328 table[apCorrInfo.fluxFlagName] = np.zeros(len(table), dtype=bool)