Coverage for python/lsst/pipe/tasks/measurePsf.py: 16%

1# This file is part of pipe_tasks.

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.

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/>.

22__all__ = ["MeasurePsfConfig", "MeasurePsfTask"]

24import numpy as np

26import lsst.afw.display as afwDisplay

27import lsst.afw.math as afwMath

28import lsst.meas.algorithms as measAlg

29import lsst.meas.algorithms.utils as maUtils

30import lsst.pex.config as pexConfig

31import lsst.pipe.base as pipeBase

32import lsst.meas.extensions.psfex.psfexPsfDeterminer # noqa: F401

33from lsst.utils.timer import timeMethod

36class NonfinitePsfShapeError(pipeBase.AlgorithmError):

37 """Raised if the radius of the fitted psf model is non-finite.

39 Parameters

40 ----------

41 psf_size : `float`

42 Fitted size of the failing PSF model.

43 """

44 def __init__(self, psf_size) -> None:

45 self._psf_size = psf_size

46 super().__init__(

47 f"Failed to determine PSF: fitter returned model with non-finite PSF size ({psf_size} pixels)."

48 )

50 @property

51 def metadata(self) -> dict:

52 return {

53 "psf_size": self._psf_size,

54 }

57class MeasurePsfConfig(pexConfig.Config):

58 starSelector = measAlg.sourceSelectorRegistry.makeField(

59 "Star selection algorithm",

60 default="objectSize"

61 )

62 makePsfCandidates = pexConfig.ConfigurableField(

63 target=measAlg.MakePsfCandidatesTask,

64 doc="Task to make psf candidates from selected stars.",

65 )

66 psfDeterminer = measAlg.psfDeterminerRegistry.makeField(

67 "PSF Determination algorithm",

68 default="psfex"

69 )

70 reserve = pexConfig.ConfigurableField(

71 target=measAlg.ReserveSourcesTask,

72 doc="Reserve sources from fitting"

73 )

75 def setDefaults(self):

76 super().setDefaults()

77 if self.psfDeterminer.name == "piff" and self.psfDeterminer["piff"].useCoordinates == "sky":

78 self.makePsfCandidates.kernelSize = 35

80 def validate(self):

81 super().validate()

83 match self.psfDeterminer.name:

84 # Perform Piff-specific validations.

85 case "piff":

86 if (

87 self.psfDeterminer["piff"].stampSize

88 and self.psfDeterminer["piff"].stampSize

89 > self.makePsfCandidates.kernelSize

90 ):

91 # Allowing this would result in a cutout size globally.

92 msg = (

93 f"PIFF stampSize={self.psfDeterminer['piff'].stampSize}"

94 f" must be >= psf candidate kernelSize={self.makePsfCandidates.kernelSize}."

95 )

96 raise pexConfig.FieldValidationError(

97 MeasurePsfConfig.makePsfCandidates, self, msg

98 )

100 model_size = self.psfDeterminer["piff"].modelSize

101 sampling_size = self.psfDeterminer["piff"].samplingSize

102 # Calculate the minimum cutout size for stars, given how large

103 # the PSF model is and the sampling size.

104 if self.psfDeterminer["piff"].useCoordinates == "sky":

105 min_kernel_size = int(

106 1.415 * model_size / sampling_size

107 ) # 1.415 = sqrt(2)

108 else:

109 min_kernel_size = int(model_size / sampling_size)

110

111 if self.makePsfCandidates.kernelSize < min_kernel_size:

112 msg = (

113 f"psf candidate kernelSize={self.makePsfCandidates.kernelSize}"

114 f" must be >= {min_kernel_size} for PIFF modelSize={model_size}."

115 )

116 raise pexConfig.FieldValidationError(

117 MeasurePsfConfig.makePsfCandidates, self, msg

118 )

119

120

121class MeasurePsfTask(pipeBase.Task):

