LSST Applications  21.0.0-172-gfb10e10a+18fedfabac,22.0.0+297cba6710,22.0.0+80564b0ff1,22.0.0+8d77f4f51a,22.0.0+a28f4c53b1,22.0.0+dcf3732eb2,22.0.1-1-g7d6de66+2a20fdde0d,22.0.1-1-g8e32f31+297cba6710,22.0.1-1-geca5380+7fa3b7d9b6,22.0.1-12-g44dc1dc+2a20fdde0d,22.0.1-15-g6a90155+515f58c32b,22.0.1-16-g9282f48+790f5f2caa,22.0.1-2-g92698f7+dcf3732eb2,22.0.1-2-ga9b0f51+7fa3b7d9b6,22.0.1-2-gd1925c9+bf4f0e694f,22.0.1-24-g1ad7a390+a9625a72a8,22.0.1-25-g5bf6245+3ad8ecd50b,22.0.1-25-gb120d7b+8b5510f75f,22.0.1-27-g97737f7+2a20fdde0d,22.0.1-32-gf62ce7b1+aa4237961e,22.0.1-4-g0b3f228+2a20fdde0d,22.0.1-4-g243d05b+871c1b8305,22.0.1-4-g3a563be+32dcf1063f,22.0.1-4-g44f2e3d+9e4ab0f4fa,22.0.1-42-gca6935d93+ba5e5ca3eb,22.0.1-5-g15c806e+85460ae5f3,22.0.1-5-g58711c4+611d128589,22.0.1-5-g75bb458+99c117b92f,22.0.1-6-g1c63a23+7fa3b7d9b6,22.0.1-6-g50866e6+84ff5a128b,22.0.1-6-g8d3140d+720564cf76,22.0.1-6-gd805d02+cc5644f571,22.0.1-8-ge5750ce+85460ae5f3,master-g6e05de7fdc+babf819c66,master-g99da0e417a+8d77f4f51a,w.2021.48
LSST Data Management Base Package
measurePsf.py
Go to the documentation of this file.
1 # This file is part of pipe_tasks.
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 
22 import lsst.afw.display as afwDisplay
23 import lsst.afw.math as afwMath
24 import lsst.meas.algorithms as measAlg
25 import lsst.meas.algorithms.utils as maUtils
26 import lsst.pex.config as pexConfig
27 import lsst.pipe.base as pipeBase
29 from lsst.utils.timer import timeMethod
30 
31 
32 class MeasurePsfConfig(pexConfig.Config):
33  starSelector = measAlg.sourceSelectorRegistry.makeField(
34  "Star selection algorithm",
35  default="objectSize"
36  )
37  makePsfCandidates = pexConfig.ConfigurableField(
38  target=measAlg.MakePsfCandidatesTask,
39  doc="Task to make psf candidates from selected stars.",
40  )
41  psfDeterminer = measAlg.psfDeterminerRegistry.makeField(
42  "PSF Determination algorithm",
43  default="psfex"
44  )
45  reserve = pexConfig.ConfigurableField(
46  target=measAlg.ReserveSourcesTask,
47  doc="Reserve sources from fitting"
48  )
49 
50 
56 
57 
58 class MeasurePsfTask(pipeBase.Task):
59  r"""!
60 @anchor MeasurePsfTask_
61 
62 @brief Measure the PSF
63 
64 @section pipe_tasks_measurePsf_Contents Contents
65 
66  - @ref pipe_tasks_measurePsf_Purpose
67  - @ref pipe_tasks_measurePsf_Initialize
68  - @ref pipe_tasks_measurePsf_IO
69  - @ref pipe_tasks_measurePsf_Config
70  - @ref pipe_tasks_measurePsf_Debug
71  - @ref pipe_tasks_measurePsf_Example
72 
73 @section pipe_tasks_measurePsf_Purpose Description
74 
75 A task that selects stars from a catalog of sources and uses those to measure the PSF.
76 
77 The star selector is a subclass of
78 @ref lsst.meas.algorithms.starSelector.BaseStarSelectorTask "lsst.meas.algorithms.BaseStarSelectorTask"
79 and the PSF determiner is a sublcass of
80 @ref lsst.meas.algorithms.psfDeterminer.BasePsfDeterminerTask "lsst.meas.algorithms.BasePsfDeterminerTask"
81 
82 @warning
83 There is no establised set of configuration parameters for these algorithms, so once you start modifying
84 parameters (as we do in @ref pipe_tasks_measurePsf_Example) your code is no longer portable.
85 
86 @section pipe_tasks_measurePsf_Initialize Task initialisation
87 
88 @copydoc \_\_init\_\_
89 
90 @section pipe_tasks_measurePsf_IO Invoking the Task
91 
92 @copydoc run
93 
94 @section pipe_tasks_measurePsf_Config Configuration parameters
95 
96 See @ref MeasurePsfConfig.
97 
98 @section pipe_tasks_measurePsf_Debug Debug variables
99 
100 The @link lsst.pipe.base.cmdLineTask.CmdLineTask command line task@endlink interface supports a
101 flag @c -d to import @b debug.py from your @c PYTHONPATH; see @ref baseDebug for more about @b debug.py files.
102 
103 <DL>
104  <DT> @c display
105  <DD> If True, display debugging plots
106  <DT> displayExposure
107  <DD> display the Exposure + spatialCells
108  <DT> displayPsfCandidates
109  <DD> show mosaic of candidates
110  <DT> showBadCandidates
111  <DD> Include bad candidates
112  <DT> displayPsfMosaic
113  <DD> show mosaic of reconstructed PSF(xy)
114  <DT> displayResiduals
115  <DD> show residuals
116  <DT> normalizeResiduals
117  <DD> Normalise residuals by object amplitude
118 </DL>
119 
120 Additionally you can enable any debug outputs that your chosen star selector and psf determiner support.
121 
122 @section pipe_tasks_measurePsf_Example A complete example of using MeasurePsfTask
123 
124 This code is in @link measurePsfTask.py@endlink in the examples directory, and can be run as @em e.g.
125 @code
126 examples/measurePsfTask.py --doDisplay
127 @endcode
128 @dontinclude measurePsfTask.py
129 
130 The example also runs SourceDetectionTask and SingleFrameMeasurementTask;
131 see @ref meas_algorithms_measurement_Example for more explanation.
132 
133 Import the tasks (there are some other standard imports; read the file to see them all):
134 
135 @skip SourceDetectionTask
136 @until MeasurePsfTask
137 
138 We need to create the tasks before processing any data as the task constructor
139 can add an extra column to the schema, but first we need an almost-empty
140 Schema:
141 
142 @skipline makeMinimalSchema
143 
144 We can now call the constructors for the tasks we need to find and characterize candidate
145 PSF stars:
146 
147 @skip SourceDetectionTask.ConfigClass
148 @until measureTask
149 
150 Note that we've chosen a minimal set of measurement plugins: we need the
151 outputs of @c base_SdssCentroid, @c base_SdssShape and @c base_CircularApertureFlux as
152 inputs to the PSF measurement algorithm, while @c base_PixelFlags identifies
153 and flags bad sources (e.g. with pixels too close to the edge) so they can be
154 removed later.
155 
156 Now we can create and configure the task that we're interested in:
157 
158 @skip MeasurePsfTask
159 @until measurePsfTask
160 
161 We're now ready to process the data (we could loop over multiple exposures/catalogues using the same
162 task objects). First create the output table:
163 
164 @skipline afwTable
165 
166 And process the image:
167 
168 @skip sources =
169 @until result
170 
171 We can then unpack and use the results:
172 
173 @skip psf
174 @until cellSet
175 
176 If you specified @c --doDisplay you can see the PSF candidates:
177 
178 @skip display
179 @until RED
180 
181 <HR>
182 
183 To investigate the @ref pipe_tasks_measurePsf_Debug, put something like
184 @code{.py}
185  import lsstDebug
186  def DebugInfo(name):
187  di = lsstDebug.getInfo(name) # N.b. lsstDebug.Info(name) would call us recursively
188 
189  if name == "lsst.pipe.tasks.measurePsf" :
190  di.display = True
191  di.displayExposure = False # display the Exposure + spatialCells
192  di.displayPsfCandidates = True # show mosaic of candidates
193  di.displayPsfMosaic = True # show mosaic of reconstructed PSF(xy)
194  di.displayResiduals = True # show residuals
195  di.showBadCandidates = True # Include bad candidates
196  di.normalizeResiduals = False # Normalise residuals by object amplitude
197 
198  return di
199 
200  lsstDebug.Info = DebugInfo
201 @endcode
202 into your debug.py file and run measurePsfTask.py with the @c --debug flag.
203  """
204  ConfigClass = MeasurePsfConfig
205  _DefaultName = "measurePsf"
206 
207  def __init__(self, schema=None, **kwargs):
208  """!Create the detection task. Most arguments are simply passed onto pipe.base.Task.
209 
210  @param schema An lsst::afw::table::Schema used to create the output lsst.afw.table.SourceCatalog
211  @param **kwargs Keyword arguments passed to lsst.pipe.base.task.Task.__init__.
212 
213  If schema is not None, 'calib_psf_candidate' and 'calib_psf_used' fields will be added to
214  identify which stars were employed in the PSF estimation.
215 
216  @note This task can add fields to the schema, so any code calling this task must ensure that
217  these fields are indeed present in the input table.
218  """
219 
220  pipeBase.Task.__init__(self, **kwargs)
221  if schema is not None:
222  self.candidateKeycandidateKey = schema.addField(
223  "calib_psf_candidate", type="Flag",
224  doc=("Flag set if the source was a candidate for PSF determination, "
225  "as determined by the star selector.")
226  )
227  self.usedKeyusedKey = schema.addField(
228  "calib_psf_used", type="Flag",
229  doc=("Flag set if the source was actually used for PSF determination, "
230  "as determined by the '%s' PSF determiner.") % self.config.psfDeterminer.name
231  )
232  else:
233  self.candidateKeycandidateKey = None
234  self.usedKeyusedKey = None
235  self.makeSubtask("starSelector")
236  self.makeSubtask("makePsfCandidates")
237  self.makeSubtask("psfDeterminer", schema=schema)
238  self.makeSubtask("reserve", columnName="calib_psf", schema=schema,
239  doc="set if source was reserved from PSF determination")
240 
241  @timeMethod
242  def run(self, exposure, sources, expId=0, matches=None):
243  """!Measure the PSF
244 
245  @param[in,out] exposure Exposure to process; measured PSF will be added.
246  @param[in,out] sources Measured sources on exposure; flag fields will be set marking
247  stars chosen by the star selector and the PSF determiner if a schema
248  was passed to the task constructor.
249  @param[in] expId Exposure id used for generating random seed.
250  @param[in] matches A list of lsst.afw.table.ReferenceMatch objects
251  (@em i.e. of lsst.afw.table.Match
252  with @c first being of type lsst.afw.table.SimpleRecord and @c second
253  type lsst.afw.table.SourceRecord --- the reference object and detected
254  object respectively) as returned by @em e.g. the AstrometryTask.
255  Used by star selectors that choose to refer to an external catalog.
256 
257  @return a pipe.base.Struct with fields:
258  - psf: The measured PSF (also set in the input exposure)
259  - cellSet: an lsst.afw.math.SpatialCellSet containing the PSF candidates
260  as returned by the psf determiner.
261  """
262  self.log.info("Measuring PSF")
263 
264  import lsstDebug
265  display = lsstDebug.Info(__name__).display
266  displayExposure = lsstDebug.Info(__name__).displayExposure # display the Exposure + spatialCells
267  displayPsfMosaic = lsstDebug.Info(__name__).displayPsfMosaic # show mosaic of reconstructed PSF(x,y)
268  displayPsfCandidates = lsstDebug.Info(__name__).displayPsfCandidates # show mosaic of candidates
269  displayResiduals = lsstDebug.Info(__name__).displayResiduals # show residuals
270  showBadCandidates = lsstDebug.Info(__name__).showBadCandidates # include bad candidates
271  normalizeResiduals = lsstDebug.Info(__name__).normalizeResiduals # normalise residuals by object peak
272 
273  #
274  # Run star selector
275  #
276  stars = self.starSelector.run(sourceCat=sources, matches=matches, exposure=exposure)
277  selectionResult = self.makePsfCandidates.run(stars.sourceCat, exposure=exposure)
278  self.log.info("PSF star selector found %d candidates", len(selectionResult.psfCandidates))
279  reserveResult = self.reserve.run(selectionResult.goodStarCat, expId=expId)
280  # Make list of psf candidates to send to the determiner (omitting those marked as reserved)
281  psfDeterminerList = [cand for cand, use
282  in zip(selectionResult.psfCandidates, reserveResult.use) if use]
283 
284  if selectionResult.psfCandidates and self.candidateKeycandidateKey is not None:
285  for cand in selectionResult.psfCandidates:
286  source = cand.getSource()
287  source.set(self.candidateKeycandidateKey, True)
288 
289  self.log.info("Sending %d candidates to PSF determiner", len(psfDeterminerList))
290 
291  if display:
292  frame = 1
293  if displayExposure:
294  disp = afwDisplay.Display(frame=frame)
295  disp.mtv(exposure, title="psf determination")
296  frame += 1
297  #
298  # Determine PSF
299  #
300  psf, cellSet = self.psfDeterminer.determinePsf(exposure, psfDeterminerList, self.metadata,
301  flagKey=self.usedKeyusedKey)
302  self.log.info("PSF determination using %d/%d stars.",
303  self.metadata.getScalar("numGoodStars"), self.metadata.getScalar("numAvailStars"))
304 
305  exposure.setPsf(psf)
306 
307  if display:
308  frame = display
309  if displayExposure:
310  disp = afwDisplay.Display(frame=frame)
311  showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=frame)
312  frame += 1
313 
314  if displayPsfCandidates: # Show a mosaic of PSF candidates
315  plotPsfCandidates(cellSet, showBadCandidates=showBadCandidates, frame=frame)
316  frame += 1
317 
318  if displayResiduals:
319  frame = plotResiduals(exposure, cellSet,
320  showBadCandidates=showBadCandidates,
321  normalizeResiduals=normalizeResiduals,
322  frame=frame)
323  if displayPsfMosaic:
324  disp = afwDisplay.Display(frame=frame)
325  maUtils.showPsfMosaic(exposure, psf, display=disp, showFwhm=True)
326  disp.scale("linear", 0, 1)
327  frame += 1
328 
329  return pipeBase.Struct(
330  psf=psf,
331  cellSet=cellSet,
332  )
333 
334  @property
335  def usesMatches(self):
336  """Return True if this task makes use of the "matches" argument to the run method"""
337  return self.starSelector.usesMatches
338 
339 #
340 # Debug code
341 #
342 
343 
344 def showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=1):
345  disp = afwDisplay.Display(frame=frame)
346  maUtils.showPsfSpatialCells(exposure, cellSet,
347  symb="o", ctype=afwDisplay.CYAN, ctypeUnused=afwDisplay.YELLOW,
348  size=4, display=disp)
349  for cell in cellSet.getCellList():
350  for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
351  status = cand.getStatus()
352  disp.dot('+', *cand.getSource().getCentroid(),
353  ctype=afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
354  afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else afwDisplay.RED)
355 
356 
357 def plotPsfCandidates(cellSet, showBadCandidates=False, frame=1):
358  stamps = []
359  for cell in cellSet.getCellList():
360  for cand in cell.begin(not showBadCandidates): # maybe include bad candidates
361  try:
362  im = cand.getMaskedImage()
363 
364  chi2 = cand.getChi2()
365  if chi2 < 1e100:
366  chi2 = "%.1f" % chi2
367  else:
368  chi2 = float("nan")
369 
370  stamps.append((im, "%d%s" %
371  (maUtils.splitId(cand.getSource().getId(), True)["objId"], chi2),
372  cand.getStatus()))
373  except Exception:
374  continue
375 
376  mos = afwDisplay.utils.Mosaic()
377  disp = afwDisplay.Display(frame=frame)
378  for im, label, status in stamps:
379  im = type(im)(im, True)
380  try:
381  im /= afwMath.makeStatistics(im, afwMath.MAX).getValue()
382  except NotImplementedError:
383  pass
384 
385  mos.append(im, label,
386  afwDisplay.GREEN if status == afwMath.SpatialCellCandidate.GOOD else
387  afwDisplay.YELLOW if status == afwMath.SpatialCellCandidate.UNKNOWN else afwDisplay.RED)
388 
389  if mos.images:
390  disp.mtv(mos.makeMosaic(), title="Psf Candidates")
391 
392 
393 def plotResiduals(exposure, cellSet, showBadCandidates=False, normalizeResiduals=True, frame=2):
394  psf = exposure.getPsf()
395  disp = afwDisplay.Display(frame=frame)
396  while True:
397  try:
398  maUtils.showPsfCandidates(exposure, cellSet, psf=psf, display=disp,
399  normalize=normalizeResiduals,
400  showBadCandidates=showBadCandidates)
401  frame += 1
402  maUtils.showPsfCandidates(exposure, cellSet, psf=psf, display=disp,
403  normalize=normalizeResiduals,
404  showBadCandidates=showBadCandidates,
405  variance=True)
406  frame += 1
407  except Exception:
408  if not showBadCandidates:
409  showBadCandidates = True
410  continue
411  break
412 
413  return frame
table::Key< int > type
Definition: Detector.cc:163
def run(self, exposure, sources, expId=0, matches=None)
Measure the PSF.
Definition: measurePsf.py:242
def __init__(self, schema=None, **kwargs)
Create the detection task.
Definition: measurePsf.py:207
Statistics makeStatistics(lsst::afw::image::Image< Pixel > const &img, lsst::afw::image::Mask< image::MaskPixel > const &msk, int const flags, StatisticsControl const &sctrl=StatisticsControl())
Handle a watered-down front-end to the constructor (no variance)
Definition: Statistics.h:359
Fit spatial kernel using approximate fluxes for candidates, and solving a linear system of equations.
def plotPsfCandidates(cellSet, showBadCandidates=False, frame=1)
Definition: measurePsf.py:357
def plotResiduals(exposure, cellSet, showBadCandidates=False, normalizeResiduals=True, frame=2)
Definition: measurePsf.py:393
def showPsfSpatialCells(exposure, cellSet, showBadCandidates, frame=1)
Definition: measurePsf.py:344