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 if schema is not None: 

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

112 self.doApCorrColumn = (name == model or model + "_apCorr" not in schema) 

113 if self.doApCorrColumn: 

114 schema.addField( 

115 name + "_apCorr", 

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

117 type=np.float64, 

118 ) 

119 schema.addField( 

120 name + "_apCorrErr", 

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

122 type=np.float64, 

123 ) 

124 else: 

125 aliases = schema.getAliasMap() 

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

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

128 

129 if name + "_flag_apCorr" not in schema: 

130 schema.addField( 

131 name + "_flag_apCorr", 

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

133 type="Flag", 

134 ) 

135 else: 

136 # Always write the aperture corrections for astropy/arrow tables. 

137 self.doApCorrColumn = True 

138 

139 

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

141 """Aperture correction configuration. 

142 """ 

143 

144 ignoreList = lsst.pex.config.ListField( 

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

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

147 dtype=str, 

148 optional=False, 

149 default=(), 

150 ) 

151 doFlagApCorrFailures = lsst.pex.config.Field( 

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

153 dtype=bool, 

154 default=True, 

155 ) 

156 proxies = lsst.pex.config.DictField( 

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

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

159 "is the algorithm supplying the corrections", 

160 keytype=str, 

161 itemtype=str, 

162 default={}, 

163 ) 

164 xColumn = lsst.pex.config.Field( 

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

166 dtype=str, 

167 default="slot_Centroid_x", 

168 ) 

169 yColumn = lsst.pex.config.Field( 

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

171 dtype=str, 

172 default="slot_Centroid_y", 

173 ) 

174 

175 

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

177 """Apply aperture corrections. 

178 

179 Parameters 

180 ---------- 

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

182 """ 

183 ConfigClass = ApplyApCorrConfig 

184 _DefaultName = "applyApCorr" 

185 

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

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

188 

189 self.apCorrInfoDict = dict() 

190 apCorrNameSet = getApCorrNameSet() 

191 ignoreSet = set(self.config.ignoreList) 

192 missingNameSet = ignoreSet - set(apCorrNameSet) 

193 if missingNameSet: 

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

195 sorted(missingNameSet)) 

196 for name in sorted(apCorrNameSet - ignoreSet): 

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

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

199 continue 

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

201 

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

203 if name in apCorrNameSet: 

204 # Already done or ignored 

205 continue 

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

207 # Silently ignore 

208 continue 

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

210 

211 def run(self, catalog, apCorrMap): 

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

213 

214 Parameters 

215 ---------- 

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

217 Catalog of sources. Will be updated in place. 

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

219 Aperture correction map 

220 

221 Notes 

222 ----- 

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

224 the effects of aperture correction. 

225 """ 

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

227 if UseNaiveFluxErr: 

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

229 else: 

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

231 "and thus over-estimates instFlux uncertainty") 

232 

233 # Wrap the task logger to a periodic logger. 

234 periodicLog = PeriodicLogger(self.log) 

235 

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

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

238 self._prepAstropyTable(catalog) 

239 xColumn = self.config.xColumn 

240 yColumn = self.config.yColumn 

241 

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

243 apCorrModel = apCorrMap.get(apCorrInfo.modelName) 

244 apCorrErrModel = apCorrMap.get(apCorrInfo.modelSigmaName) 

245 if None in (apCorrModel, apCorrErrModel): 

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

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

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

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

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

251 continue 

252 

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

254 catalog[apCorrInfo.apCorrFlagName] = True 

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

256 if self.config.doFlagApCorrFailures: 

257 oldFluxFlagState = catalog[apCorrInfo.fluxFlagName] 

258 catalog[apCorrInfo.fluxFlagName] = True 

259 

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

261 if not UseNaiveFluxErr: 

262 apCorrErr = apCorrErrModel.evaluate( 

263 catalog[xColumn], 

264 catalog[yColumn], 

265 ) 

266 else: 

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

268 

269 if apCorrInfo.doApCorrColumn: 

270 catalog[apCorrInfo.apCorrName] = apCorr 

271 catalog[apCorrInfo.apCorrErrName] = apCorrErr 

272 

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

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

275 

276 if good.sum() == 0: 

277 continue 

278 

279 instFlux = catalog[apCorrInfo.instFluxName] 

280 instFluxErr = catalog[apCorrInfo.instFluxErrName] 

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

282 if UseNaiveFluxErr: 

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

284 else: 

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

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

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

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

289 

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

291 flags[good] = False 

292 catalog[apCorrInfo.apCorrFlagName] = flags 

293 if self.config.doFlagApCorrFailures: 

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

295 fluxFlags[good] = oldFluxFlagState[good] 

296 catalog[apCorrInfo.fluxFlagName] = fluxFlags 

297 

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

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

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

301 

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

303 # log statistics on the effects of aperture correction 

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

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

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

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

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

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

310 

311 def _prepAstropyTable(self, table): 

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

313 

314 Parameters 

315 ---------- 

316 table : `astropy.table.Table` 

317 Table to prepare. 

318 

319 Returns 

320 ------- 

321 `astropy.table.Table` 

322 Prepared table. 

323 """ 

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

325 if apCorrInfo.apCorrName not in table.colnames: 

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

327 if apCorrInfo.apCorrErrName not in table.colnames: 

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

329 if apCorrInfo.apCorrFlagName not in table.colnames: 

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

331 if apCorrInfo.fluxFlagName not in table.colnames: 

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