122 """A task that selects stars from a catalog of sources and uses those to measure the PSF.

123

124 Parameters

125 ----------

126 schema : `lsst.sfw.table.Schema`

127 An `lsst.afw.table.Schema` used to create the output `lsst.afw.table.SourceCatalog`.

128 **kwargs :

129 Keyword arguments passed to lsst.pipe.base.task.Task.__init__.

130

131 Notes

132 -----

133 If schema is not None, 'calib_psf_candidate' and 'calib_psf_used' fields will be added to

134 identify which stars were employed in the PSF estimation.

135

136 This task can add fields to the schema, so any code calling this task must ensure that

137 these fields are indeed present in the input table.

138

139 The star selector is a subclass of

140 ``lsst.meas.algorithms.starSelector.BaseStarSelectorTask`` "lsst.meas.algorithms.BaseStarSelectorTask"

141 and the PSF determiner is a sublcass of

142 ``lsst.meas.algorithms.psfDeterminer.BasePsfDeterminerTask`` "lsst.meas.algorithms.BasePsfDeterminerTask"

143

144 There is no establised set of configuration parameters for these algorithms, so once you start modifying

145 parameters (as we do in @ref pipe_tasks_measurePsf_Example) your code is no longer portable.

146

147 Debugging:

148

149 .. code-block:: none

150

151 display

152 If True, display debugging plots

153 displayExposure

154 display the Exposure + spatialCells

155 displayPsfCandidates

156 show mosaic of candidates

157 showBadCandidates

158 Include bad candidates

159 displayPsfMosaic

160 show mosaic of reconstructed PSF(xy)

161 displayResiduals

162 show residuals

163 normalizeResiduals

164 Normalise residuals by object amplitude

165

166 Additionally you can enable any debug outputs that your chosen star selector and psf determiner support.

167 """

168 ConfigClass = MeasurePsfConfig

169 _DefaultName = "measurePsf"

170

171 def __init__(self, schema=None, **kwargs):

172 pipeBase.Task.__init__(self, **kwargs)

173 if schema is not None:

174 self.candidateKey = schema.addField(

175 "calib_psf_candidate", type="Flag",

176 doc=("Flag set if the source was a candidate for PSF determination, "

177 "as determined by the star selector.")

178 )

179 self.usedKey = schema.addField(

180 "calib_psf_used", type="Flag",

181 doc=("Flag set if the source was actually used for PSF determination, "

182 "as determined by the '%s' PSF determiner.") % self.config.psfDeterminer.name

183 )

184 else:

185 self.candidateKey = None

186 self.usedKey = None

187 self.makeSubtask("starSelector")

188 self.makeSubtask("makePsfCandidates")

189 self.makeSubtask("psfDeterminer", schema=schema)

190 self.makeSubtask("reserve", columnName="calib_psf", schema=schema,

191 doc="set if source was reserved from PSF determination")

192

193 @timeMethod

194 def run(self, exposure, sources, expId=0, matches=None):

195 """Measure the PSF.

196

197 Parameters

198 ----------

199 exposure : `lsst.afw.image.Exposure`

200 Exposure to process; measured PSF will be added.

201 sources : `Unknown`

202 Measured sources on exposure; flag fields will be set marking

203 stars chosen by the star selector and the PSF determiner if a schema

204 was passed to the task constructor.

205 expId : `int`, optional

206 Exposure id used for generating random seed.

207 matches : `list`, optional

208 A list of ``lsst.afw.table.ReferenceMatch`` objects

209 (i.e. of ``lsst.afw.table.Match`` with @c first being

210 of type ``lsst.afw.table.SimpleRecord`` and @c second

211 type lsst.afw.table.SourceRecord --- the reference object and detected

212 object respectively) as returned by @em e.g. the AstrometryTask.

213 Used by star selectors that choose to refer to an external catalog.

214

215 Returns

216 -------

217 measurement : `lsst.pipe.base.Struct`

218 PSF measurement as a struct with attributes:

219

220 ``psf``

221 The measured PSF (also set in the input exposure).

222 ``cellSet``

223 An `lsst.afw.math.SpatialCellSet` containing the PSF candidates

224 as returned by the psf determiner.

225

226 Raises

227 ------

228 NonfinitePsfShapeError

229 If the new PSF has NaN or Inf width.

230 """

231 self.log.info("Measuring PSF")

232

233 import lsstDebug

234 display = lsstDebug.Info(__name__).display

235 displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells

236 displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)

237 displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates # show mosaic of candidates

238 displayResiduals = lsstDebug.Info(__name__).displayResiduals # show residuals

239 showBadCandidates = lsstDebug.Info(__name__).showBadCandidates # include bad candidates

240 normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals # normalise residuals by object peak

241

242 #

243 # Run star selector

244 #

245 stars = self.starSelector.run(sourceCat=sources, matches=matches, exposure=exposure)

246 selectionResult = self.makePsfCandidates.run(stars.sourceCat, exposure=exposure)

247 self.log.info("PSF star selector found %d candidates", len(selectionResult.psfCandidates))

248 reserveResult = self.reserve.run(selectionResult.goodStarCat, expId=expId)

249 # Make list of psf candidates to send to the determiner (omitting those marked as reserved)

250 psfDeterminerList = [cand for cand, use

251 in zip(selectionResult.psfCandidates, reserveResult.use) if use]

252

253 if selectionResult.psfCandidates and self.candidateKey is not None:

254 for cand in selectionResult.psfCandidates:

