LSST Applications 27.0.0,g0265f82a02+469cd937ee,g02d81e74bb+21ad69e7e1,g1470d8bcf6+cbe83ee85a,g2079a07aa2+e67c6346a6,g212a7c68fe+04a9158687,g2305ad1205+94392ce272,g295015adf3+81dd352a9d,g2bbee38e9b+469cd937ee,g337abbeb29+469cd937ee,g3939d97d7f+72a9f7b576,g487adcacf7+71499e7cba,g50ff169b8f+5929b3527e,g52b1c1532d+a6fc98d2e7,g591dd9f2cf+df404f777f,g5a732f18d5+be83d3ecdb,g64a986408d+21ad69e7e1,g858d7b2824+21ad69e7e1,g8a8a8dda67+a6fc98d2e7,g99cad8db69+f62e5b0af5,g9ddcbc5298+d4bad12328,ga1e77700b3+9c366c4306,ga8c6da7877+71e4819109,gb0e22166c9+25ba2f69a1,gb6a65358fc+469cd937ee,gbb8dafda3b+69d3c0e320,gc07e1c2157+a98bf949bb,gc120e1dc64+615ec43309,gc28159a63d+469cd937ee,gcf0d15dbbd+72a9f7b576,gdaeeff99f8+a38ce5ea23,ge6526c86ff+3a7c1ac5f1,ge79ae78c31+469cd937ee,gee10cc3b42+a6fc98d2e7,gf1cff7945b+21ad69e7e1,gfbcc870c63+9a11dc8c8f
LSST Data Management Base Package
Loading...
Searching...
No Matches
characterizeImage.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__all__ = ["CharacterizeImageConfig", "CharacterizeImageTask"]
23
24import numpy as np
25
26from lsstDebug import getDebugFrame
27import lsst.afw.table as afwTable
28import lsst.pex.config as pexConfig
29import lsst.pipe.base as pipeBase
30import lsst.daf.base as dafBase
31import lsst.pipe.base.connectionTypes as cT
32from lsst.afw.math import BackgroundList
33from lsst.afw.table import SourceTable
34from lsst.meas.algorithms import (
35 SubtractBackgroundTask,
36 SourceDetectionTask,
37 MeasureApCorrTask,
38 MeasureApCorrError,
39 MaskStreaksTask,
40)
41from lsst.meas.algorithms.installGaussianPsf import InstallGaussianPsfTask
42from lsst.meas.astrom import displayAstrometry
43from lsst.meas.base import (
44 SingleFrameMeasurementTask,
45 ApplyApCorrTask,
46 CatalogCalculationTask,
47 IdGenerator,
48 DetectorVisitIdGeneratorConfig,
49)
50from lsst.meas.deblender import SourceDeblendTask
51import lsst.meas.extensions.shapeHSM # noqa: F401 needed for default shape plugin
52from .measurePsf import MeasurePsfTask
53from .repair import RepairTask
54from .computeExposureSummaryStats import ComputeExposureSummaryStatsTask
55from lsst.pex.exceptions import LengthError
56from lsst.utils.timer import timeMethod
57
58
59class CharacterizeImageConnections(pipeBase.PipelineTaskConnections,
60 dimensions=("instrument", "visit", "detector")):
61 exposure = cT.Input(
62 doc="Input exposure data",
63 name="postISRCCD",
64 storageClass="Exposure",
65 dimensions=["instrument", "exposure", "detector"],
66 )
67 characterized = cT.Output(
68 doc="Output characterized data.",
69 name="icExp",
70 storageClass="ExposureF",
71 dimensions=["instrument", "visit", "detector"],
72 )
73 sourceCat = cT.Output(
74 doc="Output source catalog.",
75 name="icSrc",
76 storageClass="SourceCatalog",
77 dimensions=["instrument", "visit", "detector"],
78 )
79 backgroundModel = cT.Output(
80 doc="Output background model.",
81 name="icExpBackground",
82 storageClass="Background",
83 dimensions=["instrument", "visit", "detector"],
84 )
85 outputSchema = cT.InitOutput(
86 doc="Schema of the catalog produced by CharacterizeImage",
87 name="icSrc_schema",
88 storageClass="SourceCatalog",
89 )
90
91 def adjustQuantum(self, inputs, outputs, label, dataId):
92 # Docstring inherited from PipelineTaskConnections
93 try:
94 return super().adjustQuantum(inputs, outputs, label, dataId)
95 except pipeBase.ScalarError as err:
96 raise pipeBase.ScalarError(
97 "CharacterizeImageTask can at present only be run on visits that are associated with "
98 "exactly one exposure. Either this is not a valid exposure for this pipeline, or the "
99 "snap-combination step you probably want hasn't been configured to run between ISR and "
100 "this task (as of this writing, that would be because it hasn't been implemented yet)."
101 ) from err
102
103
104class CharacterizeImageConfig(pipeBase.PipelineTaskConfig,
105 pipelineConnections=CharacterizeImageConnections):
106 """Config for CharacterizeImageTask."""
107
108 doMeasurePsf = pexConfig.Field(
109 dtype=bool,
110 default=True,
111 doc="Measure PSF? If False then for all subsequent operations use either existing PSF "
112 "model when present, or install simple PSF model when not (see installSimplePsf "
113 "config options)"
114 )
115 doWrite = pexConfig.Field(
116 dtype=bool,
117 default=True,
118 doc="Persist results?",
119 )
120 doWriteExposure = pexConfig.Field(
121 dtype=bool,
122 default=True,
123 doc="Write icExp and icExpBackground in addition to icSrc? Ignored if doWrite False.",
124 )
125 psfIterations = pexConfig.RangeField(
126 dtype=int,
127 default=2,
128 min=1,
129 doc="Number of iterations of detect sources, measure sources, "
130 "estimate PSF. If useSimplePsf is True then 2 should be plenty; "
131 "otherwise more may be wanted.",
132 )
133 background = pexConfig.ConfigurableField(
134 target=SubtractBackgroundTask,
135 doc="Configuration for initial background estimation",
136 )
137 detection = pexConfig.ConfigurableField(
138 target=SourceDetectionTask,
139 doc="Detect sources"
140 )
141 doDeblend = pexConfig.Field(
142 dtype=bool,
143 default=True,
144 doc="Run deblender input exposure"
145 )
146 deblend = pexConfig.ConfigurableField(
147 target=SourceDeblendTask,
148 doc="Split blended source into their components"
149 )
150 measurement = pexConfig.ConfigurableField(
151 target=SingleFrameMeasurementTask,
152 doc="Measure sources"
153 )
154 doApCorr = pexConfig.Field(
155 dtype=bool,
156 default=True,
157 doc="Run subtasks to measure and apply aperture corrections"
158 )
159 measureApCorr = pexConfig.ConfigurableField(
160 target=MeasureApCorrTask,
161 doc="Subtask to measure aperture corrections"
162 )
163 applyApCorr = pexConfig.ConfigurableField(
164 target=ApplyApCorrTask,
165 doc="Subtask to apply aperture corrections"
166 )
167 # If doApCorr is False, and the exposure does not have apcorrections already applied, the
168 # active plugins in catalogCalculation almost certainly should not contain the characterization plugin
169 catalogCalculation = pexConfig.ConfigurableField(
170 target=CatalogCalculationTask,
171 doc="Subtask to run catalogCalculation plugins on catalog"
172 )
173 doComputeSummaryStats = pexConfig.Field(
174 dtype=bool,
175 default=True,
176 doc="Run subtask to measure exposure summary statistics",
177 deprecated=("This subtask has been moved to CalibrateTask "
178 "with DM-30701.")
179 )
180 computeSummaryStats = pexConfig.ConfigurableField(
181 target=ComputeExposureSummaryStatsTask,
182 doc="Subtask to run computeSummaryStats on exposure",
183 deprecated=("This subtask has been moved to CalibrateTask "
184 "with DM-30701.")
185 )
186 useSimplePsf = pexConfig.Field(
187 dtype=bool,
188 default=True,
189 doc="Replace the existing PSF model with a simplified version that has the same sigma "
190 "at the start of each PSF determination iteration? Doing so makes PSF determination "
191 "converge more robustly and quickly.",
192 )
193 installSimplePsf = pexConfig.ConfigurableField(
194 target=InstallGaussianPsfTask,
195 doc="Install a simple PSF model",
196 )
197 measurePsf = pexConfig.ConfigurableField(
198 target=MeasurePsfTask,
199 doc="Measure PSF",
200 )
201 repair = pexConfig.ConfigurableField(
202 target=RepairTask,
203 doc="Remove cosmic rays",
204 )
205 requireCrForPsf = pexConfig.Field(
206 dtype=bool,
207 default=True,
208 doc="Require cosmic ray detection and masking to run successfully before measuring the PSF."
209 )
210 checkUnitsParseStrict = pexConfig.Field(
211 doc="Strictness of Astropy unit compatibility check, can be 'raise', 'warn' or 'silent'",
212 dtype=str,
213 default="raise",
214 )
215 doMaskStreaks = pexConfig.Field(
216 doc="Mask streaks",
217 default=True,
218 dtype=bool,
219 )
220 maskStreaks = pexConfig.ConfigurableField(
221 target=MaskStreaksTask,
222 doc="Subtask for masking streaks. Only used if doMaskStreaks is True. "
223 "Adds a mask plane to an exposure, with the mask plane name set by streakMaskName.",
224 )
225 idGenerator = DetectorVisitIdGeneratorConfig.make_field()
226
227 def setDefaults(self):
228 super().setDefaults()
229 # just detect bright stars; includeThresholdMultipler=10 seems large,
230 # but these are the values we have been using
231 self.detection.thresholdValue = 5.0
232 self.detection.includeThresholdMultiplier = 10.0
233 # do not deblend, as it makes a mess
234 self.doDeblend = False
235 # measure and apply aperture correction; note: measuring and applying aperture
236 # correction are disabled until the final measurement, after PSF is measured
237 self.doApCorr = True
238 # During characterization, we don't have full source measurement information,
239 # so must do the aperture correction with only psf stars, combined with the
240 # default signal-to-noise cuts in MeasureApCorrTask.
241 selector = self.measureApCorr.sourceSelector["science"]
242 selector.doUnresolved = False
243 selector.flags.good = ["calib_psf_used"]
244 selector.flags.bad = []
245
246 # minimal set of measurements needed to determine PSF
247 self.measurement.plugins.names = [
248 "base_PixelFlags",
249 "base_SdssCentroid",
250 "ext_shapeHSM_HsmSourceMoments",
251 "base_GaussianFlux",
252 "base_PsfFlux",
253 "base_CircularApertureFlux",
254 "base_ClassificationSizeExtendedness",
255 ]
256 self.measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
257
258 def validate(self):
259 if self.doApCorr and not self.measurePsf:
260 raise RuntimeError("Must measure PSF to measure aperture correction, "
261 "because flags determined by PSF measurement are used to identify "
262 "sources used to measure aperture correction")
263
264
265class CharacterizeImageTask(pipeBase.PipelineTask):
266 """Measure bright sources and use this to estimate background and PSF of
267 an exposure.
268
269 Given an exposure with defects repaired (masked and interpolated over,
270 e.g. as output by `~lsst.ip.isr.IsrTask`):
271 - detect and measure bright sources
272 - repair cosmic rays
273 - detect and mask streaks
274 - measure and subtract background
275 - measure PSF
276
277 Parameters
278 ----------
279 schema : `lsst.afw.table.Schema`, optional
280 Initial schema for icSrc catalog.
281 **kwargs
282 Additional keyword arguments.
283
284 Notes
285 -----
286 Debugging:
287 CharacterizeImageTask has a debug dictionary with the following keys:
288
289 frame
290 int: if specified, the frame of first debug image displayed (defaults to 1)
291 repair_iter
292 bool; if True display image after each repair in the measure PSF loop
293 background_iter
294 bool; if True display image after each background subtraction in the measure PSF loop
295 measure_iter
296 bool; if True display image and sources at the end of each iteration of the measure PSF loop
297 See `~lsst.meas.astrom.displayAstrometry` for the meaning of the various symbols.
298 psf
299 bool; if True display image and sources after PSF is measured;
300 this will be identical to the final image displayed by measure_iter if measure_iter is true
301 repair
302 bool; if True display image and sources after final repair
303 measure
304 bool; if True display image and sources after final measurement
305 """
306
307 ConfigClass = CharacterizeImageConfig
308 _DefaultName = "characterizeImage"
309
310 def __init__(self, schema=None, **kwargs):
311 super().__init__(**kwargs)
312
313 if schema is None:
314 schema = SourceTable.makeMinimalSchema()
315 self.schema = schema
316 self.makeSubtask("background")
317 self.makeSubtask("installSimplePsf")
318 self.makeSubtask("repair")
319 if self.config.doMaskStreaks:
320 self.makeSubtask("maskStreaks")
321 self.makeSubtask("measurePsf", schema=self.schema)
323 self.makeSubtask('detection', schema=self.schema)
324 if self.config.doDeblend:
325 self.makeSubtask("deblend", schema=self.schema)
326 self.makeSubtask('measurement', schema=self.schema, algMetadata=self.algMetadata)
327 if self.config.doApCorr:
328 self.makeSubtask('measureApCorr', schema=self.schema)
329 self.makeSubtask('applyApCorr', schema=self.schema)
330 self.makeSubtask('catalogCalculation', schema=self.schema)
331 self._initialFrame = getDebugFrame(self._display, "frame") or 1
333 self.schema.checkUnits(parse_strict=self.config.checkUnitsParseStrict)
335
336 def runQuantum(self, butlerQC, inputRefs, outputRefs):
337 inputs = butlerQC.get(inputRefs)
338 if 'idGenerator' not in inputs.keys():
339 inputs['idGenerator'] = self.config.idGenerator.apply(butlerQC.quantum.dataId)
340 outputs = self.run(**inputs)
341 butlerQC.put(outputs, outputRefs)
342
343 @timeMethod
344 def run(self, exposure, background=None, idGenerator=None):
345 """Characterize a science image.
346
347 Peforms the following operations:
348 - Iterate the following config.psfIterations times, or once if config.doMeasurePsf false:
349 - detect and measure sources and estimate PSF (see detectMeasureAndEstimatePsf for details)
350 - interpolate over cosmic rays
351 - perform final measurement
352
353 Parameters
354 ----------
355 exposure : `lsst.afw.image.ExposureF`
356 Exposure to characterize.
357 background : `lsst.afw.math.BackgroundList`, optional
358 Initial model of background already subtracted from exposure.
359 idGenerator : `lsst.meas.base.IdGenerator`, optional
360 Object that generates source IDs and provides RNG seeds.
361
362 Returns
363 -------
364 result : `lsst.pipe.base.Struct`
365 Results as a struct with attributes:
366
367 ``exposure``
368 Characterized exposure (`lsst.afw.image.ExposureF`).
369 ``sourceCat``
370 Detected sources (`lsst.afw.table.SourceCatalog`).
371 ``background``
372 Model of subtracted background (`lsst.afw.math.BackgroundList`).
373 ``psfCellSet``
374 Spatial cells of PSF candidates (`lsst.afw.math.SpatialCellSet`).
375 ``characterized``
376 Another reference to ``exposure`` for compatibility.
377 ``backgroundModel``
378 Another reference to ``background`` for compatibility.
379
380 Raises
381 ------
382 RuntimeError
383 Raised if PSF sigma is NaN.
384 """
385 self._frame = self._initialFrame # reset debug display frame
386
387 if not self.config.doMeasurePsf and not exposure.hasPsf():
388 self.log.info("CharacterizeImageTask initialized with 'simple' PSF.")
389 self.installSimplePsf.run(exposure=exposure)
390
391 if idGenerator is None:
392 idGenerator = IdGenerator()
393
394 # subtract an initial estimate of background level
395 background = self.background.run(exposure).background
396
397 psfIterations = self.config.psfIterations if self.config.doMeasurePsf else 1
398 for i in range(psfIterations):
399 dmeRes = self.detectMeasureAndEstimatePsf(
400 exposure=exposure,
401 idGenerator=idGenerator,
402 background=background,
403 )
404
405 psf = dmeRes.exposure.getPsf()
406 # Just need a rough estimate; average positions are fine
407 psfAvgPos = psf.getAveragePosition()
408 psfSigma = psf.computeShape(psfAvgPos).getDeterminantRadius()
409 psfDimensions = psf.computeImage(psfAvgPos).getDimensions()
410 medBackground = np.median(dmeRes.background.getImage().getArray())
411 self.log.info("iter %s; PSF sigma=%0.4f, dimensions=%s; median background=%0.2f",
412 i + 1, psfSigma, psfDimensions, medBackground)
413 if np.isnan(psfSigma):
414 raise RuntimeError("PSF sigma is NaN, cannot continue PSF determination.")
415
416 self.display("psf", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
417
418 # perform final repair with final PSF
419 self.repair.run(exposure=dmeRes.exposure)
420 self.display("repair", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
421
422 # mask streaks
423 if self.config.doMaskStreaks:
424 _ = self.maskStreaks.run(dmeRes.exposure)
425
426 # perform final measurement with final PSF, including measuring and applying aperture correction,
427 # if wanted
428 self.measurement.run(measCat=dmeRes.sourceCat, exposure=dmeRes.exposure,
429 exposureId=idGenerator.catalog_id)
430 if self.config.doApCorr:
431 try:
432 apCorrMap = self.measureApCorr.run(
433 exposure=dmeRes.exposure,
434 catalog=dmeRes.sourceCat,
435 ).apCorrMap
436 except MeasureApCorrError:
437 # We have failed to get a valid aperture correction map.
438 # Proceed with processing, and image will be filtered
439 # downstream.
440 dmeRes.exposure.info.setApCorrMap(None)
441 else:
442 dmeRes.exposure.info.setApCorrMap(apCorrMap)
443 self.applyApCorr.run(catalog=dmeRes.sourceCat, apCorrMap=exposure.getInfo().getApCorrMap())
444
445 self.catalogCalculation.run(dmeRes.sourceCat)
446
447 self.display("measure", exposure=dmeRes.exposure, sourceCat=dmeRes.sourceCat)
448
449 return pipeBase.Struct(
450 exposure=dmeRes.exposure,
451 sourceCat=dmeRes.sourceCat,
452 background=dmeRes.background,
453 psfCellSet=dmeRes.psfCellSet,
454
455 characterized=dmeRes.exposure,
456 backgroundModel=dmeRes.background
457 )
458
459 @timeMethod
460 def detectMeasureAndEstimatePsf(self, exposure, idGenerator, background):
461 """Perform one iteration of detect, measure, and estimate PSF.
462
463 Performs the following operations:
464
465 - if config.doMeasurePsf or not exposure.hasPsf():
466
467 - install a simple PSF model (replacing the existing one, if need be)
468
469 - interpolate over cosmic rays with keepCRs=True
470 - estimate background and subtract it from the exposure
471 - detect, deblend and measure sources, and subtract a refined background model;
472 - if config.doMeasurePsf:
473 - measure PSF
474
475 Parameters
476 ----------
477 exposure : `lsst.afw.image.ExposureF`
478 Exposure to characterize.
479 idGenerator : `lsst.meas.base.IdGenerator`
480 Object that generates source IDs and provides RNG seeds.
481 background : `lsst.afw.math.BackgroundList`, optional
482 Initial model of background already subtracted from exposure.
483
484 Returns
485 -------
486 result : `lsst.pipe.base.Struct`
487 Results as a struct with attributes:
488
489 ``exposure``
490 Characterized exposure (`lsst.afw.image.ExposureF`).
491 ``sourceCat``
492 Detected sources (`lsst.afw.table.SourceCatalog`).
493 ``background``
494 Model of subtracted background (`lsst.afw.math.BackgroundList`).
495 ``psfCellSet``
496 Spatial cells of PSF candidates (`lsst.afw.math.SpatialCellSet`).
497
498 Raises
499 ------
500 LengthError
501 Raised if there are too many CR pixels.
502 """
503 # install a simple PSF model, if needed or wanted
504 if not exposure.hasPsf() or (self.config.doMeasurePsf and self.config.useSimplePsf):
505 self.log.info("PSF estimation initialized with 'simple' PSF")
506 self.installSimplePsf.run(exposure=exposure)
507
508 # run repair, but do not interpolate over cosmic rays (do that elsewhere, with the final PSF model)
509 if self.config.requireCrForPsf:
510 self.repair.run(exposure=exposure, keepCRs=True)
511 else:
512 try:
513 self.repair.run(exposure=exposure, keepCRs=True)
514 except LengthError:
515 self.log.warning("Skipping cosmic ray detection: Too many CR pixels (max %0.f)",
516 self.config.repair.cosmicray.nCrPixelMax)
517
518 self.display("repair_iter", exposure=exposure)
519
520 if background is None:
521 background = BackgroundList()
522
523 sourceIdFactory = idGenerator.make_table_id_factory()
524 table = SourceTable.make(self.schema, sourceIdFactory)
525 table.setMetadata(self.algMetadata)
526
527 detRes = self.detection.run(table=table, exposure=exposure, doSmooth=True)
528 sourceCat = detRes.sources
529 if detRes.background:
530 for bg in detRes.background:
531 background.append(bg)
532
533 if self.config.doDeblend:
534 self.deblend.run(exposure=exposure, sources=sourceCat)
535 # We need the output catalog to be contiguous for further processing.
536 if not sourceCat.isContiguous():
537 sourceCat = sourceCat.copy(deep=True)
538
539 self.measurement.run(measCat=sourceCat, exposure=exposure, exposureId=idGenerator.catalog_id)
540
541 measPsfRes = pipeBase.Struct(cellSet=None)
542 if self.config.doMeasurePsf:
543 measPsfRes = self.measurePsf.run(exposure=exposure, sources=sourceCat,
544 expId=idGenerator.catalog_id)
545 self.display("measure_iter", exposure=exposure, sourceCat=sourceCat)
546
547 return pipeBase.Struct(
548 exposure=exposure,
549 sourceCat=sourceCat,
550 background=background,
551 psfCellSet=measPsfRes.cellSet,
552 )
553
554 def display(self, itemName, exposure, sourceCat=None):
555 """Display exposure and sources on next frame (for debugging).
556
557 Parameters
558 ----------
559 itemName : `str`
560 Name of item in ``debugInfo``.
561 exposure : `lsst.afw.image.ExposureF`
562 Exposure to display.
563 sourceCat : `lsst.afw.table.SourceCatalog`, optional
564 Catalog of sources detected on the exposure.
565 """
566 val = getDebugFrame(self._display, itemName)
567 if not val:
568 return
569
570 displayAstrometry(exposure=exposure, sourceCat=sourceCat, frame=self._frame, pause=False)
571 self._frame += 1
Class for storing ordered metadata with comments.
detectMeasureAndEstimatePsf(self, exposure, idGenerator, background)
run(self, exposure, background=None, idGenerator=None)
display(self, itemName, exposure, sourceCat=None)
runQuantum(self, butlerQC, inputRefs, outputRefs)