Coverage for python/lsst/pipe/tasks/calibrateImage.py: 13%
691 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-13 02:08 -0800
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/>.
22__all__ = ["CalibrateImageTask", "CalibrateImageConfig", "NoPsfStarsToStarsMatchError",
23 "AllCentroidsFlaggedError"]
25import numpy as np
26import requests
27import os
29from lsst.afw.geom import SpanSet
30import lsst.afw.table as afwTable
31import lsst.afw.image as afwImage
32import lsst.afw.math as afwMath
33from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache
34import lsst.meas.algorithms
35import lsst.meas.algorithms.installGaussianPsf
36import lsst.meas.algorithms.measureApCorr
37import lsst.meas.algorithms.setPrimaryFlags
38from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask
39import lsst.meas.base
40import lsst.meas.astrom
41import lsst.meas.deblender
42import lsst.meas.extensions.shapeHSM
43from lsst.obs.base import createInitialSkyWcs
44import lsst.pex.config as pexConfig
45import lsst.pipe.base as pipeBase
46from lsst.pipe.base import connectionTypes
47from lsst.utils.timer import timeMethod
49from . import measurePsf, repair, photoCal, computeExposureSummaryStats, snapCombine, diffractionSpikeMask
52class AllCentroidsFlaggedError(pipeBase.AlgorithmError):
53 """Raised when there are no valid centroids during psf fitting.
54 """
55 def __init__(self, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size):
56 msg = (f"All source centroids (out of {n_sources}) flagged during PSF fitting. "
57 "Original image PSF is likely unuseable; best-fit PSF shape parameters: "
58 f"Ixx={psf_shape_ixx}, Iyy={psf_shape_iyy}, Ixy={psf_shape_ixy}, size={psf_size}"
59 )
60 super().__init__(msg)
61 self.n_sources = n_sources
62 self.psf_shape_ixx = psf_shape_ixx
63 self.psf_shape_iyy = psf_shape_iyy
64 self.psf_shape_ixy = psf_shape_ixy
65 self.psf_size = psf_size
67 @property
68 def metadata(self):
69 return {"n_sources": self.n_sources,
70 "psf_shape_ixx": self.psf_shape_ixx,
71 "psf_shape_iyy": self.psf_shape_iyy,
72 "psf_shape_ixy": self.psf_shape_ixy,
73 "psf_size": self.psf_size
74 }
77class NoPsfStarsToStarsMatchError(pipeBase.AlgorithmError):
78 """Raised when there are no matches between the psf_stars and stars
79 catalogs.
80 """
81 def __init__(self, *, n_psf_stars, n_stars):
82 msg = (f"No psf stars out of {n_psf_stars} matched {n_stars} calib stars."
83 " Downstream processes probably won't have useful {calib_type} stars in this case."
84 " Is `star_selector` too strict or is this a bad image?")
85 super().__init__(msg)
86 self.n_psf_stars = n_psf_stars
87 self.n_stars = n_stars
89 @property
90 def metadata(self):
91 return {"n_psf_stars": self.n_psf_stars,
92 "n_stars": self.n_stars
93 }
96class CalibrateImageConnections(pipeBase.PipelineTaskConnections,
97 dimensions=("instrument", "visit", "detector")):
99 astrometry_ref_cat = connectionTypes.PrerequisiteInput(
100 doc="Reference catalog to use for astrometric calibration.",
101 name="gaia_dr3_20230707",
102 storageClass="SimpleCatalog",
103 dimensions=("skypix",),
104 deferLoad=True,
105 multiple=True,
106 )
107 photometry_ref_cat = connectionTypes.PrerequisiteInput(
108 doc="Reference catalog to use for photometric calibration.",
109 name="ps1_pv3_3pi_20170110",
110 storageClass="SimpleCatalog",
111 dimensions=("skypix",),
112 deferLoad=True,
113 multiple=True
114 )
115 camera_model = connectionTypes.PrerequisiteInput(
116 doc="Camera distortion model to use for astrometric calibration.",
117 name="astrometry_camera",
118 storageClass="Camera",
119 dimensions=("instrument", "physical_filter"),
120 isCalibration=True,
121 minimum=0, # Can fall back on WCS already attached to exposure.
122 )
123 exposures = connectionTypes.Input(
124 doc="Exposure (or two snaps) to be calibrated, and detected and measured on.",
125 name="postISRCCD",
126 storageClass="Exposure",
127 multiple=True, # to handle 1 exposure or 2 snaps
128 dimensions=["instrument", "exposure", "detector"],
129 )
131 background_flat = connectionTypes.PrerequisiteInput(
132 name="flat",
133 doc="Flat calibration frame used for background correction.",
134 storageClass="ExposureF",
135 dimensions=["instrument", "detector", "physical_filter"],
136 isCalibration=True,
137 )
138 illumination_correction = connectionTypes.PrerequisiteInput(
139 name="illuminationCorrection",
140 doc="Illumination correction frame.",
141 storageClass="Exposure",
142 dimensions=["instrument", "detector", "physical_filter"],
143 isCalibration=True,
144 )
146 # outputs
147 initial_stars_schema = connectionTypes.InitOutput(
148 doc="Schema of the output initial stars catalog.",
149 name="initial_stars_schema",
150 storageClass="SourceCatalog",
151 )
153 # TODO DM-38732: We want some kind of flag on Exposures/Catalogs to make
154 # it obvious which components had failed to be computed/persisted.
155 exposure = connectionTypes.Output(
156 doc="Photometrically calibrated, background-subtracted exposure with fitted calibrations and "
157 "summary statistics. To recover the original exposure, first add the background "
158 "(`initial_pvi_background`), and then uncalibrate (divide by `initial_photoCalib_detector`).",
159 name="initial_pvi",
160 storageClass="ExposureF",
161 dimensions=("instrument", "visit", "detector"),
162 )
163 stars = connectionTypes.Output(
164 doc="Catalog of unresolved sources detected on the calibrated exposure.",
165 name="initial_stars_detector",
166 storageClass="ArrowAstropy",
167 dimensions=["instrument", "visit", "detector"],
168 )
169 stars_footprints = connectionTypes.Output(
170 doc="Catalog of unresolved sources detected on the calibrated exposure; "
171 "includes source footprints.",
172 name="initial_stars_footprints_detector",
173 storageClass="SourceCatalog",
174 dimensions=["instrument", "visit", "detector"],
175 )
176 applied_photo_calib = connectionTypes.Output(
177 doc=(
178 "Photometric calibration that was applied to exposure's pixels. "
179 "This connection is disabled when do_calibrate_pixels=False."
180 ),
181 name="initial_photoCalib_detector",
182 storageClass="PhotoCalib",
183 dimensions=("instrument", "visit", "detector"),
184 )
185 background = connectionTypes.Output(
186 doc="Background models estimated during calibration task; calibrated to be in nJy units.",
187 name="initial_pvi_background",
188 storageClass="Background",
189 dimensions=("instrument", "visit", "detector"),
190 )
191 background_to_photometric_ratio = connectionTypes.Output(
192 doc="Ratio of a background-flattened image to a photometric-flattened image. Only persisted "
193 "if do_illumination_correction is True.",
194 name="background_to_photometric_ratio",
195 storageClass="Image",
196 dimensions=("instrument", "visit", "detector"),
197 )
199 # Optional outputs
200 mask = connectionTypes.Output(
201 doc="The mask plane of the calibrated exposure, written as a separate "
202 "output so that it can be passed to downstream tasks even if the "
203 "exposure is removed to save storage.",
204 name="preliminary_visit_mask",
205 storageClass="Mask",
206 dimensions=("instrument", "visit", "detector"),
207 )
208 psf_stars_footprints = connectionTypes.Output(
209 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination; "
210 "includes source footprints.",
211 name="initial_psf_stars_footprints_detector",
212 storageClass="SourceCatalog",
213 dimensions=["instrument", "visit", "detector"],
214 )
215 psf_stars = connectionTypes.Output(
216 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination.",
217 name="initial_psf_stars_detector",
218 storageClass="ArrowAstropy",
219 dimensions=["instrument", "visit", "detector"],
220 )
221 astrometry_matches = connectionTypes.Output(
222 doc="Source to reference catalog matches from the astrometry solver.",
223 name="initial_astrometry_match_detector",
224 storageClass="Catalog",
225 dimensions=("instrument", "visit", "detector"),
226 )
227 photometry_matches = connectionTypes.Output(
228 doc="Source to reference catalog matches from the photometry solver.",
229 name="initial_photometry_match_detector",
230 storageClass="Catalog",
231 dimensions=("instrument", "visit", "detector"),
232 )
234 def __init__(self, *, config=None):
235 super().__init__(config=config)
236 if "psf_stars" not in config.optional_outputs:
237 del self.psf_stars
238 if "psf_stars_footprints" not in config.optional_outputs:
239 del self.psf_stars_footprints
240 if "astrometry_matches" not in config.optional_outputs:
241 del self.astrometry_matches
242 if "photometry_matches" not in config.optional_outputs:
243 del self.photometry_matches
244 if "mask" not in config.optional_outputs:
245 del self.mask
246 if not config.do_calibrate_pixels:
247 del self.applied_photo_calib
248 if not config.do_illumination_correction:
249 del self.background_flat
250 del self.illumination_correction
251 del self.background_to_photometric_ratio
252 if not config.useButlerCamera:
253 del self.camera_model
256class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CalibrateImageConnections):
257 optional_outputs = pexConfig.ListField(
258 doc="Which optional outputs to save (as their connection name)?",
259 dtype=str,
260 # TODO: note somewhere to disable this for benchmarking, but should
261 # we always have it on for production runs?
262 default=["psf_stars", "psf_stars_footprints", "astrometry_matches", "photometry_matches", "mask"],
263 optional=False
264 )
266 # To generate catalog ids consistently across subtasks.
267 id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()
269 snap_combine = pexConfig.ConfigurableField(
270 target=snapCombine.SnapCombineTask,
271 doc="Task to combine two snaps to make one exposure.",
272 )
274 # subtasks used during psf characterization
275 install_simple_psf = pexConfig.ConfigurableField(
276 target=lsst.meas.algorithms.installGaussianPsf.InstallGaussianPsfTask,
277 doc="Task to install a simple PSF model into the input exposure to use "
278 "when detecting bright sources for PSF estimation.",
279 )
280 psf_repair = pexConfig.ConfigurableField(
281 target=repair.RepairTask,
282 doc="Task to repair cosmic rays on the exposure before PSF determination.",
283 )
284 psf_subtract_background = pexConfig.ConfigurableField(
285 target=lsst.meas.algorithms.SubtractBackgroundTask,
286 doc="Task to perform intial background subtraction, before first detection pass.",
287 )
288 psf_detection = pexConfig.ConfigurableField(
289 target=lsst.meas.algorithms.SourceDetectionTask,
290 doc="Task to detect sources for PSF determination."
291 )
292 do_adaptive_threshold_detection = pexConfig.Field(
293 dtype=bool,
294 default=True,
295 doc="Implement the adaptive detection thresholding approach?",
296 )
297 psf_adaptive_threshold_detection = pexConfig.ConfigurableField(
298 target=AdaptiveThresholdDetectionTask,
299 doc="Task to adaptively detect sources for PSF determination.",
300 )
301 psf_source_measurement = pexConfig.ConfigurableField(
302 target=lsst.meas.base.SingleFrameMeasurementTask,
303 doc="Task to measure sources to be used for psf estimation."
304 )
305 psf_measure_psf = pexConfig.ConfigurableField(
306 target=measurePsf.MeasurePsfTask,
307 doc="Task to measure the psf on bright sources."
308 )
309 psf_normalized_calibration_flux = pexConfig.ConfigurableField(
310 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
311 doc="Task to normalize the calibration flux (e.g. compensated tophats) "
312 "for the bright stars used for psf estimation.",
313 )
315 # TODO DM-39203: we can remove aperture correction from this task once we
316 # are using the shape-based star/galaxy code.
317 measure_aperture_correction = pexConfig.ConfigurableField(
318 target=lsst.meas.algorithms.measureApCorr.MeasureApCorrTask,
319 doc="Task to compute the aperture correction from the bright stars."
320 )
322 # subtasks used during star measurement
323 star_background = pexConfig.ConfigurableField(
324 target=lsst.meas.algorithms.SubtractBackgroundTask,
325 doc="Task to perform final background subtraction, just before photoCal.",
326 )
327 star_detection = pexConfig.ConfigurableField(
328 target=lsst.meas.algorithms.SourceDetectionTask,
329 doc="Task to detect stars to return in the output catalog."
330 )
331 star_sky_sources = pexConfig.ConfigurableField(
332 target=lsst.meas.algorithms.SkyObjectsTask,
333 doc="Task to generate sky sources ('empty' regions where there are no detections).",
334 )
335 do_downsample_footprints = pexConfig.Field(
336 dtype=bool,
337 default=False,
338 doc="Downsample footprints prior to deblending to optimize speed?",
339 )
340 downsample_max_footprints = pexConfig.Field(
341 dtype=int,
342 default=1000,
343 doc="Maximum number of non-sky-source footprints to use if do_downsample_footprints is True,",
344 )
345 star_deblend = pexConfig.ConfigurableField(
346 target=lsst.meas.deblender.SourceDeblendTask,
347 doc="Split blended sources into their components."
348 )
349 star_measurement = pexConfig.ConfigurableField(
350 target=lsst.meas.base.SingleFrameMeasurementTask,
351 doc="Task to measure stars to return in the output catalog."
352 )
353 star_normalized_calibration_flux = pexConfig.ConfigurableField(
354 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
355 doc="Task to apply the normalization for calibration fluxes (e.g. compensated tophats) "
356 "for the final output star catalog.",
357 )
358 star_apply_aperture_correction = pexConfig.ConfigurableField(
359 target=lsst.meas.base.ApplyApCorrTask,
360 doc="Task to apply aperture corrections to the selected stars."
361 )
362 star_catalog_calculation = pexConfig.ConfigurableField(
363 target=lsst.meas.base.CatalogCalculationTask,
364 doc="Task to compute extendedness values on the star catalog, "
365 "for the star selector to remove extended sources."
366 )
367 star_set_primary_flags = pexConfig.ConfigurableField(
368 target=lsst.meas.algorithms.setPrimaryFlags.SetPrimaryFlagsTask,
369 doc="Task to add isPrimary to the catalog."
370 )
371 star_selector = lsst.meas.algorithms.sourceSelectorRegistry.makeField(
372 default="science",
373 doc="Task to select reliable stars to use for calibration."
374 )
376 # final calibrations and statistics
377 astrometry = pexConfig.ConfigurableField(
378 target=lsst.meas.astrom.AstrometryTask,
379 doc="Task to perform astrometric calibration to fit a WCS.",
380 )
381 astrometry_ref_loader = pexConfig.ConfigField(
382 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
383 doc="Configuration of reference object loader for astrometric fit.",
384 )
385 photometry = pexConfig.ConfigurableField(
386 target=photoCal.PhotoCalTask,
387 doc="Task to perform photometric calibration to fit a PhotoCalib.",
388 )
389 photometry_ref_loader = pexConfig.ConfigField(
390 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
391 doc="Configuration of reference object loader for photometric fit.",
392 )
393 doMaskDiffractionSpikes = pexConfig.Field(
394 dtype=bool,
395 default=False,
396 doc="If True, load the photometric reference catalog again but select "
397 "only bright stars. Use the bright star catalog to set the SPIKE "
398 "mask for regions likely contaminated by diffraction spikes.",
399 )
400 diffractionSpikeMask = pexConfig.ConfigurableField(
401 target=diffractionSpikeMask.DiffractionSpikeMaskTask,
402 doc="Task to identify and mask the diffraction spikes of bright stars.",
403 )
405 compute_summary_stats = pexConfig.ConfigurableField(
406 target=computeExposureSummaryStats.ComputeExposureSummaryStatsTask,
407 doc="Task to to compute summary statistics on the calibrated exposure."
408 )
410 do_illumination_correction = pexConfig.Field(
411 dtype=bool,
412 default=False,
413 doc="If True, apply the illumination correction. This assumes that the "
414 "input image has already been flat-fielded such that it is suitable "
415 "for background subtraction.",
416 )
418 do_calibrate_pixels = pexConfig.Field(
419 dtype=bool,
420 default=True,
421 doc=(
422 "If True, apply the photometric calibration to the image pixels "
423 "and background model, and attach an identity PhotoCalib to "
424 "the output image to reflect this. If False`, leave the image "
425 "and background uncalibrated and attach the PhotoCalib that maps "
426 "them to physical units."
427 )
428 )
430 do_include_astrometric_errors = pexConfig.Field(
431 dtype=bool,
432 default=True,
433 doc="If True, include astrometric errors in the output catalog.",
434 )
436 run_sattle = pexConfig.Field(
437 dtype=bool,
438 default=False,
439 doc="If True, the sattle service will populate a cache for later use "
440 "in ip_diffim.detectAndMeasure alert verification."
441 )
443 sattle_historical = pexConfig.Field(
444 dtype=bool,
445 default=False,
446 doc="If re-running a pipeline that requires sattle, this should be set "
447 "to True. This will populate sattle's cache with the historic data "
448 "closest in time to the exposure.",
449 )
451 useButlerCamera = pexConfig.Field(
452 dtype=bool,
453 default=False,
454 doc="If True, use a camera distortion model generated elsewhere in "
455 "the pipeline combined with the telescope boresight as a starting "
456 "point for fitting the WCS, instead of using the WCS attached to "
457 "the exposure, which is generated from the boresight and the "
458 "camera model from the obs_* package."
459 )
461 def setDefaults(self):
462 super().setDefaults()
464 # Use a very broad PSF here, to throughly reject CRs.
465 # TODO investigation: a large initial psf guess may make stars look
466 # like CRs for very good seeing images.
467 self.install_simple_psf.fwhm = 4
469 # S/N>=50 sources for PSF determination, but detection to S/N=10.
470 # The thresholdValue sets the minimum flux in a pixel to be included
471 # in the footprint, while peaks are only detected when they are above
472 # thresholdValue * includeThresholdMultiplier. The low thresholdValue
473 # ensures that the footprints are large enough for the noise replacer
474 # to mask out faint undetected neighbors that are not to be measured.
475 self.psf_detection.thresholdValue = 10.0
476 self.psf_detection.includeThresholdMultiplier = 5.0
477 # TODO investigation: Probably want False here, but that may require
478 # tweaking the background spatial scale, to make it small enough to
479 # prevent extra peaks in the wings of bright objects.
480 self.psf_detection.doTempLocalBackground = False
481 self.psf_detection.reEstimateBackground = False
483 # Minimal measurement plugins for PSF determination.
484 # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use
485 # shapeHSM/moments for star/galaxy separation.
486 # TODO DM-39203: we can remove aperture correction from this task
487 # once we are using the shape-based star/galaxy code.
488 self.psf_source_measurement.plugins = ["base_PixelFlags",
489 "base_SdssCentroid",
490 "ext_shapeHSM_HsmSourceMoments",
491 "base_CircularApertureFlux",
492 "base_GaussianFlux",
493 "base_PsfFlux",
494 "base_CompensatedTophatFlux",
495 ]
496 self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
497 # Only measure apertures we need for PSF measurement.
498 self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
499 self.psf_source_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
500 # TODO DM-40843: Remove this line once this is the psfex default.
501 self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \
502 "base_CircularApertureFlux_12_0_instFlux"
504 # No extendeness information available: we need the aperture
505 # corrections to determine that.
506 self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False
507 self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]
508 self.measure_aperture_correction.sourceSelector["science"].flags.bad = []
510 # Detection for good S/N for astrometry/photometry and other
511 # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks.
512 self.star_detection.reEstimateBackground = False
513 self.star_detection.thresholdValue = 5.0
514 self.star_detection.includeThresholdMultiplier = 2.0
515 self.star_measurement.plugins = ["base_PixelFlags",
516 "base_SdssCentroid",
517 "ext_shapeHSM_HsmSourceMoments",
518 "ext_shapeHSM_HsmPsfMoments",
519 "base_GaussianFlux",
520 "base_PsfFlux",
521 "base_CircularApertureFlux",
522 "base_ClassificationSizeExtendedness",
523 "base_CompensatedTophatFlux",
524 ]
525 self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
526 self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
527 # Only measure the apertures we need for star selection.
528 self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
529 self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
531 # We measure and apply the normalization aperture correction with
532 # the psf_normalized_calibration_flux task, and we only apply the
533 # normalization aperture correction for the full list of stars.
534 self.star_normalized_calibration_flux.do_measure_ap_corr = False
536 # Select stars with reliable measurements and no bad flags.
537 self.star_selector["science"].doFlags = True
538 self.star_selector["science"].doUnresolved = True
539 self.star_selector["science"].doSignalToNoise = True
540 self.star_selector["science"].signalToNoise.minimum = 10.0
541 # Keep sky sources in the output catalog, even though they aren't
542 # wanted for calibration.
543 self.star_selector["science"].doSkySources = True
544 # Set the flux and error fields
545 self.star_selector["science"].signalToNoise.fluxField = "slot_CalibFlux_instFlux"
546 self.star_selector["science"].signalToNoise.errField = "slot_CalibFlux_instFluxErr"
548 # Use the affine WCS fitter (assumes we have a good camera geometry).
549 self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)
550 # phot_g_mean is the primary Gaia band for all input bands.
551 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"
553 # Only reject sky sources; we already selected good stars.
554 self.astrometry.sourceSelector["science"].doFlags = True
555 self.astrometry.sourceSelector["science"].flags.good = [] # ["calib_psf_candidate"]
556 # self.astrometry.sourceSelector["science"].flags.bad = []
557 self.astrometry.sourceSelector["science"].doUnresolved = False
558 self.astrometry.sourceSelector["science"].doIsolated = False
559 self.astrometry.sourceSelector["science"].doRequirePrimary = False
560 self.photometry.match.sourceSelection.doFlags = True
561 self.photometry.match.sourceSelection.flags.bad = ["sky_source"]
562 # Unset the (otherwise reasonable, but we've already made the
563 # selections we want above) selection settings in PhotoCalTask.
564 self.photometry.match.sourceSelection.doRequirePrimary = False
565 self.photometry.match.sourceSelection.doUnresolved = False
567 def validate(self):
568 super().validate()
570 # Ensure that the normalization calibration flux tasks
571 # are configured correctly.
572 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:
573 msg = ("psf_normalized_calibration_flux task must be configured with do_measure_ap_corr=True "
574 "or else the normalization and calibration flux will not be properly measured.")
575 raise pexConfig.FieldValidationError(
576 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,
577 )
578 if self.star_normalized_calibration_flux.do_measure_ap_corr:
579 msg = ("star_normalized_calibration_flux task must be configured with do_measure_ap_corr=False "
580 "to apply the previously measured normalization to the full catalog of calibration "
581 "fluxes.")
582 raise pexConfig.FieldValidationError(
583 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,
584 )
586 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,
587 # because they'd be running too early to pick up the fitted PhotoCalib
588 # and WCS.
589 if "base_LocalWcs" in self.psf_source_measurement.plugins.names:
590 raise pexConfig.FieldValidationError(
591 CalibrateImageConfig.psf_source_measurement,
592 self,
593 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
594 )
595 if "base_LocalWcs" in self.star_measurement.plugins.names:
596 raise pexConfig.FieldValidationError(
597 CalibrateImageConfig.star_measurement,
598 self,
599 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
600 )
601 if "base_LocalPhotoCalib" in self.psf_source_measurement.plugins.names:
602 raise pexConfig.FieldValidationError(
603 CalibrateImageConfig.psf_source_measurement,
604 self,
605 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
606 "as it would be run before the photometry fit."
607 )
608 if "base_LocalPhotoCalib" in self.star_measurement.plugins.names:
609 raise pexConfig.FieldValidationError(
610 CalibrateImageConfig.star_measurement,
611 self,
612 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
613 "as it would be run before the photometry fit."
614 )
616 # Check for illumination correction and background consistency.
617 if self.do_illumination_correction:
618 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:
619 raise pexConfig.FieldValidationError(
620 CalibrateImageConfig.psf_subtract_background,
621 self,
622 "CalibrateImageTask.psf_subtract_background must be configured with "
623 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
624 )
625 if self.psf_detection.reEstimateBackground:
626 if not self.psf_detection.doApplyFlatBackgroundRatio:
627 raise pexConfig.FieldValidationError(
628 CalibrateImageConfig.psf_detection,
629 self,
630 "CalibrateImageTask.psf_detection background must be configured with "
631 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
632 )
633 if self.star_detection.reEstimateBackground:
634 if not self.star_detection.doApplyFlatBackgroundRatio:
635 raise pexConfig.FieldValidationError(
636 CalibrateImageConfig.star_detection,
637 self,
638 "CalibrateImageTask.star_detection background must be configured with "
639 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
640 )
641 if not self.star_background.doApplyFlatBackgroundRatio:
642 raise pexConfig.FieldValidationError(
643 CalibrateImageConfig.star_background,
644 self,
645 "CalibrateImageTask.star_background background must be configured with "
646 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
647 )
649 if self.run_sattle:
650 if not os.getenv("SATTLE_URI_BASE"):
651 raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self,
652 "Sattle requested but URI environment variable not set.")
654 if not self.do_adaptive_threshold_detection:
655 if not self.psf_detection.reEstimateBackground:
656 raise pexConfig.FieldValidationError(
657 CalibrateImageConfig.psf_detection,
658 self,
659 "If not using the adaptive threshold detection implementation (i.e. "
660 "config.do_adaptive_threshold_detection = False), CalibrateImageTask.psf_detection "
661 "background must be configured with "
662 "reEstimateBackground = True to maintain original behavior."
663 )
664 if not self.star_detection.reEstimateBackground:
665 raise pexConfig.FieldValidationError(
666 CalibrateImageConfig.psf_detection,
667 self,
668 "If not using the adaptive threshold detection implementation "
669 "(i.e. config.do_adaptive_threshold_detection = False), "
670 "CalibrateImageTask.star_detection background must be configured "
671 "with reEstimateBackground = True to maintain original behavior."
672 )
675class CalibrateImageTask(pipeBase.PipelineTask):
676 """Compute the PSF, aperture corrections, astrometric and photometric
677 calibrations, and summary statistics for a single science exposure, and
678 produce a catalog of brighter stars that were used to calibrate it.
680 Parameters
681 ----------
682 initial_stars_schema : `lsst.afw.table.Schema`
683 Schema of the initial_stars output catalog.
684 """
685 _DefaultName = "calibrateImage"
686 ConfigClass = CalibrateImageConfig
688 def __init__(self, initial_stars_schema=None, **kwargs):
689 super().__init__(**kwargs)
690 self.makeSubtask("snap_combine")
692 # PSF determination subtasks
693 self.makeSubtask("install_simple_psf")
694 self.makeSubtask("psf_repair")
695 self.makeSubtask("psf_subtract_background")
696 self.psf_schema = afwTable.SourceTable.makeMinimalSchema()
697 self.psf_schema.addField(
698 'psf_max_value',
699 type=np.float32,
700 doc="PSF max value.",
701 )
702 afwTable.CoordKey.addErrorFields(self.psf_schema)
703 self.makeSubtask("psf_detection", schema=self.psf_schema)
704 self.makeSubtask("psf_source_measurement", schema=self.psf_schema)
705 self.makeSubtask("psf_adaptive_threshold_detection")
706 self.makeSubtask("psf_measure_psf", schema=self.psf_schema)
707 self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema)
709 self.makeSubtask("measure_aperture_correction", schema=self.psf_schema)
710 self.makeSubtask("astrometry", schema=self.psf_schema)
712 # star measurement subtasks
713 if initial_stars_schema is None:
714 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()
716 # These fields let us track which sources were used for psf modeling,
717 # astrometric fitting, and aperture correction calculations.
718 self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved",
719 "calib_astrometry_used",
720 # TODO DM-39203:
721 # these can be removed once apcorr is gone.
722 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",
723 "apcorr_base_PsfFlux_used",)
724 for field in self.psf_fields:
725 item = self.psf_schema.find(field)
726 initial_stars_schema.addField(item.getField())
727 id_type = self.psf_schema["id"].asField().getTypeString()
728 psf_max_value_type = self.psf_schema['psf_max_value'].asField().getTypeString()
729 initial_stars_schema.addField("psf_id",
730 type=id_type,
731 doc="id of this source in psf_stars; 0 if there is no match.")
732 initial_stars_schema.addField("psf_max_value",
733 type=psf_max_value_type,
734 doc="Maximum value in the star image used to train PSF.")
736 afwTable.CoordKey.addErrorFields(initial_stars_schema)
737 self.makeSubtask("star_background")
738 self.makeSubtask("star_detection", schema=initial_stars_schema)
739 self.makeSubtask("star_sky_sources", schema=initial_stars_schema)
740 self.makeSubtask("star_deblend", schema=initial_stars_schema)
741 self.makeSubtask("star_measurement", schema=initial_stars_schema)
742 self.makeSubtask("star_normalized_calibration_flux", schema=initial_stars_schema)
744 self.makeSubtask("star_apply_aperture_correction", schema=initial_stars_schema)
745 self.makeSubtask("star_catalog_calculation", schema=initial_stars_schema)
746 self.makeSubtask("star_set_primary_flags", schema=initial_stars_schema, isSingleFrame=True)
747 self.makeSubtask("star_selector")
748 self.makeSubtask("photometry", schema=initial_stars_schema)
749 if self.config.doMaskDiffractionSpikes:
750 self.makeSubtask("diffractionSpikeMask")
751 self.makeSubtask("compute_summary_stats")
753 # The final catalog will have calibrated flux columns, which we add to
754 # the init-output schema by calibrating our zero-length catalog with an
755 # arbitrary dummy PhotoCalib. We also use this schema to initialze
756 # the stars catalog in order to ensure it's the same even when we hit
757 # an error (and write partial outputs) before calibrating the catalog
758 # - note that calibrateCatalog will happily reuse existing output
759 # columns.
760 dummy_photo_calib = afwImage.PhotoCalib(1.0, 0, bbox=lsst.geom.Box2I())
761 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(
762 afwTable.SourceCatalog(initial_stars_schema)
763 )
765 def runQuantum(self, butlerQC, inputRefs, outputRefs):
766 inputs = butlerQC.get(inputRefs)
767 exposures = inputs.pop("exposures")
769 id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)
771 astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
772 dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],
773 refCats=inputs.pop("astrometry_ref_cat"),
774 name=self.config.connections.astrometry_ref_cat,
775 config=self.config.astrometry_ref_loader, log=self.log)
776 self.astrometry.setRefObjLoader(astrometry_loader)
778 photometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
779 dataIds=[ref.datasetRef.dataId for ref in inputRefs.photometry_ref_cat],
780 refCats=inputs.pop("photometry_ref_cat"),
781 name=self.config.connections.photometry_ref_cat,
782 config=self.config.photometry_ref_loader, log=self.log)
783 self.photometry.match.setRefObjLoader(photometry_loader)
785 if self.config.doMaskDiffractionSpikes:
786 # Use the same photometry reference catalog for the bright star
787 # mask.
788 self.diffractionSpikeMask.setRefObjLoader(photometry_loader)
790 if self.config.do_illumination_correction:
791 background_flat = inputs.pop("background_flat")
792 illumination_correction = inputs.pop("illumination_correction")
793 else:
794 background_flat = None
795 illumination_correction = None
797 camera_model = None
798 if self.config.useButlerCamera:
799 if "camera_model" in inputs:
800 camera_model = inputs.pop("camera_model")
801 else:
802 self.log.warning("useButlerCamera=True, but camera is not available for filter %s. The "
803 "astrometry fit will use the WCS already attached to the exposure.",
804 exposures[0].filter.bandLabel)
806 # This should not happen with a properly configured execution context.
807 assert not inputs, "runQuantum got more inputs than expected"
809 # Specify the fields that `annotate` needs below, to ensure they
810 # exist, even as None.
811 result = pipeBase.Struct(
812 exposure=None,
813 stars_footprints=None,
814 psf_stars_footprints=None,
815 background_to_photometric_ratio=None,
816 )
817 try:
818 self.run(
819 exposures=exposures,
820 result=result,
821 id_generator=id_generator,
822 background_flat=background_flat,
823 illumination_correction=illumination_correction,
824 camera_model=camera_model,
825 )
826 except pipeBase.AlgorithmError as e:
827 error = pipeBase.AnnotatedPartialOutputsError.annotate(
828 e,
829 self,
830 result.exposure,
831 result.psf_stars_footprints,
832 result.stars_footprints,
833 log=self.log
834 )
835 butlerQC.put(result, outputRefs)
836 raise error from e
838 butlerQC.put(result, outputRefs)
840 @timeMethod
841 def run(
842 self,
843 *,
844 exposures,
845 id_generator=None,
846 result=None,
847 background_flat=None,
848 illumination_correction=None,
849 camera_model=None,
850 ):
851 """Find stars and perform psf measurement, then do a deeper detection
852 and measurement and calibrate astrometry and photometry from that.
854 Parameters
855 ----------
856 exposures : `lsst.afw.image.Exposure` or \
857 `list` [`lsst.afw.image.Exposure`]
858 Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter.
859 Modified in-place during processing if only one is passed.
860 If two exposures are passed, treat them as snaps and combine
861 before doing further processing.
862 id_generator : `lsst.meas.base.IdGenerator`, optional
863 Object that generates source IDs and provides random seeds.
864 result : `lsst.pipe.base.Struct`, optional
865 Result struct that is modified to allow saving of partial outputs
866 for some failure conditions. If the task completes successfully,
867 this is also returned.
868 background_flat : `lsst.afw.image.Exposure`, optional
869 Background flat-field image.
870 illumination_correction : `lsst.afw.image.Exposure`, optional
871 Illumination correction image.
872 camera_model : `lsst.afw.cameraGeom.Camera`, optional
873 Camera to be used if constructing updated WCS.
875 Returns
876 -------
877 result : `lsst.pipe.base.Struct`
878 Results as a struct with attributes:
880 ``exposure``
881 Calibrated exposure, with pixels in nJy units.
882 (`lsst.afw.image.Exposure`)
883 ``stars``
884 Stars that were used to calibrate the exposure, with
885 calibrated fluxes and magnitudes.
886 (`astropy.table.Table`)
887 ``stars_footprints``
888 Footprints of stars that were used to calibrate the exposure.
889 (`lsst.afw.table.SourceCatalog`)
890 ``psf_stars``
891 Stars that were used to determine the image PSF.
892 (`astropy.table.Table`)
893 ``psf_stars_footprints``
894 Footprints of stars that were used to determine the image PSF.
895 (`lsst.afw.table.SourceCatalog`)
896 ``background``
897 Background that was fit to the exposure when detecting
898 ``stars``. (`lsst.afw.math.BackgroundList`)
899 ``applied_photo_calib``
900 Photometric calibration that was fit to the star catalog and
901 applied to the exposure. (`lsst.afw.image.PhotoCalib`)
902 This is `None` if ``config.do_calibrate_pixels`` is `False`.
903 ``astrometry_matches``
904 Reference catalog stars matches used in the astrometric fit.
905 (`list` [`lsst.afw.table.ReferenceMatch`] or
906 `lsst.afw.table.BaseCatalog`).
907 ``photometry_matches``
908 Reference catalog stars matches used in the photometric fit.
909 (`list` [`lsst.afw.table.ReferenceMatch`] or
910 `lsst.afw.table.BaseCatalog`).
911 ``mask``
912 Copy of the mask plane of `exposure`.
913 (`lsst.afw.image.Mask`)
914 """
915 if result is None:
916 result = pipeBase.Struct()
917 if id_generator is None:
918 id_generator = lsst.meas.base.IdGenerator()
920 result.exposure = self.snap_combine.run(exposures).exposure
921 self._recordMaskedPixelFractions(result.exposure)
922 self.log.info("Initial PhotoCalib: %s", result.exposure.getPhotoCalib())
924 result.exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = False
926 # Check input image processing.
927 if self.config.do_illumination_correction:
928 if not result.exposure.metadata.get("LSST ISR FLAT APPLIED", False):
929 raise pipeBase.InvalidQuantumError(
930 "Cannot use do_illumination_correction with an image that has not had a flat applied",
931 )
933 if camera_model:
934 if camera_model.get(result.exposure.detector.getId()):
935 self.log.info("Updating WCS with the provided camera model.")
936 self._update_wcs_with_camera_model(result.exposure, camera_model)
937 else:
938 self.log.warning(
939 "useButlerCamera=True, but detector %s is not available in the provided camera. The "
940 "astrometry fit will use the WCS already attached to the exposure.",
941 result.exposure.detector.getId())
943 result.background = None
944 summary_stat_catalog = None
945 # Some exposure components are set to initial placeholder objects
946 # while we try to bootstrap them. If we fail before we fit for them,
947 # we want to reset those components to None so the placeholders don't
948 # masquerade as the real thing.
949 have_fit_psf = False
950 have_fit_astrometry = False
951 have_fit_photometry = False
952 try:
953 result.background_to_photometric_ratio = self._apply_illumination_correction(
954 result.exposure,
955 background_flat,
956 illumination_correction,
957 )
959 result.psf_stars_footprints, result.background, _, adaptive_det_res_struct = self._compute_psf(
960 result.exposure,
961 id_generator,
962 background_to_photometric_ratio=result.background_to_photometric_ratio,
963 )
964 have_fit_psf = True
966 # Check if all centroids have been flagged. This should happen
967 # after the PSF is fit and recorded because even a junky PSF
968 # model is informative.
969 if result.psf_stars_footprints["slot_Centroid_flag"].all():
970 psf_shape = result.exposure.psf.computeShape(result.exposure.psf.getAveragePosition())
971 raise AllCentroidsFlaggedError(
972 n_sources=len(result.psf_stars_footprints),
973 psf_shape_ixx=psf_shape.getIxx(),
974 psf_shape_iyy=psf_shape.getIyy(),
975 psf_shape_ixy=psf_shape.getIxy(),
976 psf_size=psf_shape.getDeterminantRadius(),
977 )
979 if result.background is None:
980 result.background = afwMath.BackgroundList()
982 self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)
983 result.psf_stars = result.psf_stars_footprints.asAstropy()
984 # Run astrometry using PSF candidate stars.
985 # Update "the psf_stars" source cooordinates with the current wcs.
986 afwTable.updateSourceCoords(
987 result.exposure.wcs,
988 sourceList=result.psf_stars_footprints,
989 include_covariance=self.config.do_include_astrometric_errors,
990 )
991 astrometry_matches, astrometry_meta = self._fit_astrometry(
992 result.exposure, result.psf_stars_footprints
993 )
994 num_astrometry_matches = len(astrometry_matches)
995 self.metadata["astrometry_matches_count"] = num_astrometry_matches
996 if "astrometry_matches" in self.config.optional_outputs:
997 result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches,
998 astrometry_meta)
999 result.psf_stars = result.psf_stars_footprints.asAstropy()
1001 if self.config.do_adaptive_threshold_detection:
1002 self._remeasure_star_background(
1003 result,
1004 background_to_photometric_ratio=result.background_to_photometric_ratio,
1005 )
1007 # Run the stars_detection subtask for the photometric calibration.
1008 result.stars_footprints = self._find_stars(
1009 result.exposure,
1010 result.background,
1011 id_generator,
1012 background_to_photometric_ratio=result.background_to_photometric_ratio,
1013 adaptive_det_res_struct=adaptive_det_res_struct,
1014 num_astrometry_matches=num_astrometry_matches,
1015 )
1016 psf = result.exposure.getPsf()
1017 psfSigma = psf.computeShape(result.exposure.getBBox().getCenter()).getDeterminantRadius()
1018 self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints,
1019 psfSigma=psfSigma)
1021 # Update the "stars" source cooordinates with the current wcs.
1022 afwTable.updateSourceCoords(
1023 result.exposure.wcs,
1024 sourceList=result.stars_footprints,
1025 include_covariance=self.config.do_include_astrometric_errors,
1026 )
1028 summary_stat_catalog = result.stars_footprints
1029 result.stars = result.stars_footprints.asAstropy()
1030 self.metadata["star_count"] = np.sum(~result.stars["sky_source"])
1032 # Validate the astrometric fit. Send in the stars_footprints
1033 # catalog so that its coords get set to NaN if the fit is deemed
1034 # a failure.
1035 self.astrometry.check(result.exposure, result.stars_footprints, len(astrometry_matches))
1036 result.stars = result.stars_footprints.asAstropy()
1037 have_fit_astrometry = True
1039 result.stars_footprints, photometry_matches, \
1040 photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints)
1041 have_fit_photometry = True
1042 self.metadata["photometry_matches_count"] = len(photometry_matches)
1043 # fit_photometry returns a new catalog, so we need a new astropy
1044 # table view.
1045 result.stars = result.stars_footprints.asAstropy()
1046 # summary stats don't make use of the calibrated fluxes, but we
1047 # might as well use the best catalog we've got in case that
1048 # changes, and help the old one get garbage-collected.
1049 summary_stat_catalog = result.stars_footprints
1050 if "photometry_matches" in self.config.optional_outputs:
1051 result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,
1052 photometry_meta)
1053 if self.config.doMaskDiffractionSpikes:
1054 self.diffractionSpikeMask.run(result.exposure)
1055 if "mask" in self.config.optional_outputs:
1056 result.mask = result.exposure.mask.clone()
1057 except pipeBase.AlgorithmError:
1058 if not have_fit_psf:
1059 result.exposure.setPsf(None)
1060 if not have_fit_astrometry:
1061 result.exposure.setWcs(None)
1062 if not have_fit_photometry:
1063 result.exposure.setPhotoCalib(None)
1064 # Summary stat calculations can handle missing components
1065 # gracefully, but we want to run them as late as possible (but
1066 # still before we calibrate pixels, if we do that at all).
1067 # So we run them after we succeed or if we get an AlgorithmError.
1068 # We intentionally don't use 'finally' here because we don't
1069 # want to run them if we get some other kind of error.
1070 self._summarize(result.exposure, summary_stat_catalog, result.background)
1071 raise
1072 else:
1073 self._summarize(result.exposure, summary_stat_catalog, result.background)
1075 if self.config.do_calibrate_pixels:
1076 self._apply_photometry(
1077 result.exposure,
1078 result.background,
1079 background_to_photometric_ratio=result.background_to_photometric_ratio,
1080 )
1081 result.applied_photo_calib = photo_calib
1082 else:
1083 result.applied_photo_calib = None
1085 if self.config.run_sattle:
1086 # send boresight and timing information to sattle so the cache
1087 # is populated by the time we reach ip_diffim detectAndMeasure.
1088 try:
1089 populate_sattle_visit_cache(result.exposure.getInfo().getVisitInfo(),
1090 historical=self.config.sattle_historical)
1091 self.log.info('Successfully triggered load of sattle visit cache')
1092 except requests.exceptions.HTTPError:
1093 self.log.exception("Sattle visit cache update failed; continuing with image processing")
1095 return result
1097 def _apply_illumination_correction(self, exposure, background_flat, illumination_correction):
1098 """Apply the illumination correction to a background-flattened image.
1100 Parameters
1101 ----------
1102 exposure : `lsst.afw.image.Exposure`
1103 Exposure to convert to a photometric-flattened image.
1104 background_flat : `lsst.afw.image.Exposure`
1105 Flat image that had previously been applied to exposure.
1106 illumination_correction : `lsst.afw.image.Exposure`
1107 Illumination correction image to convert to photometric-flattened
1108 image.
1110 Returns
1111 -------
1112 background_to_photometric_ratio : `lsst.afw.image.Image`
1113 Ratio image to convert a photometric-flattened image to/from
1114 a background-flattened image. Will be None if task not
1115 configured to use the illumination correction.
1116 """
1117 if not self.config.do_illumination_correction:
1118 return None
1120 # From a raw image to a background-flattened image, we have:
1121 # bfi = image / background_flat
1122 # From a raw image to a photometric-flattened image, we have:
1123 # pfi = image / reference_flux_flat
1124 # pfi = image / (dome_flat * illumination_correction),
1125 # where the illumination correction contains the jacobian
1126 # of the wcs, converting to fluence units.
1127 # Currently background_flat == dome_flat, so we have for the
1128 # "background_to_photometric_ratio", the ratio of the background-
1129 # flattened image to the photometric-flattened image:
1130 # bfi / pfi = illumination_correction.
1132 background_to_photometric_ratio = illumination_correction.image.clone()
1134 # Dividing the ratio will convert a background-flattened image to
1135 # a photometric-flattened image.
1136 exposure.maskedImage /= background_to_photometric_ratio
1138 exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = True
1140 return background_to_photometric_ratio
1142 def _compute_psf(self, exposure, id_generator, background_to_photometric_ratio=None):
1143 """Find bright sources detected on an exposure and fit a PSF model to
1144 them, repairing likely cosmic rays before detection.
1146 Repair, detect, measure, and compute PSF twice, to ensure the PSF
1147 model does not include contributions from cosmic rays.
1149 Parameters
1150 ----------
1151 exposure : `lsst.afw.image.Exposure`
1152 Exposure to detect and measure bright stars on.
1153 id_generator : `lsst.meas.base.IdGenerator`
1154 Object that generates source IDs and provides random seeds.
1155 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1156 Image to convert photometric-flattened image to
1157 background-flattened image.
1159 Returns
1160 -------
1161 sources : `lsst.afw.table.SourceCatalog`
1162 Catalog of detected bright sources.
1163 background : `lsst.afw.math.BackgroundList`
1164 Background that was fit to the exposure during detection.
1165 cell_set : `lsst.afw.math.SpatialCellSet`
1166 PSF candidates returned by the psf determiner.
1167 """
1168 def log_psf(msg, addToMetadata=False):
1169 """Log the parameters of the psf and background, with a prepended
1170 message. There is also the option to add the PSF sigma to the task
1171 metadata.
1173 Parameters
1174 ----------
1175 msg : `str`
1176 Message to prepend the log info with.
1177 addToMetadata : `bool`, optional
1178 Whether to add the final psf sigma value to the task
1179 metadata (the default is False).
1180 """
1181 position = exposure.psf.getAveragePosition()
1182 sigma = exposure.psf.computeShape(position).getDeterminantRadius()
1183 dimensions = exposure.psf.computeImage(position).getDimensions()
1184 if background is not None:
1185 median_background = np.median(background.getImage().array)
1186 else:
1187 median_background = 0.0
1188 self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",
1189 msg, sigma, dimensions, median_background)
1190 if addToMetadata:
1191 self.metadata["final_psf_sigma"] = sigma
1193 self.log.info("First pass detection with Guassian PSF FWHM=%s pixels",
1194 self.config.install_simple_psf.fwhm)
1195 self.install_simple_psf.run(exposure=exposure)
1197 background = self.psf_subtract_background.run(
1198 exposure=exposure,
1199 backgroundToPhotometricRatio=background_to_photometric_ratio,
1200 ).background
1201 log_psf("Initial PSF:")
1202 self.psf_repair.run(exposure=exposure, keepCRs=True)
1204 table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())
1205 if not self.config.do_adaptive_threshold_detection:
1206 adaptive_det_res_struct = None
1207 # Re-estimate the background during this detection step, so that
1208 # measurement uses the most accurate background-subtraction.
1209 detections = self.psf_detection.run(
1210 table=table,
1211 exposure=exposure,
1212 background=background,
1213 backgroundToPhotometricRatio=background_to_photometric_ratio,
1214 )
1215 else:
1216 initialThreshold = self.config.psf_detection.thresholdValue
1217 initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier
1218 adaptive_det_res_struct = self.psf_adaptive_threshold_detection.run(
1219 table, exposure,
1220 initialThreshold=initialThreshold,
1221 initialThresholdMultiplier=initialThresholdMultiplier,
1222 )
1223 detections = adaptive_det_res_struct.detections
1225 self.metadata["initial_psf_positive_footprint_count"] = detections.numPos
1226 self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg
1227 self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks
1228 self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks
1229 self.psf_source_measurement.run(detections.sources, exposure)
1230 psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
1232 # This 2nd round of PSF fitting has been deemed unnecessary (and
1233 # sometimes even causing harm), so it is being skipped for the
1234 # adaptive threshold logic, but is left here for now to maintain
1235 # consistency with the previous logic for any users wanting to
1236 # revert to it.
1237 if not self.config.do_adaptive_threshold_detection:
1238 # Replace the initial PSF with something simpler for the second
1239 # repair/detect/measure/measure_psf step: this can help it
1240 # converge.
1241 self.install_simple_psf.run(exposure=exposure)
1243 log_psf("Rerunning with simple PSF:")
1244 # TODO investigation: Should we only re-run repair here, to use the
1245 # new PSF? Maybe we *do* need to re-run measurement with PsfFlux,
1246 # to use the fitted PSF?
1247 # TODO investigation: do we need a separate measurement task here
1248 # for the post-psf_measure_psf step, since we only want to do
1249 # PsfFlux and GaussianFlux *after* we have a PSF? Maybe that's not
1250 # relevant once DM-39203 is merged?
1251 self.psf_repair.run(exposure=exposure, keepCRs=True)
1252 # Re-estimate the background during this detection step, so that
1253 # measurement uses the most accurate background-subtraction.
1254 detections = self.psf_detection.run(
1255 table=table,
1256 exposure=exposure,
1257 background=background,
1258 backgroundToPhotometricRatio=background_to_photometric_ratio,
1259 )
1260 self.psf_source_measurement.run(detections.sources, exposure)
1261 psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
1263 self.metadata["simple_psf_positive_footprint_count"] = detections.numPos
1264 self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg
1265 self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks
1266 self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks
1267 log_psf("Final PSF:", addToMetadata=True)
1269 # Final repair with final PSF, removing cosmic rays this time.
1270 self.psf_repair.run(exposure=exposure)
1271 # Final measurement with the CRs removed.
1272 self.psf_source_measurement.run(detections.sources, exposure)
1274 # PSF is set on exposure; candidates are returned to use for
1275 # calibration flux normalization and aperture corrections.
1276 return detections.sources, background, psf_result.cellSet, adaptive_det_res_struct
1278 def _measure_aperture_correction(self, exposure, bright_sources):
1279 """Measure and set the ApCorrMap on the Exposure, using
1280 previously-measured bright sources.
1282 This function first normalizes the calibration flux and then
1283 the full set of aperture corrections are measured relative
1284 to this normalized calibration flux.
1286 Parameters
1287 ----------
1288 exposure : `lsst.afw.image.Exposure`
1289 Exposure to set the ApCorrMap on.
1290 bright_sources : `lsst.afw.table.SourceCatalog`
1291 Catalog of detected bright sources; modified to include columns
1292 necessary for point source determination for the aperture
1293 correction calculation.
1294 """
1295 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(
1296 exposure=exposure,
1297 catalog=bright_sources,
1298 ).ap_corr_map
1300 ap_corr_map = self.measure_aperture_correction.run(exposure, bright_sources).apCorrMap
1302 # Need to merge the aperture correction map from the normalization.
1303 for key in norm_ap_corr_map:
1304 ap_corr_map[key] = norm_ap_corr_map[key]
1306 exposure.info.setApCorrMap(ap_corr_map)
1308 def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None,
1309 adaptive_det_res_struct=None, num_astrometry_matches=None):
1310 """Detect stars on an exposure that has a PSF model, and measure their
1311 PSF, circular aperture, compensated gaussian fluxes.
1313 Parameters
1314 ----------
1315 exposure : `lsst.afw.image.Exposure`
1316 Exposure to detect and measure stars on.
1317 background : `lsst.afw.math.BackgroundList`
1318 Background that was fit to the exposure during detection;
1319 modified in-place during subsequent detection.
1320 id_generator : `lsst.meas.base.IdGenerator`
1321 Object that generates source IDs and provides random seeds.
1322 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1323 Image to convert photometric-flattened image to
1324 background-flattened image.
1326 Returns
1327 -------
1328 stars : `SourceCatalog`
1329 Sources that are very likely to be stars, with a limited set of
1330 measurements performed on them.
1331 """
1332 table = afwTable.SourceTable.make(self.initial_stars_schema.schema,
1333 id_generator.make_table_id_factory())
1335 maxAdaptiveDetIter = 8
1336 threshFactor = 0.2
1337 if adaptive_det_res_struct is not None:
1338 for adaptiveDetIter in range(maxAdaptiveDetIter):
1339 adaptiveDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
1340 adaptiveDetectionConfig.reEstimateBackground = False
1341 adaptiveDetectionConfig.includeThresholdMultiplier = 2.0
1342 psfThreshold = (
1343 adaptive_det_res_struct.thresholdValue*adaptive_det_res_struct.includeThresholdMultiplier
1344 )
1345 adaptiveDetectionConfig.thresholdValue = max(
1346 self.config.star_detection.thresholdValue,
1347 threshFactor*psfThreshold/adaptiveDetectionConfig.includeThresholdMultiplier
1348 )
1349 self.log.info("Using adaptive threshold detection (nIter = %d) with "
1350 "thresholdValue = %.2f and multiplier = %.1f",
1351 adaptiveDetIter, adaptiveDetectionConfig.thresholdValue,
1352 adaptiveDetectionConfig.includeThresholdMultiplier)
1353 adaptiveDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
1354 config=adaptiveDetectionConfig
1355 )
1356 detections = adaptiveDetectionTask.run(
1357 table=table,
1358 exposure=exposure,
1359 background=background,
1360 backgroundToPhotometricRatio=background_to_photometric_ratio,
1361 )
1362 nFootprint = len(detections.sources)
1363 nPeak = 0
1364 nIsolated = 0
1365 for src in detections.sources:
1366 nPeakSrc = len(src.getFootprint().getPeaks())
1367 if nPeakSrc == 1:
1368 nIsolated += 1
1369 nPeak += nPeakSrc
1370 minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches))
1371 if nFootprint > 0:
1372 self.log.info("Adaptive threshold detection _find_stars nIter: %d; "
1373 "nPeak/nFootprint = %.2f (max is 800); nIsolated = %d (min is %.1f).",
1374 adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated)
1375 if nPeak/nFootprint > 800 or nIsolated < minIsolated:
1376 threshFactor = max(0.01, 1.5*threshFactor)
1377 self.log.warning("nPeak/nFootprint = %.2f (max is 800); nIsolated = %d "
1378 "(min is %.1f).", nPeak/nFootprint, nIsolated, minIsolated)
1379 else:
1380 break
1381 else:
1382 threshFactor *= 0.75
1383 self.log.warning("No footprints detected on image. Decreasing threshold "
1384 "factor to %.2f. and rerunning.", threshFactor)
1385 else:
1386 # Re-estimate the background during this detection step, so that
1387 # measurement uses the most accurate background-subtraction.
1388 detections = self.star_detection.run(
1389 table=table,
1390 exposure=exposure,
1391 background=background,
1392 backgroundToPhotometricRatio=background_to_photometric_ratio,
1393 )
1394 sources = detections.sources
1395 self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources)
1397 n_sky_sources = np.sum(sources["sky_source"])
1398 if (self.config.do_downsample_footprints
1399 and (len(sources) - n_sky_sources) > self.config.downsample_max_footprints):
1400 if exposure.info.id is None:
1401 self.log.warning("Exposure does not have a proper id; using 0 seed for downsample.")
1402 seed = 0
1403 else:
1404 seed = exposure.info.id & 0xFFFFFFFF
1406 gen = np.random.RandomState(seed)
1408 # Isolate the sky sources before downsampling.
1409 indices = np.arange(len(sources))[~sources["sky_source"]]
1410 indices = gen.choice(
1411 indices,
1412 size=self.config.downsample_max_footprints,
1413 replace=False,
1414 )
1415 skyIndices, = np.where(sources["sky_source"])
1416 indices = np.concatenate((indices, skyIndices))
1418 self.log.info("Downsampling from %d to %d non-sky-source footprints.", len(sources), len(indices))
1420 sel = np.zeros(len(sources), dtype=bool)
1421 sel[indices] = True
1422 sources = sources[sel]
1424 # TODO investigation: Could this deblender throw away blends of
1425 # non-PSF sources?
1426 self.star_deblend.run(exposure=exposure, sources=sources)
1427 # The deblender may not produce a contiguous catalog; ensure
1428 # contiguity for subsequent tasks.
1429 if not sources.isContiguous():
1430 sources = sources.copy(deep=True)
1432 # Measure everything, and use those results to select only stars.
1433 self.star_measurement.run(sources, exposure)
1434 self.metadata["post_deblend_source_count"] = np.sum(~sources["sky_source"])
1435 self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"])
1436 self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"])
1438 # Run the normalization calibration flux task to apply the
1439 # normalization correction to create normalized
1440 # calibration fluxes.
1441 self.star_normalized_calibration_flux.run(exposure=exposure, catalog=sources)
1442 self.star_apply_aperture_correction.run(sources, exposure.apCorrMap)
1443 self.star_catalog_calculation.run(sources)
1444 self.star_set_primary_flags.run(sources)
1446 result = self.star_selector.run(sources)
1447 # The star selector may not produce a contiguous catalog.
1448 if not result.sourceCat.isContiguous():
1449 return result.sourceCat.copy(deep=True)
1450 else:
1451 return result.sourceCat
1453 def _match_psf_stars(self, psf_stars, stars, psfSigma=None):
1454 """Match calibration stars to psf stars, to identify which were psf
1455 candidates, and which were used or reserved during psf measurement
1456 and the astrometric fit.
1458 Parameters
1459 ----------
1460 psf_stars : `lsst.afw.table.SourceCatalog`
1461 PSF candidate stars that were sent to the psf determiner and
1462 used in the astrometric fit. Used to populate psf and astrometry
1463 related flag fields.
1464 stars : `lsst.afw.table.SourceCatalog`
1465 Stars that will be used for calibration; psf-related fields will
1466 be updated in-place.
1468 Notes
1469 -----
1470 This code was adapted from CalibrateTask.copyIcSourceFields().
1471 """
1472 control = afwTable.MatchControl()
1473 matchRadius = 3.0 if psfSigma is None else max(3.0, psfSigma) # in pixels
1474 # Return all matched objects, to separate blends.
1475 control.findOnlyClosest = False
1476 matches = afwTable.matchXy(psf_stars, stars, matchRadius, control)
1477 deblend_key = stars.schema["deblend_nChild"].asKey()
1478 matches = [m for m in matches if m[1].get(deblend_key) == 0]
1480 # Because we had to allow multiple matches to handle parents, we now
1481 # need to prune to the best (closest) matches.
1482 # Closest matches is a dict of psf_stars source ID to Match record
1483 # (psf_stars source, sourceCat source, distance in pixels).
1484 best = {}
1485 for match_psf, match_stars, d in matches:
1486 match = best.get(match_psf.getId())
1487 if match is None or d <= match[2]:
1488 best[match_psf.getId()] = (match_psf, match_stars, d)
1489 matches = list(best.values())
1490 # We'll use this to construct index arrays into each catalog.
1491 ids = np.array([(match_psf.getId(), match_stars.getId()) for match_psf, match_stars, d in matches]).T
1493 if (n_matches := len(matches)) == 0:
1494 raise NoPsfStarsToStarsMatchError(n_psf_stars=len(psf_stars), n_stars=len(stars))
1496 self.log.info("%d psf/astrometry stars out of %d matched %d calib stars",
1497 n_matches, len(psf_stars), len(stars))
1498 self.metadata["matched_psf_star_count"] = n_matches
1500 # Check that no stars sources are listed twice; we already know
1501 # that each match has a unique psf_stars id, due to using as the key
1502 # in best above.
1503 n_unique = len(set(m[1].getId() for m in matches))
1504 if n_unique != n_matches:
1505 self.log.warning("%d psf_stars matched only %d stars", n_matches, n_unique)
1507 # The indices of the IDs, so we can update the flag fields as arrays.
1508 idx_psf_stars = np.searchsorted(psf_stars["id"], ids[0])
1509 idx_stars = np.searchsorted(stars["id"], ids[1])
1510 for field in self.psf_fields:
1511 result = np.zeros(len(stars), dtype=bool)
1512 result[idx_stars] = psf_stars[field][idx_psf_stars]
1513 stars[field] = result
1514 stars['psf_id'][idx_stars] = psf_stars['id'][idx_psf_stars]
1515 stars['psf_max_value'][idx_stars] = psf_stars['psf_max_value'][idx_psf_stars]
1517 def _fit_astrometry(self, exposure, stars):
1518 """Fit an astrometric model to the data and return the reference
1519 matches used in the fit, and the fitted WCS.
1521 Parameters
1522 ----------
1523 exposure : `lsst.afw.image.Exposure`
1524 Exposure that is being fit, to get PSF and other metadata from.
1525 Modified to add the fitted skyWcs.
1526 stars : `SourceCatalog`
1527 Good stars selected for use in calibration, with RA/Dec coordinates
1528 computed from the pixel positions and fitted WCS.
1530 Returns
1531 -------
1532 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1533 Reference/stars matches used in the fit.
1534 """
1535 result = self.astrometry.run(stars, exposure)
1536 return result.matches, result.matchMeta
1538 def _fit_photometry(self, exposure, stars):
1539 """Fit a photometric model to the data and return the reference
1540 matches used in the fit, and the fitted PhotoCalib.
1542 Parameters
1543 ----------
1544 exposure : `lsst.afw.image.Exposure`
1545 Exposure that is being fit, to get PSF and other metadata from.
1546 Has the fit `lsst.afw.image.PhotoCalib` attached, with pixel values
1547 unchanged.
1548 stars : `lsst.afw.table.SourceCatalog`
1549 Good stars selected for use in calibration.
1550 background : `lsst.afw.math.BackgroundList`
1551 Background that was fit to the exposure during detection of the
1552 above stars.
1554 Returns
1555 -------
1556 calibrated_stars : `lsst.afw.table.SourceCatalog`
1557 Star catalog with flux/magnitude columns computed from the fitted
1558 photoCalib (instFlux columns are retained as well).
1559 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1560 Reference/stars matches used in the fit.
1561 matchMeta : `lsst.daf.base.PropertyList`
1562 Metadata needed to unpersist matches, as returned by the matcher.
1563 photo_calib : `lsst.afw.image.PhotoCalib`
1564 Photometric calibration that was fit to the star catalog.
1565 """
1566 result = self.photometry.run(exposure, stars)
1567 calibrated_stars = result.photoCalib.calibrateCatalog(stars)
1568 exposure.setPhotoCalib(result.photoCalib)
1569 return calibrated_stars, result.matches, result.matchMeta, result.photoCalib
1571 def _apply_photometry(self, exposure, background, background_to_photometric_ratio=None):
1572 """Apply the photometric model attached to the exposure to the
1573 exposure's pixels and an associated background model.
1575 Parameters
1576 ----------
1577 exposure : `lsst.afw.image.Exposure`
1578 Exposure with the target `lsst.afw.image.PhotoCalib` attached.
1579 On return, pixel values will be calibrated and an identity
1580 photometric transform will be attached.
1581 background : `lsst.afw.math.BackgroundList`
1582 Background model to convert to nanojansky units in place.
1583 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1584 Image to convert photometric-flattened image to
1585 background-flattened image.
1586 """
1587 photo_calib = exposure.getPhotoCalib()
1588 exposure.maskedImage = photo_calib.calibrateImage(exposure.maskedImage)
1589 identity = afwImage.PhotoCalib(1.0,
1590 photo_calib.getCalibrationErr(),
1591 bbox=exposure.getBBox())
1592 exposure.setPhotoCalib(identity)
1593 exposure.metadata["BUNIT"] = "nJy"
1595 assert photo_calib._isConstant, \
1596 "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that."
1598 for bg in background:
1599 # The statsImage is a view, but we can't assign to a function call
1600 # in python.
1601 binned_image = bg[0].getStatsImage()
1602 binned_image *= photo_calib.getCalibrationMean()
1604 def _summarize(self, exposure, stars, background):
1605 """Compute summary statistics on the exposure and update in-place the
1606 calibrations attached to it.
1608 Parameters
1609 ----------
1610 exposure : `lsst.afw.image.Exposure`
1611 Exposure that was calibrated, to get PSF and other metadata from.
1612 Should be in instrumental units with the photometric calibration
1613 attached.
1614 Modified to contain the computed summary statistics.
1615 stars : `SourceCatalog`
1616 Good stars selected used in calibration.
1617 background : `lsst.afw.math.BackgroundList`
1618 Background that was fit to the exposure during detection of the
1619 above stars. Should be in instrumental units.
1620 """
1621 summary = self.compute_summary_stats.run(exposure, stars, background)
1622 exposure.info.setSummaryStats(summary)
1624 def _recordMaskedPixelFractions(self, exposure):
1625 """Record the fraction of all the pixels in an exposure
1626 that are masked with a given flag. Each fraction is
1627 recorded in the task metadata. One record per flag type.
1629 Parameters
1630 ----------
1631 exposure : `lsst.afw.image.ExposureF`
1632 The target exposure to calculate masked pixel fractions for.
1633 """
1635 mask = exposure.mask
1636 maskPlanes = list(mask.getMaskPlaneDict().keys())
1637 for maskPlane in maskPlanes:
1638 self.metadata[f"{maskPlane.lower()}_mask_fraction"] = (
1639 evaluateMaskFraction(mask, maskPlane)
1640 )
1642 def _update_wcs_with_camera_model(self, exposure, cameraModel):
1643 """Replace the existing WCS with one generated using the pointing
1644 and the input camera model.
1646 Parameters
1647 ----------
1648 exposure : `lsst.afw.image.ExposureF`
1649 The exposure whose WCS will be updated.
1650 cameraModel : `lsst.afw.cameraGeom.Camera`
1651 Camera to be used when constructing updated WCS.
1652 """
1653 detector = cameraModel[exposure.detector.getId()]
1654 newWcs = createInitialSkyWcs(exposure.visitInfo, detector)
1655 exposure.setWcs(newWcs)
1657 def _compute_mask_fraction(self, mask, mask_planes, bad_mask_planes):
1658 """Evaluate the fraction of masked pixels in a (set of) mask plane(s).
1660 Parameters
1661 ----------
1662 mask : `lsst.afw.image.Mask`
1663 The mask on which to evaluate the fraction.
1664 mask_planes : `list`, `str`
1665 The mask planes for which to evaluate the fraction.
1666 bad_mask_planes : `list`, `str`
1667 The mask planes to exclude from the computation.
1669 Returns
1670 -------
1671 detected_fraction : `float`
1672 The calculated fraction of masked pixels
1673 """
1674 bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes)
1675 n_good_pix = np.sum(mask.array & bad_pixel_mask == 0)
1676 if n_good_pix == 0:
1677 detected_fraction = float("nan")
1678 return detected_fraction
1679 detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes)
1680 n_detected_pix = np.sum((mask.array & detected_pixel_mask != 0)
1681 & (mask.array & bad_pixel_mask == 0))
1682 detected_fraction = n_detected_pix/n_good_pix
1683 return detected_fraction
1685 def _compute_per_amp_fraction(self, exposure, detected_fraction, mask_planes, bad_mask_planes):
1686 """Evaluate the maximum per-amplifier fraction of masked pixels.
1688 Parameters
1689 ----------
1690 exposure : `lsst.afw.image.ExposureF`
1691 The exposure on which to compute the per-amp masked fraction.
1692 detected_fraction : `float`
1693 The current detected_fraction of the ``mask_planes`` for the
1694 full detector.
1695 mask_planes : `list`, `str`
1696 The mask planes for which to evaluate the fraction.
1697 bad_mask_planes : `list`, `str`
1698 The mask planes to exclude from the computation.
1700 Returns
1701 -------
1702 n_above_max_per_amp : `int`
1703 The number of amplifiers with masked fractions above a maximum
1704 value (set by the current full-detector ``detected_fraction``).
1705 highest_detected_fraction_per_amp : `float`
1706 The highest value of the per-amplifier fraction of masked pixels.
1707 no_zero_det_amps : `bool`
1708 A boolean representing whether any of the amplifiers has zero
1709 masked pixels.
1710 """
1711 highest_detected_fraction_per_amp = -9.99
1712 n_above_max_per_amp = 0
1713 n_no_zero_det_amps = 0
1714 no_zero_det_amps = True
1715 amps = exposure.detector.getAmplifiers()
1716 if amps is not None:
1717 for ia, amp in enumerate(amps):
1718 amp_bbox = amp.getBBox()
1719 exp_bbox = exposure.getBBox()
1720 if not exp_bbox.contains(amp_bbox):
1721 self.log.info("Bounding box of amplifier (%s) does not fit in exposure's "
1722 "bounding box (%s). Skipping...", amp_bbox, exp_bbox)
1723 continue
1724 sub_image = exposure.subset(amp.getBBox())
1725 detected_fraction_amp = self._compute_mask_fraction(sub_image.mask,
1726 mask_planes,
1727 bad_mask_planes)
1728 self.log.debug("Current detected fraction for amplifier %s = %.3f",
1729 amp.getName(), detected_fraction_amp)
1730 if detected_fraction_amp < 0.002:
1731 n_no_zero_det_amps += 1
1732 if n_no_zero_det_amps > 2:
1733 no_zero_det_amps = False
1734 break
1735 highest_detected_fraction_per_amp = max(detected_fraction_amp,
1736 highest_detected_fraction_per_amp)
1737 if highest_detected_fraction_per_amp > min(0.998, max(0.8, 3.0*detected_fraction)):
1738 n_above_max_per_amp += 1
1739 if n_above_max_per_amp > 2:
1740 break
1741 else:
1742 self.log.info("No amplifier object for detector %d, so skipping per-amp "
1743 "detection fraction checks.", exposure.detector.getId())
1744 return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps
1746 def _remeasure_star_background(self, result, background_to_photometric_ratio=None):
1747 """Remeasure the exposure's background with iterative adaptive
1748 threshold detection.
1750 Parameters
1751 ----------
1752 result : `lsst.pipe.base.Struct`
1753 The current state of the result Strut from the run method of
1754 CalibrateImageTask. Will be modified in place.
1755 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1756 Image to convert photometric-flattened image to
1757 background-flattened image.
1759 Returns
1760 -------
1761 result : `lsst.pipe.base.Struct`
1762 The modified result Struct with the new background subtracted.
1763 """
1764 # Restore the previously measured backgroud and remeasure it
1765 # using an adaptive threshold detection iteration to ensure a
1766 # "Goldilocks Zone" for the fraction of detected pixels.
1767 backgroundOrig = result.background.clone()
1768 median_background = np.median(backgroundOrig.getImage().array)
1769 self.log.warning("Original median_background = %.2f", median_background)
1770 result.exposure.image.array += result.background.getImage().array
1771 result.background = afwMath.BackgroundList()
1773 origMask = result.exposure.mask.clone()
1774 bad_mask_planes = ["BAD", "EDGE", "NO_DATA"]
1775 detected_mask_planes = ["DETECTED", "DETECTED_NEGATIVE"]
1776 detected_fraction_orig = self._compute_mask_fraction(result.exposure.mask,
1777 detected_mask_planes,
1778 bad_mask_planes)
1779 minDetFracForFinalBg = 0.02
1780 maxDetFracForFinalBg = 0.93
1781 # Dilate the current detected mask planes and don't clear
1782 # them in the detection step.
1783 nPixToDilate = 10
1784 maxAdaptiveDetIter = 10
1785 for dilateIter in range(maxAdaptiveDetIter):
1786 dilatedMask = origMask.clone()
1787 for maskName in detected_mask_planes:
1788 # Compute the grown detection mask plane using SpanSet
1789 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)
1790 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)
1791 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)
1792 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())
1793 # Clear the detected mask plane
1794 detectedMask = dilatedMask.getMaskPlane(maskName)
1795 dilatedMask.clearMaskPlane(detectedMask)
1796 # Set the mask plane to the dilated one
1797 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)
1799 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,
1800 detected_mask_planes,
1801 bad_mask_planes)
1802 if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1:
1803 break
1804 else:
1805 nPixToDilate -= 1
1807 result.exposure.mask = dilatedMask
1808 self.log.warning("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f",
1809 detected_fraction_orig, detected_fraction_dilated)
1810 n_above_max_per_amp = -99
1811 highest_detected_fraction_per_amp = float("nan")
1813 # Check the per-amplifier detection fractions.
1814 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \
1815 self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated,
1816 detected_mask_planes, bad_mask_planes)
1817 self.log.warning("Dilated mask: n_above_max_per_amp = %d, "
1818 "highest_detected_fraction_per_amp = %.3f",
1819 n_above_max_per_amp, highest_detected_fraction_per_amp)
1821 detected_fraction = 1.0
1822 nFootprintTemp = 1e12
1823 starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
1824 starBackgroundDetectionConfig.doTempLocalBackground = False
1825 starBackgroundDetectionConfig.nSigmaToGrow = 70.0
1826 starBackgroundDetectionConfig.reEstimateBackground = False
1827 starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0
1828 starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background)
1829 starBackgroundDetectionConfig.thresholdType = "pixel_stdev" # "stdev"
1831 n_above_max_per_amp = -99
1832 highest_detected_fraction_per_amp = float("nan")
1833 doCheckPerAmpDetFraction = True
1835 maxIter = 40
1836 for nIter in range(maxIter):
1837 currentThresh = starBackgroundDetectionConfig.thresholdValue
1838 if detected_fraction > maxDetFracForFinalBg:
1839 starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh
1840 if nFootprintTemp < 3 and detected_fraction > 0.9*maxDetFracForFinalBg:
1841 starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh
1842 if n_above_max_per_amp > 1:
1843 starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh
1844 if detected_fraction < minDetFracForFinalBg:
1845 starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh
1846 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
1847 config=starBackgroundDetectionConfig)
1848 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)
1849 tempDetections = starBackgroundDetectionTask.run(
1850 table=table, exposure=result.exposure, clearMask=True)
1851 result.exposure.mask |= dilatedMask
1852 nFootprintTemp = len(tempDetections.sources)
1853 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
1854 bad_mask_planes)
1855 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
1856 "DETECTED_NEGATIVE in star_background_detection = %.3f "
1857 "(max is %.3f; min is %.3f)",
1858 nIter, starBackgroundDetectionConfig.thresholdValue,
1859 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg)
1861 n_amp = len(result.exposure.detector.getAmplifiers())
1862 if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg:
1863 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \
1864 self._compute_per_amp_fraction(result.exposure, detected_fraction,
1865 detected_mask_planes, bad_mask_planes)
1867 if not no_zero_det_amps:
1868 starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh
1870 if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg
1871 and n_above_max_per_amp < int(0.75*n_amp)
1872 and no_zero_det_amps):
1873 if (n_above_max_per_amp < max(1, int(0.15*n_amp))
1874 or detected_fraction < 0.85*maxDetFracForFinalBg):
1875 break
1876 else:
1877 self.log.warning("Making small tweak....")
1878 starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh
1879 self.log.warning("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp))
1881 self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f "
1882 "(highest per amp section = %.5f)",
1883 detected_fraction, highest_detected_fraction_per_amp)
1885 if detected_fraction > maxDetFracForFinalBg:
1886 result.exposure.mask = dilatedMask
1887 self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE "
1888 "was too large in star_background_detection = %.3f (max = %.3f). "
1889 "Reverting to dilated mask from PSF detection...",
1890 detected_fraction, maxDetFracForFinalBg)
1891 star_background = self.star_background.run(
1892 exposure=result.exposure,
1893 backgroundToPhotometricRatio=background_to_photometric_ratio,
1894 ).background
1895 result.background.append(star_background[0])
1897 # Perform one more round of background subtraction that is just an
1898 # overall pedestal (order = 0). This is intended to account for
1899 # any potential gross oversubtraction imposed by the higher-order
1900 # subtraction.
1901 # Dilate DETECTED mask a bit more if it's below 50% detected.
1902 nPixToDilate = 2
1903 if detected_fraction < 0.5:
1904 dilatedMask = result.exposure.mask.clone()
1905 for maskName in detected_mask_planes:
1906 # Compute the grown detection mask plane using SpanSet
1907 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)
1908 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)
1909 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)
1910 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())
1911 # Clear the detected mask plane
1912 detectedMask = dilatedMask.getMaskPlane(maskName)
1913 dilatedMask.clearMaskPlane(detectedMask)
1914 # Set the mask plane to the dilated one
1915 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)
1917 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,
1918 detected_mask_planes,
1919 bad_mask_planes)
1920 result.exposure.mask = dilatedMask
1921 self.log.debug("detected_fraction_orig = %.3f detected_fraction_dilated = %.3f",
1922 detected_fraction_orig, detected_fraction_dilated)
1924 pedestalBackgroundConfig = lsst.meas.algorithms.SubtractBackgroundConfig()
1925 pedestalBackgroundConfig.statisticsProperty = "MEDIAN"
1926 pedestalBackgroundConfig.approxOrderX = 0
1927 pedestalBackgroundConfig.binSize = 64
1928 pedestalBackgroundTask = lsst.meas.algorithms.SubtractBackgroundTask(config=pedestalBackgroundConfig)
1929 pedestalBackgroundList = pedestalBackgroundTask.run(
1930 exposure=result.exposure,
1931 background=result.background,
1932 backgroundToPhotometricRatio=background_to_photometric_ratio,
1933 ).background
1934 # Isolate the final pedestal background to log the computed value
1935 pedestalBackground = afwMath.BackgroundList()
1936 pedestalBackground.append(pedestalBackgroundList[1])
1937 pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0]
1938 self.log.warning("Subtracted pedestalBackgroundLevel = %.4f", pedestalBackgroundLevel)
1940 # Clear detected mask plane before final round of detection
1941 mask = result.exposure.mask
1942 for mp in detected_mask_planes:
1943 if mp not in mask.getMaskPlaneDict():
1944 mask.addMaskPlane(mp)
1945 mask &= ~mask.getPlaneBitMask(detected_mask_planes)
1947 return result