255 source = cand.getSource()

256 source.set(self.candidateKey, True)

257

258 self.log.info("Sending %d candidates to PSF determiner", len(psfDeterminerList))

259

260 if display:

261 frame = 1

262 if displayExposure:

263 disp = afwDisplay.Display(frame=frame)

264 disp.mtv(exposure, title="psf determination")

265 frame += 1

266 #

267 # Determine PSF

268 #

269 psf, cellSet = self.psfDeterminer.determinePsf(exposure, psfDeterminerList, self.metadata,

270 flagKey=self.usedKey)

271 self.log.info("PSF determination using %d/%d stars.",

272 self.metadata.getScalar("numGoodStars"), self.metadata.getScalar("numAvailStars"))

273 if not np.isfinite((psfSize := psf.computeShape(psf.getAveragePosition()).getDeterminantRadius())):

274 raise NonfinitePsfShapeError(psf_size=psfSize)

275 else:

276 self.log.info("Fitted PSF size: %f pixels", psfSize)

277

278 exposure.setPsf(psf)

279

280 if display:

281 frame = display

282 if displayExposure:

283 disp = afwDisplay.Display(frame=frame)

284 showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=frame)

285 frame += 1

286

287 if displayPsfCandidates: # Show a mosaic of PSF candidates

288 plotPsfCandidates(cellSet, showBadCandidates=showBadCandidates, frame=frame)

289 frame += 1

290

291 if displayResiduals:

292 frame = plotResiduals(exposure, cellSet,

293 showBadCandidates=showBadCandidates,

294 normalizeResiduals=normalizeResiduals,

295 frame=frame)

296 if displayPsfMosaic:

297 disp = afwDisplay.Display(frame=frame)

298 maUtils.showPsfMosaic(exposure, psf, display=disp, showFwhm=True)

299 disp.scale("linear", 0, 1)

300 frame += 1

301

302 return pipeBase.Struct(

303 psf=psf,

304 cellSet=cellSet,

305 )

306

307 @property

308 def usesMatches(self):

309 """Return True if this task makes use of the "matches" argument to the run method"""

310 return self.starSelector.usesMatches

311

312#

313# Debug code

314#

315

316

317def showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=1):

318 disp = afwDisplay.Display(frame=frame)

319 maUtils.showPsfSpatialCells(exposure, cellSet,

320 symb="o", ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,

321 size=4, display=disp)

322 for cell in cellSet.getCellList():

323 for cand in cell.begin(not showBadCandidates): # maybe include bad candidates

324 status = cand.getStatus()

325 disp.dot('+', *cand.getSource().getCentroid(),

326 ctype=afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else

327 afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else afwDisplay.RED)

328

329

330def plotPsfCandidates(cellSet, showBadCandidates=False, frame=1):

331 stamps = []

332 for cell in cellSet.getCellList():

333 for cand in cell.begin(not showBadCandidates): # maybe include bad candidates

334 try:

335 im = cand.getMaskedImage()

336

337 chi2 = cand.getChi2()

338 if chi2 < 1e100:

339 chi2 = "%.1f" % chi2

340 else:

341 chi2 = float("nan")

342

343 stamps.append((im, "%d%s" %

344 (maUtils.splitId(cand.getSource().getId(), True)["objId"], chi2),

345 cand.getStatus()))

346 except Exception:

347 continue

348

349 mos = afwDisplay.utils.Mosaic()

350 disp = afwDisplay.Display(frame=frame)

351 for im, label, status in stamps:

352 im = type(im)(im, True)

353 try:

354 im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()

355 except NotImplementedError:

356 pass

357

358 mos.append(im, label,

359 afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else

360 afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else afwDisplay.RED)

361

362 if mos.images:

363 disp.mtv(mos.makeMosaic(), title="Psf Candidates")

364

365

366def plotResiduals(exposure, cellSet, showBadCandidates=False, normalizeResiduals=True, frame=2):

367 psf = exposure.getPsf()

368 disp = afwDisplay.Display(frame=frame)

369 while True:

370 try:

371 maUtils.showPsfCandidates(exposure, cellSet, psf=psf, display=disp,

372 normalize=normalizeResiduals,

373 showBadCandidates=showBadCandidates)

374 frame += 1

375 maUtils.showPsfCandidates(exposure, cellSet, psf=psf, display=disp,

376 normalize=normalizeResiduals,

377 showBadCandidates=showBadCandidates,

378 variance=True)

379 frame += 1

380 except Exception:

381 if not showBadCandidates:

382 showBadCandidates = True

383 continue

384 break

385

386 return frame

Coverage for python / lsst / pipe / tasks / measurePsf.py: 16%

154 statements