Coverage for python / lsst / pipe / tasks / calibrateImage.py: 12%
835 statements
« prev ^ index » next coverage.py v7.14.0, created at 2026-05-13 08:48 +0000
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 math
26import numpy as np
27import requests
28import os
30from lsst.geom import SpherePoint, degrees
31from lsst.afw.geom import SpanSet
32import lsst.afw.table as afwTable
33import lsst.afw.image as afwImage
34import lsst.afw.math as afwMath
35from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache
36import lsst.meas.algorithms
37import lsst.meas.algorithms.installGaussianPsf
38import lsst.meas.algorithms.measureApCorr
39import lsst.meas.algorithms.setPrimaryFlags
40from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask
41from lsst.meas.algorithms.computeRoughPsfShapelets import ComputeRoughPsfShapeletsTask
42import lsst.meas.base
43import lsst.meas.astrom
44import lsst.meas.deblender
45import lsst.meas.extensions.shapeHSM
46from lsst.obs.base.utils import createInitialSkyWcsFromBoresight
47import lsst.pex.config as pexConfig
48import lsst.pipe.base as pipeBase
49from lsst.pipe.base import connectionTypes
50from lsst.utils.timer import timeMethod
52from . import measurePsf, repair, photoCal, computeExposureSummaryStats, snapCombine, diffractionSpikeMask
55class AllCentroidsFlaggedError(pipeBase.AlgorithmError):
56 """Raised when there are no valid centroids during psf fitting.
57 """
58 def __init__(self, shapelets_iq_score, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size):
59 msg = (f"All source centroids (out of {n_sources}) flagged during PSF fitting. "
60 "Original image PSF is likely unusable; best-fit PSF shape parameters: "
61 f"Ixx={psf_shape_ixx}, Iyy={psf_shape_iyy}, Ixy={psf_shape_ixy}, size={psf_size}"
62 )
63 super().__init__(msg)
64 self.shapelets_iq_score = shapelets_iq_score
65 self.n_sources = n_sources
66 self.psf_shape_ixx = psf_shape_ixx
67 self.psf_shape_iyy = psf_shape_iyy
68 self.psf_shape_ixy = psf_shape_ixy
69 self.psf_size = psf_size
71 @property
72 def metadata(self):
73 return {"shapelets_iq_score": self.shapelets_iq_score,
74 "n_sources": self.n_sources,
75 "psf_shape_ixx": self.psf_shape_ixx,
76 "psf_shape_iyy": self.psf_shape_iyy,
77 "psf_shape_ixy": self.psf_shape_ixy,
78 "psf_size": self.psf_size
79 }
82class NoPsfStarsToStarsMatchError(pipeBase.AlgorithmError):
83 """Raised when there are no matches between the psf_stars and stars
84 catalogs.
85 """
86 def __init__(self, *, n_psf_stars, n_stars):
87 msg = (f"No psf stars out of {n_psf_stars} matched {n_stars} calib stars."
88 " Downstream processes probably won't have useful {calib_type} stars in this case."
89 " Is `star_selector` too strict or is this a bad image?")
90 super().__init__(msg)
91 self.n_psf_stars = n_psf_stars
92 self.n_stars = n_stars
94 @property
95 def metadata(self):
96 return {"n_psf_stars": self.n_psf_stars,
97 "n_stars": self.n_stars
98 }
101class CalibrateImageConnections(pipeBase.PipelineTaskConnections,
102 dimensions=("instrument", "visit", "detector")):
104 astrometry_ref_cat = connectionTypes.PrerequisiteInput(
105 doc="Reference catalog to use for astrometric calibration.",
106 name="gaia_dr3_20230707",
107 storageClass="SimpleCatalog",
108 dimensions=("skypix",),
109 deferLoad=True,
110 multiple=True,
111 )
112 photometry_ref_cat = connectionTypes.PrerequisiteInput(
113 doc="Reference catalog to use for photometric calibration.",
114 name="ps1_pv3_3pi_20170110",
115 storageClass="SimpleCatalog",
116 dimensions=("skypix",),
117 deferLoad=True,
118 multiple=True
119 )
120 camera_model = connectionTypes.PrerequisiteInput(
121 doc="Camera distortion model to use for astrometric calibration.",
122 name="astrometry_camera",
123 storageClass="Camera",
124 dimensions=("instrument", "physical_filter"),
125 isCalibration=True,
126 minimum=0, # Can fall back on WCS already attached to exposure.
127 )
128 exposures = connectionTypes.Input(
129 doc="Exposure (or two snaps) to be calibrated, and detected and measured on.",
130 name="postISRCCD",
131 storageClass="Exposure",
132 multiple=True, # to handle 1 exposure or 2 snaps
133 dimensions=["instrument", "exposure", "detector"],
134 )
136 background_flat = connectionTypes.PrerequisiteInput(
137 name="flat",
138 doc="Flat calibration frame used for background correction.",
139 storageClass="ExposureF",
140 dimensions=["instrument", "detector", "physical_filter"],
141 isCalibration=True,
142 )
143 illumination_correction = connectionTypes.PrerequisiteInput(
144 name="illuminationCorrection",
145 doc="Illumination correction frame.",
146 storageClass="Exposure",
147 dimensions=["instrument", "detector", "physical_filter"],
148 isCalibration=True,
149 )
151 # outputs
152 initial_stars_schema = connectionTypes.InitOutput(
153 doc="Schema of the output initial stars catalog.",
154 name="initial_stars_schema",
155 storageClass="SourceCatalog",
156 )
158 # TODO DM-38732: We want some kind of flag on Exposures/Catalogs to make
159 # it obvious which components had failed to be computed/persisted.
160 exposure = connectionTypes.Output(
161 doc="Photometrically calibrated, background-subtracted exposure with fitted calibrations and "
162 "summary statistics. To recover the original exposure, first add the background "
163 "(`initial_pvi_background`), and then uncalibrate (divide by `initial_photoCalib_detector`).",
164 name="initial_pvi",
165 storageClass="ExposureF",
166 dimensions=("instrument", "visit", "detector"),
167 )
168 stars = connectionTypes.Output(
169 doc="Catalog of unresolved sources detected on the calibrated exposure.",
170 name="initial_stars_detector",
171 storageClass="ArrowAstropy",
172 dimensions=["instrument", "visit", "detector"],
173 )
174 stars_footprints = connectionTypes.Output(
175 doc="Catalog of unresolved sources detected on the calibrated exposure; "
176 "includes source footprints.",
177 name="initial_stars_footprints_detector",
178 storageClass="SourceCatalog",
179 dimensions=["instrument", "visit", "detector"],
180 )
181 applied_photo_calib = connectionTypes.Output(
182 doc=(
183 "Photometric calibration that was applied to exposure's pixels. "
184 "This connection is disabled when do_calibrate_pixels=False."
185 ),
186 name="initial_photoCalib_detector",
187 storageClass="PhotoCalib",
188 dimensions=("instrument", "visit", "detector"),
189 )
190 background = connectionTypes.Output(
191 doc="Background models estimated during calibration task; calibrated to be in nJy units.",
192 name="initial_pvi_background",
193 storageClass="Background",
194 dimensions=("instrument", "visit", "detector"),
195 )
196 background_to_photometric_ratio = connectionTypes.Output(
197 doc="Ratio of a background-flattened image to a photometric-flattened image. Only persisted "
198 "if do_illumination_correction is True.",
199 name="background_to_photometric_ratio",
200 storageClass="Image",
201 dimensions=("instrument", "visit", "detector"),
202 )
204 # Optional outputs
205 mask = connectionTypes.Output(
206 doc="The mask plane of the calibrated exposure, written as a separate "
207 "output so that it can be passed to downstream tasks even if the "
208 "exposure is removed to save storage.",
209 name="preliminary_visit_mask",
210 storageClass="Mask",
211 dimensions=("instrument", "visit", "detector"),
212 )
213 psf_stars_footprints = connectionTypes.Output(
214 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination; "
215 "includes source footprints.",
216 name="initial_psf_stars_footprints_detector",
217 storageClass="SourceCatalog",
218 dimensions=["instrument", "visit", "detector"],
219 )
220 psf_stars = connectionTypes.Output(
221 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination.",
222 name="initial_psf_stars_detector",
223 storageClass="ArrowAstropy",
224 dimensions=["instrument", "visit", "detector"],
225 )
226 astrometry_matches = connectionTypes.Output(
227 doc="Source to reference catalog matches from the astrometry solver.",
228 name="initial_astrometry_match_detector",
229 storageClass="Catalog",
230 dimensions=("instrument", "visit", "detector"),
231 )
232 photometry_matches = connectionTypes.Output(
233 doc="Source to reference catalog matches from the photometry solver.",
234 name="initial_photometry_match_detector",
235 storageClass="Catalog",
236 dimensions=("instrument", "visit", "detector"),
237 )
239 def __init__(self, *, config=None):
240 super().__init__(config=config)
241 if "psf_stars" not in config.optional_outputs:
242 del self.psf_stars
243 if "psf_stars_footprints" not in config.optional_outputs:
244 del self.psf_stars_footprints
245 if "astrometry_matches" not in config.optional_outputs:
246 del self.astrometry_matches
247 if "photometry_matches" not in config.optional_outputs:
248 del self.photometry_matches
249 if "mask" not in config.optional_outputs:
250 del self.mask
251 if not config.do_calibrate_pixels:
252 del self.applied_photo_calib
253 if not config.do_illumination_correction:
254 del self.background_flat
255 del self.illumination_correction
256 del self.background_to_photometric_ratio
257 if not config.useButlerCamera:
258 del self.camera_model
261class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CalibrateImageConnections):
262 optional_outputs = pexConfig.ListField(
263 doc="Which optional outputs to save (as their connection name)?",
264 dtype=str,
265 # TODO: note somewhere to disable this for benchmarking, but should
266 # we always have it on for production runs?
267 default=["psf_stars", "psf_stars_footprints", "astrometry_matches", "photometry_matches", "mask"],
268 optional=False
269 )
271 # To generate catalog ids consistently across subtasks.
272 id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()
274 snap_combine = pexConfig.ConfigurableField(
275 target=snapCombine.SnapCombineTask,
276 doc="Task to combine two snaps to make one exposure.",
277 )
279 # subtasks used during psf characterization
280 install_simple_psf = pexConfig.ConfigurableField(
281 target=lsst.meas.algorithms.installGaussianPsf.InstallGaussianPsfTask,
282 doc="Task to install a simple PSF model into the input exposure to use "
283 "when detecting bright sources for PSF estimation.",
284 )
285 compute_shapelet_decomposition = pexConfig.ConfigurableField(
286 target=ComputeRoughPsfShapeletsTask,
287 doc="Task to compute a rough shapelet expansion of the PSF from a set "
288 "of high S/N detections.",
289 )
290 do_compute_shapelet_decomposition = pexConfig.Field(
291 dtype=bool,
292 default=True,
293 doc="Compute a rough shapelet expansion of the PSF stars?",
294 )
295 psf_repair = pexConfig.ConfigurableField(
296 target=repair.RepairTask,
297 doc="Task to repair cosmic rays on the exposure before PSF determination.",
298 )
299 psf_subtract_background = pexConfig.ConfigurableField(
300 target=lsst.meas.algorithms.SubtractBackgroundTask,
301 doc="Task to perform initial background subtraction, before first detection pass.",
302 )
303 psf_detection = pexConfig.ConfigurableField(
304 target=lsst.meas.algorithms.SourceDetectionTask,
305 doc="Task to detect sources for PSF determination."
306 )
307 do_adaptive_threshold_detection = pexConfig.Field(
308 dtype=bool,
309 default=True,
310 doc="Implement the adaptive detection thresholding approach?",
311 )
312 do_remeasure_star_background = pexConfig.Field(
313 dtype=bool,
314 default=True,
315 doc="Do iterative star background measurement (ignored if do_adaptive_threshold_detection is False).",
316 )
317 psf_adaptive_threshold_detection = pexConfig.ConfigurableField(
318 target=AdaptiveThresholdDetectionTask,
319 doc="Task to adaptively detect sources for PSF determination.",
320 )
321 psf_source_measurement = pexConfig.ConfigurableField(
322 target=lsst.meas.base.SingleFrameMeasurementTask,
323 doc="Task to measure sources to be used for psf estimation."
324 )
325 psf_measure_psf = pexConfig.ConfigurableField(
326 target=measurePsf.MeasurePsfTask,
327 doc="Task to measure the psf on bright sources."
328 )
329 psf_normalized_calibration_flux = pexConfig.ConfigurableField(
330 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
331 doc="Task to normalize the calibration flux (e.g. compensated tophats) "
332 "for the bright stars used for psf estimation.",
333 )
335 # TODO DM-39203: we can remove aperture correction from this task once we
336 # are using the shape-based star/galaxy code.
337 measure_aperture_correction = pexConfig.ConfigurableField(
338 target=lsst.meas.algorithms.measureApCorr.MeasureApCorrTask,
339 doc="Task to compute the aperture correction from the bright stars."
340 )
342 # subtasks used during star measurement
343 star_background = pexConfig.ConfigurableField(
344 target=lsst.meas.algorithms.SubtractBackgroundTask,
345 doc="Task to perform final background subtraction, just before photoCal.",
346 )
347 star_background_min_footprints = pexConfig.Field(
348 dtype=int,
349 default=3,
350 doc="Minimum number of footprints in the detection mask for star_background measurement. "
351 "This number will get adjusted to the fraction config.star_background_peak_fraction "
352 "of the detected peaks if that number is larger. If the number of footprints is less "
353 "than the minimum, the detection threshold is iteratively increased until the "
354 "threshold is met.",
355 )
356 star_background_peak_fraction = pexConfig.Field(
357 dtype=float,
358 default=0.01,
359 doc="The minimum number of footprints in the detection mask for star_background measurement. "
360 "gets set to the maximum of this fraction of the detected peaks and the value set in "
361 "config.star_background_min_footprints. If the number of footprints is less than the "
362 "current minimum set, the detection threshold is iteratively increased until the "
363 "threshold is met.",
364 )
365 star_detection = pexConfig.ConfigurableField(
366 target=lsst.meas.algorithms.SourceDetectionTask,
367 doc="Task to detect stars to return in the output catalog."
368 )
369 star_sky_sources = pexConfig.ConfigurableField(
370 target=lsst.meas.algorithms.SkyObjectsTask,
371 doc="Task to generate sky sources ('empty' regions where there are no detections).",
372 )
373 do_downsample_footprints = pexConfig.Field(
374 dtype=bool,
375 default=False,
376 doc="Downsample footprints prior to deblending to optimize speed?",
377 )
378 downsample_max_footprints = pexConfig.Field(
379 dtype=int,
380 default=1000,
381 doc="Maximum number of non-sky-source footprints to use if do_downsample_footprints is True,",
382 )
383 star_deblend = pexConfig.ConfigurableField(
384 target=lsst.meas.deblender.SourceDeblendTask,
385 doc="Split blended sources into their components."
386 )
387 star_measurement = pexConfig.ConfigurableField(
388 target=lsst.meas.base.SingleFrameMeasurementTask,
389 doc="Task to measure stars to return in the output catalog."
390 )
391 star_normalized_calibration_flux = pexConfig.ConfigurableField(
392 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
393 doc="Task to apply the normalization for calibration fluxes (e.g. compensated tophats) "
394 "for the final output star catalog.",
395 )
396 star_apply_aperture_correction = pexConfig.ConfigurableField(
397 target=lsst.meas.base.ApplyApCorrTask,
398 doc="Task to apply aperture corrections to the selected stars."
399 )
400 star_catalog_calculation = pexConfig.ConfigurableField(
401 target=lsst.meas.base.CatalogCalculationTask,
402 doc="Task to compute extendedness values on the star catalog, "
403 "for the star selector to remove extended sources."
404 )
405 star_set_primary_flags = pexConfig.ConfigurableField(
406 target=lsst.meas.algorithms.setPrimaryFlags.SetPrimaryFlagsTask,
407 doc="Task to add isPrimary to the catalog."
408 )
409 star_selector = lsst.meas.algorithms.sourceSelectorRegistry.makeField(
410 default="science",
411 doc="Task to select reliable stars to use for calibration."
412 )
414 # final calibrations and statistics
415 astrometry = pexConfig.ConfigurableField(
416 target=lsst.meas.astrom.AstrometryTask,
417 doc="Task to perform astrometric calibration to fit a WCS.",
418 )
419 astrometry_ref_loader = pexConfig.ConfigField(
420 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
421 doc="Configuration of reference object loader for astrometric fit.",
422 )
423 photometry = pexConfig.ConfigurableField(
424 target=photoCal.PhotoCalTask,
425 doc="Task to perform photometric calibration to fit a PhotoCalib.",
426 )
427 photometry_ref_loader = pexConfig.ConfigField(
428 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
429 doc="Configuration of reference object loader for photometric fit.",
430 )
431 doMaskDiffractionSpikes = pexConfig.Field(
432 dtype=bool,
433 default=False,
434 doc="If True, load the photometric reference catalog again but select "
435 "only bright stars. Use the bright star catalog to set the SPIKE "
436 "mask for regions likely contaminated by diffraction spikes.",
437 )
438 diffractionSpikeMask = pexConfig.ConfigurableField(
439 target=diffractionSpikeMask.DiffractionSpikeMaskTask,
440 doc="Task to identify and mask the diffraction spikes of bright stars.",
441 )
443 compute_summary_stats = pexConfig.ConfigurableField(
444 target=computeExposureSummaryStats.ComputeExposureSummaryStatsTask,
445 doc="Task to to compute summary statistics on the calibrated exposure."
446 )
448 do_illumination_correction = pexConfig.Field(
449 dtype=bool,
450 default=False,
451 doc="If True, apply the illumination correction. This assumes that the "
452 "input image has already been flat-fielded such that it is suitable "
453 "for background subtraction.",
454 )
456 do_calibrate_pixels = pexConfig.Field(
457 dtype=bool,
458 default=True,
459 doc=(
460 "If True, apply the photometric calibration to the image pixels "
461 "and background model, and attach an identity PhotoCalib to "
462 "the output image to reflect this. If False`, leave the image "
463 "and background uncalibrated and attach the PhotoCalib that maps "
464 "them to physical units."
465 )
466 )
468 do_include_astrometric_errors = pexConfig.Field(
469 dtype=bool,
470 default=True,
471 doc="If True, include astrometric errors in the output catalog.",
472 )
474 background_stats_ignored_pixel_masks = pexConfig.ListField(
475 dtype=str,
476 default=["SAT", "SUSPECT", "SPIKE"],
477 doc="Pixel mask flags to ignore when calculating post-sky-subtraction "
478 "background statistics. These are added to those ignored by the "
479 "meas.algorithms.SubtractBackgroundConfig algorithm."
480 )
482 run_sattle = pexConfig.Field(
483 dtype=bool,
484 default=False,
485 doc="If True, the sattle service will populate a cache for later use "
486 "in ip_diffim.detectAndMeasure alert verification."
487 )
489 sattle_historical = pexConfig.Field(
490 dtype=bool,
491 default=False,
492 doc="If re-running a pipeline that requires sattle, this should be set "
493 "to True. This will populate sattle's cache with the historic data "
494 "closest in time to the exposure.",
495 )
497 useButlerExposureRegion = pexConfig.Field(
498 dtype=bool,
499 default=False,
500 doc="If True, use the exposure region stored in the butler as the "
501 "region for reference catalog trimming. If False, the reference loader "
502 "trimming will be set by the loader (typically from a padded version of "
503 "the detector's bounding box + current WCS)."
504 )
506 useButlerCamera = pexConfig.Field(
507 dtype=bool,
508 default=False,
509 doc="If True, use a camera distortion model generated elsewhere in "
510 "the pipeline combined with the telescope boresight as a starting "
511 "point for fitting the WCS, instead of using the WCS attached to "
512 "the exposure, which is generated from the boresight and the "
513 "camera model from the obs_* package."
514 )
516 background_stats_flux_column = pexConfig.Field(
517 dtype=str,
518 default="base_CircularApertureFlux_12_0_flux",
519 doc="Column used to generate post-subtracted background stats."
520 )
522 def setDefaults(self):
523 super().setDefaults()
525 # Use a very broad PSF here, to thoroughly reject CRs.
526 # TODO investigation: a large initial psf guess may make stars look
527 # like CRs for very good seeing images.
528 self.install_simple_psf.fwhm = 4
530 # S/N>=50 sources for PSF determination, but detection to S/N=10.
531 # The thresholdValue sets the minimum flux in a pixel to be included
532 # in the footprint, while peaks are only detected when they are above
533 # thresholdValue * includeThresholdMultiplier. The low thresholdValue
534 # ensures that the footprints are large enough for the noise replacer
535 # to mask out faint undetected neighbors that are not to be measured.
536 self.psf_detection.thresholdValue = 10.0
537 self.psf_detection.includeThresholdMultiplier = 5.0
538 # TODO investigation: Probably want False here, but that may require
539 # tweaking the background spatial scale, to make it small enough to
540 # prevent extra peaks in the wings of bright objects.
541 self.psf_detection.doTempLocalBackground = False
542 self.psf_detection.reEstimateBackground = False
544 # Minimal measurement plugins for PSF determination.
545 # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use
546 # shapeHSM/moments for star/galaxy separation.
547 # TODO DM-39203: we can remove aperture correction from this task
548 # once we are using the shape-based star/galaxy code.
549 self.psf_source_measurement.plugins = ["base_PixelFlags",
550 "base_SdssCentroid",
551 "ext_shapeHSM_HsmSourceMoments",
552 "base_CircularApertureFlux",
553 "base_GaussianFlux",
554 "base_PsfFlux",
555 "base_CompensatedTophatFlux",
556 ]
557 self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
558 # Only measure apertures we need for PSF measurement.
559 self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
560 self.psf_source_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
561 # TODO DM-40843: Remove this line once this is the psfex default.
562 self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \
563 "base_CircularApertureFlux_12_0_instFlux"
565 # No extendedness information available: we need the aperture
566 # corrections to determine that.
567 self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False
568 self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]
569 self.measure_aperture_correction.sourceSelector["science"].flags.bad = []
571 # Detection for good S/N for astrometry/photometry and other
572 # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks.
573 self.star_detection.reEstimateBackground = False
574 self.star_detection.thresholdValue = 5.0
575 self.star_detection.includeThresholdMultiplier = 2.0
576 self.star_measurement.plugins = ["base_PixelFlags",
577 "base_SdssCentroid",
578 "ext_shapeHSM_HsmSourceMoments",
579 "ext_shapeHSM_HsmPsfMoments",
580 "base_GaussianFlux",
581 "base_PsfFlux",
582 "base_CircularApertureFlux",
583 "base_ClassificationSizeExtendedness",
584 "base_CompensatedTophatFlux",
585 ]
586 self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
587 self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
588 # Only measure the apertures we need for star selection.
589 self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
590 self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
592 # We measure and apply the normalization aperture correction with
593 # the psf_normalized_calibration_flux task, and we only apply the
594 # normalization aperture correction for the full list of stars.
595 self.star_normalized_calibration_flux.do_measure_ap_corr = False
597 # Select stars with reliable measurements and no bad flags.
598 self.star_selector["science"].doFlags = True
599 self.star_selector["science"].doUnresolved = True
600 self.star_selector["science"].doSignalToNoise = True
601 self.star_selector["science"].signalToNoise.minimum = 10.0
602 # Keep sky sources in the output catalog, even though they aren't
603 # wanted for calibration.
604 self.star_selector["science"].doSkySources = True
605 # Set the flux and error fields
606 self.star_selector["science"].signalToNoise.fluxField = "slot_CalibFlux_instFlux"
607 self.star_selector["science"].signalToNoise.errField = "slot_CalibFlux_instFluxErr"
609 # Use the affine WCS fitter (assumes we have a good camera geometry).
610 self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)
611 # phot_g_mean is the primary Gaia band for all input bands.
612 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"
614 # Only reject sky sources; we already selected good stars.
615 self.astrometry.sourceSelector["science"].doFlags = True
616 self.astrometry.sourceSelector["science"].flags.good = [] # ["calib_psf_candidate"]
617 # self.astrometry.sourceSelector["science"].flags.bad = []
618 self.astrometry.sourceSelector["science"].doUnresolved = False
619 self.astrometry.sourceSelector["science"].doIsolated = False
620 self.astrometry.sourceSelector["science"].doRequirePrimary = False
621 self.photometry.match.sourceSelection.doFlags = True
622 self.photometry.match.sourceSelection.flags.bad = ["sky_source"]
623 # Unset the (otherwise reasonable, but we've already made the
624 # selections we want above) selection settings in PhotoCalTask.
625 self.photometry.match.sourceSelection.doRequirePrimary = False
626 self.photometry.match.sourceSelection.doUnresolved = False
628 def validate(self):
629 super().validate()
631 # Ensure that the normalization calibration flux tasks
632 # are configured correctly.
633 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:
634 msg = ("psf_normalized_calibration_flux task must be configured with do_measure_ap_corr=True "
635 "or else the normalization and calibration flux will not be properly measured.")
636 raise pexConfig.FieldValidationError(
637 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,
638 )
639 if self.star_normalized_calibration_flux.do_measure_ap_corr:
640 msg = ("star_normalized_calibration_flux task must be configured with do_measure_ap_corr=False "
641 "to apply the previously measured normalization to the full catalog of calibration "
642 "fluxes.")
643 raise pexConfig.FieldValidationError(
644 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,
645 )
647 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,
648 # because they'd be running too early to pick up the fitted PhotoCalib
649 # and WCS.
650 if "base_LocalWcs" in self.psf_source_measurement.plugins.names:
651 raise pexConfig.FieldValidationError(
652 CalibrateImageConfig.psf_source_measurement,
653 self,
654 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
655 )
656 if "base_LocalWcs" in self.star_measurement.plugins.names:
657 raise pexConfig.FieldValidationError(
658 CalibrateImageConfig.star_measurement,
659 self,
660 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
661 )
662 if "base_LocalPhotoCalib" in self.psf_source_measurement.plugins.names:
663 raise pexConfig.FieldValidationError(
664 CalibrateImageConfig.psf_source_measurement,
665 self,
666 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
667 "as it would be run before the photometry fit."
668 )
669 if "base_LocalPhotoCalib" in self.star_measurement.plugins.names:
670 raise pexConfig.FieldValidationError(
671 CalibrateImageConfig.star_measurement,
672 self,
673 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
674 "as it would be run before the photometry fit."
675 )
677 # Check for illumination correction and background consistency.
678 if self.do_illumination_correction:
679 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:
680 raise pexConfig.FieldValidationError(
681 CalibrateImageConfig.psf_subtract_background,
682 self,
683 "CalibrateImageTask.psf_subtract_background must be configured with "
684 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
685 )
686 if self.psf_detection.reEstimateBackground:
687 if not self.psf_detection.doApplyFlatBackgroundRatio:
688 raise pexConfig.FieldValidationError(
689 CalibrateImageConfig.psf_detection,
690 self,
691 "CalibrateImageTask.psf_detection background must be configured with "
692 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
693 )
694 if self.star_detection.reEstimateBackground:
695 if not self.star_detection.doApplyFlatBackgroundRatio:
696 raise pexConfig.FieldValidationError(
697 CalibrateImageConfig.star_detection,
698 self,
699 "CalibrateImageTask.star_detection background must be configured with "
700 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
701 )
702 if not self.star_background.doApplyFlatBackgroundRatio:
703 raise pexConfig.FieldValidationError(
704 CalibrateImageConfig.star_background,
705 self,
706 "CalibrateImageTask.star_background background must be configured with "
707 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
708 )
710 if self.run_sattle:
711 if not os.getenv("SATTLE_URI_BASE"):
712 raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self,
713 "Sattle requested but URI environment variable not set.")
715 if not self.do_adaptive_threshold_detection:
716 if not self.psf_detection.reEstimateBackground:
717 raise pexConfig.FieldValidationError(
718 CalibrateImageConfig.psf_detection,
719 self,
720 "If not using the adaptive threshold detection implementation (i.e. "
721 "config.do_adaptive_threshold_detection = False), CalibrateImageTask.psf_detection "
722 "background must be configured with "
723 "reEstimateBackground = True to maintain original behavior."
724 )
725 if not self.star_detection.reEstimateBackground:
726 raise pexConfig.FieldValidationError(
727 CalibrateImageConfig.psf_detection,
728 self,
729 "If not using the adaptive threshold detection implementation "
730 "(i.e. config.do_adaptive_threshold_detection = False), "
731 "CalibrateImageTask.star_detection background must be configured "
732 "with reEstimateBackground = True to maintain original behavior."
733 )
736class CalibrateImageTask(pipeBase.PipelineTask):
737 """Compute the PSF, aperture corrections, astrometric and photometric
738 calibrations, and summary statistics for a single science exposure, and
739 produce a catalog of brighter stars that were used to calibrate it.
741 Parameters
742 ----------
743 initial_stars_schema : `lsst.afw.table.Schema`
744 Schema of the initial_stars output catalog.
745 """
746 _DefaultName = "calibrateImage"
747 ConfigClass = CalibrateImageConfig
749 def __init__(self, initial_stars_schema=None, **kwargs):
750 super().__init__(**kwargs)
751 self.makeSubtask("snap_combine")
753 # PSF determination subtasks
754 self.makeSubtask("install_simple_psf")
755 self.makeSubtask("psf_repair")
756 self.makeSubtask("psf_subtract_background")
757 self.psf_schema = afwTable.SourceTable.makeMinimalSchema()
758 self.psf_schema.addField(
759 'psf_max_value',
760 type=np.float32,
761 doc="PSF max value.",
762 )
763 afwTable.CoordKey.addErrorFields(self.psf_schema)
764 if self.config.do_compute_shapelet_decomposition:
765 self.makeSubtask("compute_shapelet_decomposition", schema=self.psf_schema)
766 self.makeSubtask("psf_detection", schema=self.psf_schema)
767 self.makeSubtask("psf_source_measurement", schema=self.psf_schema)
768 self.makeSubtask("psf_adaptive_threshold_detection")
769 self.makeSubtask("psf_measure_psf", schema=self.psf_schema)
770 self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema)
772 self.makeSubtask("measure_aperture_correction", schema=self.psf_schema)
773 self.makeSubtask("astrometry", schema=self.psf_schema)
775 # star measurement subtasks
776 if initial_stars_schema is None:
777 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()
779 # These fields let us track which sources were used for psf modeling,
780 # astrometric fitting, and aperture correction calculations.
781 self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved",
782 "calib_astrometry_used",
783 # TODO DM-39203:
784 # these can be removed once apcorr is gone.
785 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",
786 "apcorr_base_PsfFlux_used",)
787 for field in self.psf_fields:
788 item = self.psf_schema.find(field)
789 initial_stars_schema.addField(item.getField())
790 id_type = self.psf_schema["id"].asField().getTypeString()
791 psf_max_value_type = self.psf_schema['psf_max_value'].asField().getTypeString()
792 initial_stars_schema.addField("psf_id",
793 type=id_type,
794 doc="id of this source in psf_stars; 0 if there is no match.")
795 initial_stars_schema.addField("psf_max_value",
796 type=psf_max_value_type,
797 doc="Maximum value in the star image used to train PSF.")
799 afwTable.CoordKey.addErrorFields(initial_stars_schema)
800 self.makeSubtask("star_background")
801 self.makeSubtask("star_detection", schema=initial_stars_schema)
802 self.makeSubtask("star_sky_sources", schema=initial_stars_schema)
803 self.makeSubtask("star_deblend", schema=initial_stars_schema)
804 self.makeSubtask("star_measurement", schema=initial_stars_schema)
805 self.makeSubtask("star_normalized_calibration_flux", schema=initial_stars_schema)
807 self.makeSubtask("star_apply_aperture_correction", schema=initial_stars_schema)
808 self.makeSubtask("star_catalog_calculation", schema=initial_stars_schema)
809 self.makeSubtask("star_set_primary_flags", schema=initial_stars_schema, isSingleFrame=True)
810 self.makeSubtask("star_selector")
811 self.makeSubtask("photometry", schema=initial_stars_schema)
812 if self.config.doMaskDiffractionSpikes:
813 self.makeSubtask("diffractionSpikeMask")
814 self.makeSubtask("compute_summary_stats")
816 # The final catalog will have calibrated flux columns, which we add to
817 # the init-output schema by calibrating our zero-length catalog with an
818 # arbitrary dummy PhotoCalib. We also use this schema to initialize
819 # the stars catalog in order to ensure it's the same even when we hit
820 # an error (and write partial outputs) before calibrating the catalog
821 # - note that calibrateCatalog will happily reuse existing output
822 # columns.
823 dummy_photo_calib = afwImage.PhotoCalib(1.0, 0, bbox=lsst.geom.Box2I())
824 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(
825 afwTable.SourceCatalog(initial_stars_schema)
826 )
828 def runQuantum(self, butlerQC, inputRefs, outputRefs):
829 inputs = butlerQC.get(inputRefs)
830 exposures = inputs.pop("exposures")
832 exposure_record = inputRefs.exposures[0].dataId.records["exposure"]
833 if self.config.useButlerExposureRegion:
834 exposure_region = butlerQC.quantum.dataId.region
835 else:
836 exposure_region = None
838 id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)
840 astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
841 dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],
842 refCats=inputs.pop("astrometry_ref_cat"),
843 name=self.config.connections.astrometry_ref_cat,
844 config=self.config.astrometry_ref_loader, log=self.log)
845 self.astrometry.setRefObjLoader(astrometry_loader)
847 photometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
848 dataIds=[ref.datasetRef.dataId for ref in inputRefs.photometry_ref_cat],
849 refCats=inputs.pop("photometry_ref_cat"),
850 name=self.config.connections.photometry_ref_cat,
851 config=self.config.photometry_ref_loader, log=self.log)
852 self.photometry.match.setRefObjLoader(photometry_loader)
854 if self.config.doMaskDiffractionSpikes:
855 # Use the same photometry reference catalog for the bright star
856 # mask.
857 self.diffractionSpikeMask.setRefObjLoader(photometry_loader)
859 if self.config.do_illumination_correction:
860 background_flat = inputs.pop("background_flat")
861 illumination_correction = inputs.pop("illumination_correction")
862 else:
863 background_flat = None
864 illumination_correction = None
866 camera_model = None
867 if self.config.useButlerCamera:
868 if "camera_model" in inputs:
869 camera_model = inputs.pop("camera_model")
870 else:
871 self.log.warning("useButlerCamera=True, but camera is not available for filter %s. The "
872 "astrometry fit will use the WCS already attached to the exposure.",
873 exposures[0].filter.bandLabel)
875 # This should not happen with a properly configured execution context.
876 assert not inputs, "runQuantum got more inputs than expected"
878 # Specify the fields that `annotate` needs below, to ensure they
879 # exist, even as None.
880 result = pipeBase.Struct(
881 exposure=None,
882 stars_footprints=None,
883 psf_stars_footprints=None,
884 background_to_photometric_ratio=None,
885 )
886 try:
887 self.run(
888 exposures=exposures,
889 result=result,
890 id_generator=id_generator,
891 background_flat=background_flat,
892 illumination_correction=illumination_correction,
893 camera_model=camera_model,
894 exposure_record=exposure_record,
895 exposure_region=exposure_region,
896 )
897 except pipeBase.AlgorithmError as e:
898 error = pipeBase.AnnotatedPartialOutputsError.annotate(
899 e,
900 self,
901 result.exposure,
902 result.psf_stars_footprints,
903 result.stars_footprints,
904 log=self.log
905 )
906 butlerQC.put(result, outputRefs)
907 raise error from e
909 butlerQC.put(result, outputRefs)
911 @timeMethod
912 def run(
913 self,
914 *,
915 exposures,
916 id_generator=None,
917 result=None,
918 background_flat=None,
919 illumination_correction=None,
920 camera_model=None,
921 exposure_record=None,
922 exposure_region=None,
923 ):
924 """Find stars and perform psf measurement, then do a deeper detection
925 and measurement and calibrate astrometry and photometry from that.
927 Parameters
928 ----------
929 exposures : `lsst.afw.image.Exposure` or \
930 `list` [`lsst.afw.image.Exposure`]
931 Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter.
932 Modified in-place during processing if only one is passed.
933 If two exposures are passed, treat them as snaps and combine
934 before doing further processing.
935 id_generator : `lsst.meas.base.IdGenerator`, optional
936 Object that generates source IDs and provides random seeds.
937 result : `lsst.pipe.base.Struct`, optional
938 Result struct that is modified to allow saving of partial outputs
939 for some failure conditions. If the task completes successfully,
940 this is also returned.
941 background_flat : `lsst.afw.image.Exposure`, optional
942 Background flat-field image.
943 illumination_correction : `lsst.afw.image.Exposure`, optional
944 Illumination correction image.
945 camera_model : `lsst.afw.cameraGeom.Camera`, optional
946 Camera to be used if constructing updated WCS.
947 exposure_record : `lsst.daf.butler.DimensionRecord`, optional
948 Butler metadata for the ``exposure`` dimension. Used if
949 constructing an updated WCS instead of the boresight and
950 orientation angle from the visit info.
951 exposure_region : `lsst.sphgeom.Region`, optional
952 The exposure region to use for the reference catalog filtering.
953 If `None`, this region will be set as a padded bbox + current WCS
954 of the exposure.
956 Returns
957 -------
958 result : `lsst.pipe.base.Struct`
959 Results as a struct with attributes:
961 ``exposure``
962 Calibrated exposure, with pixels in nJy units.
963 (`lsst.afw.image.Exposure`)
964 ``stars``
965 Stars that were used to calibrate the exposure, with
966 calibrated fluxes and magnitudes.
967 (`astropy.table.Table`)
968 ``stars_footprints``
969 Footprints of stars that were used to calibrate the exposure.
970 (`lsst.afw.table.SourceCatalog`)
971 ``psf_stars``
972 Stars that were used to determine the image PSF.
973 (`astropy.table.Table`)
974 ``psf_stars_footprints``
975 Footprints of stars that were used to determine the image PSF.
976 (`lsst.afw.table.SourceCatalog`)
977 ``background``
978 Background that was fit to the exposure when detecting
979 ``stars``. (`lsst.afw.math.BackgroundList`)
980 ``applied_photo_calib``
981 Photometric calibration that was fit to the star catalog and
982 applied to the exposure. (`lsst.afw.image.PhotoCalib`)
983 This is `None` if ``config.do_calibrate_pixels`` is `False`.
984 ``astrometry_matches``
985 Reference catalog stars matches used in the astrometric fit.
986 (`list` [`lsst.afw.table.ReferenceMatch`] or
987 `lsst.afw.table.BaseCatalog`).
988 ``photometry_matches``
989 Reference catalog stars matches used in the photometric fit.
990 (`list` [`lsst.afw.table.ReferenceMatch`] or
991 `lsst.afw.table.BaseCatalog`).
992 ``mask``
993 Copy of the mask plane of `exposure`.
994 (`lsst.afw.image.Mask`)
995 """
996 if result is None:
997 result = pipeBase.Struct()
998 if id_generator is None:
999 id_generator = lsst.meas.base.IdGenerator()
1001 result.exposure = self.snap_combine.run(exposures).exposure
1002 self.log.info("Initial PhotoCalib: %s", result.exposure.getPhotoCalib())
1004 result.exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = False
1006 # Check input image processing.
1007 if self.config.do_illumination_correction:
1008 if not result.exposure.metadata.get("LSST ISR FLAT APPLIED", False):
1009 raise pipeBase.InvalidQuantumError(
1010 "Cannot use do_illumination_correction with an image that has not had a flat applied",
1011 )
1013 # Update the exposure pointing to that stored in the butler record
1014 # (regardless of a camera model update). This is to pick up any
1015 # updates to the pointing stored in the butler record that may not
1016 # be reflected in the exposure metadata (headers).
1017 if camera_model:
1018 self._update_wcs_with_camera_model(
1019 result.exposure, camera_model, exposure_record=exposure_record
1020 )
1021 elif exposure_record is not None:
1022 self._update_wcs_to_exposure_record(result.exposure, exposure_record)
1024 result.background = None
1025 result.background_to_photometric_ratio = None
1026 summary_stat_catalog = None
1027 # Some exposure components are set to initial placeholder objects
1028 # while we try to bootstrap them. If we fail before we fit for them,
1029 # we want to reset those components to None so the placeholders don't
1030 # masquerade as the real thing.
1031 have_fit_psf = False
1032 have_fit_astrometry = False
1033 have_fit_photometry = False
1034 try:
1035 result.background_to_photometric_ratio = self._apply_illumination_correction(
1036 result.exposure,
1037 background_flat,
1038 illumination_correction,
1039 )
1040 compute_psf_struct = self._compute_psf(
1041 result,
1042 id_generator,
1043 )
1044 result.psf_stars_footprints = compute_psf_struct.detections_sources
1045 adaptive_det_res_struct = compute_psf_struct.adaptive_det_res_struct
1047 have_fit_psf = True
1049 if self.config.do_adaptive_threshold_detection:
1050 metadata_name = "psf_adaptive_threshold_value"
1051 result.exposure.metadata[metadata_name.upper()] = adaptive_det_res_struct.thresholdValue
1052 self.metadata[metadata_name] = adaptive_det_res_struct.thresholdValue
1053 metadata_name = "psf_adaptive_include_threshold_multiplier"
1054 result.exposure.metadata[metadata_name.upper()] = \
1055 adaptive_det_res_struct.includeThresholdMultiplier
1056 self.metadata[metadata_name] = adaptive_det_res_struct.includeThresholdMultiplier
1058 # Check if all centroids have been flagged. This should happen
1059 # after the PSF is fit and recorded because even a junky PSF
1060 # model is informative.
1061 if result.psf_stars_footprints["slot_Centroid_flag"].all():
1062 psf_shape = result.exposure.psf.computeShape(result.exposure.psf.getAveragePosition())
1063 raise AllCentroidsFlaggedError(
1064 shapelets_iq_score=result.exposure.info.getSummaryStats().shapeletsIqScore,
1065 n_sources=len(result.psf_stars_footprints),
1066 psf_shape_ixx=psf_shape.getIxx(),
1067 psf_shape_iyy=psf_shape.getIyy(),
1068 psf_shape_ixy=psf_shape.getIxy(),
1069 psf_size=psf_shape.getDeterminantRadius(),
1070 )
1072 if result.background is None:
1073 result.background = afwMath.BackgroundList()
1075 self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)
1076 result.psf_stars = result.psf_stars_footprints.asAstropy()
1077 # Run astrometry using PSF candidate stars.
1078 # Update "the psf_stars" source coordinates with the current wcs.
1079 afwTable.updateSourceCoords(
1080 result.exposure.wcs,
1081 sourceList=result.psf_stars_footprints,
1082 include_covariance=self.config.do_include_astrometric_errors,
1083 )
1084 astrometry_matches, astrometry_meta = self._fit_astrometry(
1085 result.exposure, result.psf_stars_footprints, exposure_region=exposure_region
1086 )
1087 num_astrometry_matches = len(astrometry_matches)
1088 if "astrometry_matches" in self.config.optional_outputs:
1089 result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches,
1090 astrometry_meta)
1091 result.psf_stars = result.psf_stars_footprints.asAstropy()
1093 # Add diffraction spike mask here so it can be used (i.e. avoided)
1094 # in the adaptive background estimation.
1095 if self.config.doMaskDiffractionSpikes:
1096 self.diffractionSpikeMask.run(result.exposure)
1098 self.metadata['adaptive_threshold_value'] = float("nan")
1099 if self.config.do_remeasure_star_background and self.config.do_adaptive_threshold_detection:
1100 self._remeasure_star_background(
1101 result,
1102 background_to_photometric_ratio=result.background_to_photometric_ratio,
1103 )
1105 # Run the stars_detection subtask for the photometric calibration.
1106 result.stars_footprints = self._find_stars(
1107 result.exposure,
1108 result.background,
1109 id_generator,
1110 background_to_photometric_ratio=result.background_to_photometric_ratio,
1111 adaptive_det_res_struct=adaptive_det_res_struct,
1112 num_astrometry_matches=num_astrometry_matches,
1113 )
1114 psf = result.exposure.getPsf()
1115 psfSigma = psf.computeShape(result.exposure.getBBox().getCenter()).getDeterminantRadius()
1116 self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints,
1117 psfSigma=psfSigma)
1119 # Update the "stars" source coordinates with the current wcs.
1120 afwTable.updateSourceCoords(
1121 result.exposure.wcs,
1122 sourceList=result.stars_footprints,
1123 include_covariance=self.config.do_include_astrometric_errors,
1124 )
1126 summary_stat_catalog = result.stars_footprints
1127 result.stars = result.stars_footprints.asAstropy()
1128 self.metadata["star_count"] = np.sum(~result.stars["sky_source"])
1130 # Validate the astrometric fit. Send in the stars_footprints
1131 # catalog so that its coords get set to NaN if the fit is deemed
1132 # a failure.
1133 self.astrometry.check(result.exposure, result.stars_footprints, len(astrometry_matches))
1134 result.stars = result.stars_footprints.asAstropy()
1135 have_fit_astrometry = True
1137 result.stars_footprints, photometry_matches, \
1138 photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints)
1139 have_fit_photometry = True
1140 self.metadata["photometry_matches_count"] = len(photometry_matches)
1141 # fit_photometry returns a new catalog, so we need a new astropy
1142 # table view.
1143 result.stars = result.stars_footprints.asAstropy()
1144 # summary stats don't make use of the calibrated fluxes, but we
1145 # might as well use the best catalog we've got in case that
1146 # changes, and help the old one get garbage-collected.
1147 summary_stat_catalog = result.stars_footprints
1148 if "photometry_matches" in self.config.optional_outputs:
1149 result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,
1150 photometry_meta)
1151 if "mask" in self.config.optional_outputs:
1152 result.mask = result.exposure.mask.clone()
1153 except pipeBase.AlgorithmError:
1154 if not have_fit_psf:
1155 result.exposure.setPsf(None)
1156 if not have_fit_astrometry:
1157 result.exposure.setWcs(None)
1158 if not have_fit_photometry:
1159 result.exposure.setPhotoCalib(None)
1160 # Summary stat calculations can handle missing components
1161 # gracefully, but we want to run them as late as possible (but
1162 # still before we calibrate pixels, if we do that at all).
1163 # So we run them after we succeed or if we get an AlgorithmError.
1164 # We intentionally don't use 'finally' here because we don't
1165 # want to run them if we get some other kind of error.
1166 self._summarize(result.exposure, summary_stat_catalog, result.background,
1167 summary=result.exposure.info.getSummaryStats())
1168 raise
1169 else:
1170 self._summarize(result.exposure, summary_stat_catalog, result.background,
1171 summary=result.exposure.info.getSummaryStats())
1173 # Output post-subtraction background stats to task metadata:
1174 # Specify the pixels flags to ignore, starting with those ignored
1175 # by the subtraction.
1176 bgIgnoredPixelMasks = lsst.meas.algorithms.SubtractBackgroundConfig().ignoredPixelMask.list()
1177 for maskName in self.config.background_stats_ignored_pixel_masks:
1178 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
1179 and maskName not in bgIgnoredPixelMasks):
1180 bgIgnoredPixelMasks += [maskName]
1182 statsCtrl = afwMath.StatisticsControl()
1183 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(bgIgnoredPixelMasks))
1185 statObj = afwMath.makeStatistics(result.exposure.getMaskedImage(), afwMath.MEDIAN, statsCtrl)
1186 median_bg, _ = statObj.getResult(afwMath.MEDIAN)
1187 self.metadata["bg_subtracted_skyPixel_instFlux_median"] = median_bg
1189 statObj = afwMath.makeStatistics(result.exposure.getMaskedImage(), afwMath.STDEV, statsCtrl)
1190 stdev_bg, _ = statObj.getResult(afwMath.STDEV)
1191 self.metadata["bg_subtracted_skyPixel_instFlux_stdev"] = stdev_bg
1193 self.metadata["bg_subtracted_skySource_flux_median"] = (
1194 np.median(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column])
1195 )
1196 self.metadata["bg_subtracted_skySource_flux_stdev"] = (
1197 np.std(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column])
1198 )
1200 if self.config.do_calibrate_pixels:
1201 self._apply_photometry(
1202 result.exposure,
1203 result.background,
1204 background_to_photometric_ratio=result.background_to_photometric_ratio,
1205 )
1206 result.applied_photo_calib = photo_calib
1207 else:
1208 result.applied_photo_calib = None
1210 self._recordMaskedPixelFractions(result.exposure)
1212 if self.config.run_sattle:
1213 # send boresight and timing information to sattle so the cache
1214 # is populated by the time we reach ip_diffim detectAndMeasure.
1215 try:
1216 populate_sattle_visit_cache(result.exposure.getInfo().getVisitInfo(),
1217 historical=self.config.sattle_historical)
1218 self.log.info('Successfully triggered load of sattle visit cache')
1219 except requests.exceptions.HTTPError:
1220 self.log.exception("Sattle visit cache update failed; continuing with image processing")
1222 return result
1224 def _apply_illumination_correction(self, exposure, background_flat, illumination_correction):
1225 """Apply the illumination correction to a background-flattened image.
1227 Parameters
1228 ----------
1229 exposure : `lsst.afw.image.Exposure`
1230 Exposure to convert to a photometric-flattened image.
1231 background_flat : `lsst.afw.image.Exposure`
1232 Flat image that had previously been applied to exposure.
1233 illumination_correction : `lsst.afw.image.Exposure`
1234 Illumination correction image to convert to photometric-flattened
1235 image.
1237 Returns
1238 -------
1239 background_to_photometric_ratio : `lsst.afw.image.Image`
1240 Ratio image to convert a photometric-flattened image to/from
1241 a background-flattened image. Will be None if task not
1242 configured to use the illumination correction.
1243 """
1244 if not self.config.do_illumination_correction:
1245 return None
1247 # From a raw image to a background-flattened image, we have:
1248 # bfi = image / background_flat
1249 # From a raw image to a photometric-flattened image, we have:
1250 # pfi = image / reference_flux_flat
1251 # pfi = image / (dome_flat * illumination_correction),
1252 # where the illumination correction contains the jacobian
1253 # of the wcs, converting to fluence units.
1254 # Currently background_flat == dome_flat, so we have for the
1255 # "background_to_photometric_ratio", the ratio of the background-
1256 # flattened image to the photometric-flattened image:
1257 # bfi / pfi = illumination_correction.
1259 background_to_photometric_ratio = illumination_correction.image.clone()
1261 # Dividing the ratio will convert a background-flattened image to
1262 # a photometric-flattened image.
1263 exposure.maskedImage /= background_to_photometric_ratio
1265 exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = True
1267 return background_to_photometric_ratio
1269 def _compute_psf(self, result, id_generator):
1270 """Find bright sources detected on an exposure and fit a PSF model to
1271 them, repairing likely cosmic rays before detection.
1273 Repair, detect, measure, and compute PSF twice, to ensure the PSF
1274 model does not include contributions from cosmic rays.
1276 Parameters
1277 ----------
1278 result : `lsst.pipe.base.Struct`
1279 Result struct that is modified to allow saving of partial outputs
1280 for some failure conditions. Should contain at least the following
1281 attributes:
1283 - exposure : `lsst.afw.image.Exposure`
1284 Exposure to detect and measure bright stars on.
1285 - background : `lsst.afw.math.BackgroundList` | `None`
1286 Background that was fit to the exposure during detection.
1287 - background_to_photometric_ratio : `lsst.afw.image.Image` | `None`
1288 Image to convert photometric-flattened image to
1289 background-flattened image.
1290 id_generator : `lsst.meas.base.IdGenerator`
1291 Object that generates source IDs and provides random seeds.
1293 Returns
1294 -------
1295 sources : `lsst.afw.table.SourceCatalog`
1296 Catalog of detected bright sources.
1297 cell_set : `lsst.afw.math.SpatialCellSet`
1298 PSF candidates returned by the psf determiner.
1299 adaptive_det_res_struct : `lsst.pipe.base.Struct`
1300 Result struct from the adaptive threshold detection.
1302 Notes
1303 -----
1304 This method modifies the exposure, background and
1305 background_to_photometric_ratio attributes of the result struct
1306 in-place.
1307 """
1308 def log_psf(msg, addToMetadata=False):
1309 """Log the parameters of the psf and background, with a prepended
1310 message. There is also the option to add the PSF sigma to the task
1311 metadata.
1313 Parameters
1314 ----------
1315 msg : `str`
1316 Message to prepend the log info with.
1317 addToMetadata : `bool`, optional
1318 Whether to add the final psf sigma value to the task
1319 metadata (the default is False).
1320 """
1321 position = result.exposure.psf.getAveragePosition()
1322 sigma = result.exposure.psf.computeShape(position).getDeterminantRadius()
1323 dimensions = result.exposure.psf.computeImage(position).getDimensions()
1324 if result.background is not None:
1325 median_background = np.median(result.background.getImage().array)
1326 else:
1327 median_background = 0.0
1328 self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",
1329 msg, sigma, dimensions, median_background)
1330 if addToMetadata:
1331 self.metadata["final_psf_sigma"] = sigma
1333 self.log.info("First pass detection with Gaussian PSF FWHM=%s pixels",
1334 self.config.install_simple_psf.fwhm)
1335 self.install_simple_psf.run(exposure=result.exposure)
1337 result.background = self.psf_subtract_background.run(
1338 exposure=result.exposure,
1339 backgroundToPhotometricRatio=result.background_to_photometric_ratio,
1340 ).background
1341 log_psf("Initial PSF:")
1342 self.psf_repair.run(exposure=result.exposure, keepCRs=True)
1344 table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())
1345 if not self.config.do_adaptive_threshold_detection:
1346 adaptive_det_res_struct = None
1347 # Re-estimate the background during this detection step, so that
1348 # measurement uses the most accurate background-subtraction.
1349 detections = self.psf_detection.run(
1350 table=table,
1351 exposure=result.exposure,
1352 background=result.background,
1353 backgroundToPhotometricRatio=result.background_to_photometric_ratio,
1354 )
1355 else:
1356 initialThreshold = self.config.psf_detection.thresholdValue
1357 initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier
1358 adaptive_det_res_struct = self.psf_adaptive_threshold_detection.run(
1359 table, result.exposure,
1360 initialThreshold=initialThreshold,
1361 initialThresholdMultiplier=initialThresholdMultiplier,
1362 )
1363 detections = adaptive_det_res_struct.detections
1365 self.metadata["initial_psf_positive_footprint_count"] = detections.numPos
1366 self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg
1367 self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks
1368 self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks
1370 if self.config.do_compute_shapelet_decomposition:
1371 # Compute a rough shapelet decomposition on the PSF stars to assess
1372 # focus (IQ).
1373 self.log.info("Computing rough shapelet decomposition for PSF star candidates.")
1374 seed = id_generator.catalog_id & 0xFFFFFFFF
1375 shapelets_results = self.compute_shapelet_decomposition.run(
1376 masked_image=result.exposure.getMaskedImage(), catalog=detections.sources, seed=seed
1377 )
1378 result.exposure.info.setSummaryStats(shapelets_results.shapelets_summary)
1379 detections.sources = shapelets_results.catalog
1381 self.psf_source_measurement.run(detections.sources, result.exposure)
1383 if self.config.do_compute_shapelet_decomposition:
1384 # If the source measurement succeeds, add a shapelets IQ score that
1385 # includes power from the median centroid offset between those used
1386 # in the shapelet decomposition and the measured and populated in
1387 # the centroid slot values (note that for very non-Gaussian
1388 # sources, this will often have reverted to the peak position).
1389 self.log.info("Computing shapelets_iq_score that includes power from the "
1390 "median centroid offset.")
1391 centroid_offset_scale = 0.5*np.sqrt(
1392 self.config.compute_shapelet_decomposition.max_footprint_area)
1393 shapelets_results = self.compute_shapelet_decomposition.compute_shapelets_iq_score(
1394 shapelets_results=shapelets_results,
1395 catalog=detections.sources,
1396 centroid_offset_scale=centroid_offset_scale,
1397 shapelets_base_name=self.compute_shapelet_decomposition.shapelets_base_name,
1398 log=self.log,
1399 )
1400 result.exposure.info.setSummaryStats(shapelets_results.shapelets_summary)
1402 psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)
1404 # This 2nd round of PSF fitting has been deemed unnecessary (and
1405 # sometimes even causing harm), so it is being skipped for the
1406 # adaptive threshold logic, but is left here for now to maintain
1407 # consistency with the previous logic for any users wanting to
1408 # revert to it.
1409 if not self.config.do_adaptive_threshold_detection:
1410 # Replace the initial PSF with something simpler for the second
1411 # repair/detect/measure/measure_psf step: this can help it
1412 # converge.
1413 self.install_simple_psf.run(exposure=result.exposure)
1415 log_psf("Rerunning with simple PSF:")
1416 # TODO investigation: Should we only re-run repair here, to use the
1417 # new PSF? Maybe we *do* need to re-run measurement with PsfFlux,
1418 # to use the fitted PSF?
1419 # TODO investigation: do we need a separate measurement task here
1420 # for the post-psf_measure_psf step, since we only want to do
1421 # PsfFlux and GaussianFlux *after* we have a PSF? Maybe that's not
1422 # relevant once DM-39203 is merged?
1423 self.psf_repair.run(exposure=result.exposure, keepCRs=True)
1424 # Re-estimate the background during this detection step, so that
1425 # measurement uses the most accurate background-subtraction.
1426 detections = self.psf_detection.run(
1427 table=table,
1428 exposure=result.exposure,
1429 background=result.background,
1430 backgroundToPhotometricRatio=result.background_to_photometric_ratio,
1431 )
1432 self.psf_source_measurement.run(detections.sources, result.exposure)
1433 psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)
1434 self.metadata["simple_psf_positive_footprint_count"] = detections.numPos
1435 self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg
1436 self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks
1437 self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks
1438 log_psf("Final PSF:", addToMetadata=True)
1440 # Final repair with final PSF, removing cosmic rays this time.
1441 self.psf_repair.run(exposure=result.exposure)
1442 # Final measurement with the CRs removed.
1443 self.psf_source_measurement.run(detections.sources, result.exposure)
1445 # PSF is set on exposure; candidates are returned to use for
1446 # calibration flux normalization and aperture corrections.
1447 return pipeBase.Struct(
1448 detections_sources=detections.sources,
1449 psf_result_cellSet=psf_result.cellSet,
1450 adaptive_det_res_struct=adaptive_det_res_struct,
1451 )
1453 def _measure_aperture_correction(self, exposure, bright_sources):
1454 """Measure and set the ApCorrMap on the Exposure, using
1455 previously-measured bright sources.
1457 This function first normalizes the calibration flux and then
1458 the full set of aperture corrections are measured relative
1459 to this normalized calibration flux.
1461 Parameters
1462 ----------
1463 exposure : `lsst.afw.image.Exposure`
1464 Exposure to set the ApCorrMap on.
1465 bright_sources : `lsst.afw.table.SourceCatalog`
1466 Catalog of detected bright sources; modified to include columns
1467 necessary for point source determination for the aperture
1468 correction calculation.
1469 """
1470 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(
1471 exposure=exposure,
1472 catalog=bright_sources,
1473 ).ap_corr_map
1475 ap_corr_map = self.measure_aperture_correction.run(exposure, bright_sources).apCorrMap
1477 # Need to merge the aperture correction map from the normalization.
1478 for key in norm_ap_corr_map:
1479 ap_corr_map[key] = norm_ap_corr_map[key]
1481 exposure.info.setApCorrMap(ap_corr_map)
1483 def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None,
1484 adaptive_det_res_struct=None, num_astrometry_matches=None):
1485 """Detect stars on an exposure that has a PSF model, and measure their
1486 PSF, circular aperture, compensated gaussian fluxes.
1488 Parameters
1489 ----------
1490 exposure : `lsst.afw.image.Exposure`
1491 Exposure to detect and measure stars on.
1492 background : `lsst.afw.math.BackgroundList`
1493 Background that was fit to the exposure during detection;
1494 modified in-place during subsequent detection.
1495 id_generator : `lsst.meas.base.IdGenerator`
1496 Object that generates source IDs and provides random seeds.
1497 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1498 Image to convert photometric-flattened image to
1499 background-flattened image.
1501 Returns
1502 -------
1503 stars : `SourceCatalog`
1504 Sources that are very likely to be stars, with a limited set of
1505 measurements performed on them.
1506 """
1507 table = afwTable.SourceTable.make(self.initial_stars_schema.schema,
1508 id_generator.make_table_id_factory())
1510 maxAdaptiveDetIter = 8
1511 threshFactor = 0.2
1512 if adaptive_det_res_struct is not None:
1513 for adaptiveDetIter in range(maxAdaptiveDetIter):
1514 adaptiveDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
1515 if self.config.do_remeasure_star_background:
1516 adaptiveDetectionConfig.reEstimateBackground = False
1517 else:
1518 adaptiveDetectionConfig.reEstimateBackground = True
1519 adaptiveDetectionConfig.includeThresholdMultiplier = 2.0
1520 psfThreshold = (
1521 adaptive_det_res_struct.thresholdValue*adaptive_det_res_struct.includeThresholdMultiplier
1522 )
1523 adaptiveDetectionConfig.thresholdValue = max(
1524 self.config.star_detection.thresholdValue,
1525 threshFactor*psfThreshold/adaptiveDetectionConfig.includeThresholdMultiplier
1526 )
1527 self.log.info("Using adaptive threshold detection (nIter = %d) with "
1528 "thresholdValue = %.2f and multiplier = %.1f",
1529 adaptiveDetIter, adaptiveDetectionConfig.thresholdValue,
1530 adaptiveDetectionConfig.includeThresholdMultiplier)
1531 adaptiveDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
1532 config=adaptiveDetectionConfig
1533 )
1534 detections = adaptiveDetectionTask.run(
1535 table=table,
1536 exposure=exposure,
1537 background=background,
1538 backgroundToPhotometricRatio=background_to_photometric_ratio,
1539 )
1540 nFootprint = len(detections.sources)
1541 nPeak = 0
1542 nIsolated = 0
1543 for src in detections.sources:
1544 nPeakSrc = len(src.getFootprint().getPeaks())
1545 if nPeakSrc == 1:
1546 nIsolated += 1
1547 nPeak += nPeakSrc
1548 minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches))
1549 if nFootprint > 0:
1550 self.log.info("Adaptive threshold detection _find_stars nIter: %d; "
1551 "nPeak/nFootprint = %.2f (max is 800); nIsolated = %d (min is %.1f).",
1552 adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated)
1553 if nPeak/nFootprint > 800 or nIsolated < minIsolated:
1554 threshFactor = max(0.01, 1.5*threshFactor)
1555 self.log.warning("nPeak/nFootprint = %.2f (max is 800); nIsolated = %d "
1556 "(min is %.1f).", nPeak/nFootprint, nIsolated, minIsolated)
1557 else:
1558 break
1559 else:
1560 threshFactor *= 0.75
1561 self.log.warning("No footprints detected on image. Decreasing threshold "
1562 "factor to %.2f. and rerunning.", threshFactor)
1563 else:
1564 # Re-estimate the background during this detection step, so that
1565 # measurement uses the most accurate background-subtraction.
1566 detections = self.star_detection.run(
1567 table=table,
1568 exposure=exposure,
1569 background=background,
1570 backgroundToPhotometricRatio=background_to_photometric_ratio,
1571 )
1572 sources = detections.sources
1573 self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources)
1575 n_sky_sources = np.sum(sources["sky_source"])
1576 if (self.config.do_downsample_footprints
1577 and (len(sources) - n_sky_sources) > self.config.downsample_max_footprints):
1578 if exposure.info.id is None:
1579 self.log.warning("Exposure does not have a proper id; using 0 seed for downsample.")
1580 seed = 0
1581 else:
1582 seed = exposure.info.id & 0xFFFFFFFF
1584 gen = np.random.RandomState(seed)
1586 # Isolate the sky sources before downsampling.
1587 indices = np.arange(len(sources))[~sources["sky_source"]]
1588 indices = gen.choice(
1589 indices,
1590 size=self.config.downsample_max_footprints,
1591 replace=False,
1592 )
1593 skyIndices, = np.where(sources["sky_source"])
1594 indices = np.concatenate((indices, skyIndices))
1596 self.log.info("Downsampling from %d to %d non-sky-source footprints.", len(sources), len(indices))
1598 sel = np.zeros(len(sources), dtype=bool)
1599 sel[indices] = True
1600 sources = sources[sel]
1602 # TODO investigation: Could this deblender throw away blends of
1603 # non-PSF sources?
1604 self.star_deblend.run(exposure=exposure, sources=sources)
1605 # The deblender may not produce a contiguous catalog; ensure
1606 # contiguity for subsequent tasks.
1607 if not sources.isContiguous():
1608 sources = sources.copy(deep=True)
1610 # Measure everything, and use those results to select only stars.
1611 self.star_measurement.run(sources, exposure)
1612 self.metadata["post_deblend_source_count"] = np.sum(~sources["sky_source"])
1613 self.metadata["failed_deblend_source_count"] = np.sum(
1614 ~sources["sky_source"] & sources["deblend_failed"]
1615 )
1616 self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"])
1617 self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"])
1619 # Run the normalization calibration flux task to apply the
1620 # normalization correction to create normalized
1621 # calibration fluxes.
1622 self.star_normalized_calibration_flux.run(exposure=exposure, catalog=sources)
1623 self.star_apply_aperture_correction.run(sources, exposure.apCorrMap)
1624 self.star_catalog_calculation.run(sources)
1625 self.star_set_primary_flags.run(sources)
1627 result = self.star_selector.run(sources)
1628 # The star selector may not produce a contiguous catalog.
1629 if not result.sourceCat.isContiguous():
1630 return result.sourceCat.copy(deep=True)
1631 else:
1632 return result.sourceCat
1634 def _match_psf_stars(self, psf_stars, stars, psfSigma=None):
1635 """Match calibration stars to psf stars, to identify which were psf
1636 candidates, and which were used or reserved during psf measurement
1637 and the astrometric fit.
1639 Parameters
1640 ----------
1641 psf_stars : `lsst.afw.table.SourceCatalog`
1642 PSF candidate stars that were sent to the psf determiner and
1643 used in the astrometric fit. Used to populate psf and astrometry
1644 related flag fields.
1645 stars : `lsst.afw.table.SourceCatalog`
1646 Stars that will be used for calibration; psf-related fields will
1647 be updated in-place.
1649 Notes
1650 -----
1651 This code was adapted from CalibrateTask.copyIcSourceFields().
1652 """
1653 control = afwTable.MatchControl()
1654 matchRadius = 3.0 if psfSigma is None else max(3.0, psfSigma) # in pixels
1655 # Return all matched objects, to separate blends.
1656 control.findOnlyClosest = False
1657 matches = afwTable.matchXy(psf_stars, stars, matchRadius, control)
1658 deblend_key = stars.schema["deblend_nChild"].asKey()
1659 matches = [m for m in matches if m[1].get(deblend_key) == 0]
1661 # Because we had to allow multiple matches to handle parents, we now
1662 # need to prune to the best (closest) matches.
1663 # Closest matches is a dict of psf_stars source ID to Match record
1664 # (psf_stars source, sourceCat source, distance in pixels).
1665 best = {}
1666 for match_psf, match_stars, d in matches:
1667 match = best.get(match_psf.getId())
1668 if match is None or d <= match[2]:
1669 best[match_psf.getId()] = (match_psf, match_stars, d)
1670 matches = list(best.values())
1671 # We'll use this to construct index arrays into each catalog.
1672 ids = np.array([(match_psf.getId(), match_stars.getId()) for match_psf, match_stars, d in matches]).T
1674 if (n_matches := len(matches)) == 0:
1675 raise NoPsfStarsToStarsMatchError(n_psf_stars=len(psf_stars), n_stars=len(stars))
1677 self.log.info("%d psf/astrometry stars out of %d matched %d calib stars",
1678 n_matches, len(psf_stars), len(stars))
1679 self.metadata["matched_psf_star_count"] = n_matches
1681 # Check that no stars sources are listed twice; we already know
1682 # that each match has a unique psf_stars id, due to using as the key
1683 # in best above.
1684 n_unique = len(set(m[1].getId() for m in matches))
1685 if n_unique != n_matches:
1686 self.log.warning("%d psf_stars matched only %d stars", n_matches, n_unique)
1688 # The indices of the IDs, so we can update the flag fields as arrays.
1689 idx_psf_stars = np.searchsorted(psf_stars["id"], ids[0])
1690 idx_stars = np.searchsorted(stars["id"], ids[1])
1691 for field in self.psf_fields:
1692 result = np.zeros(len(stars), dtype=bool)
1693 result[idx_stars] = psf_stars[field][idx_psf_stars]
1694 stars[field] = result
1695 stars['psf_id'][idx_stars] = psf_stars['id'][idx_psf_stars]
1696 stars['psf_max_value'][idx_stars] = psf_stars['psf_max_value'][idx_psf_stars]
1698 def _fit_astrometry(self, exposure, stars, exposure_region=None):
1699 """Fit an astrometric model to the data and return the reference
1700 matches used in the fit, and the fitted WCS.
1702 Parameters
1703 ----------
1704 exposure : `lsst.afw.image.Exposure`
1705 Exposure that is being fit, to get PSF and other metadata from.
1706 Modified to add the fitted skyWcs.
1707 stars : `SourceCatalog`
1708 Good stars selected for use in calibration, with RA/Dec coordinates
1709 computed from the pixel positions and fitted WCS.
1710 exposure_region : `lsst.sphgeom.Region`, optional
1711 The exposure region to use for the reference catalog filtering.
1712 If `None`, this region will be set as a padded bbox + current WCS
1713 of the exposure.
1715 Returns
1716 -------
1717 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1718 Reference/stars matches used in the fit.
1719 """
1720 initialWcs = exposure.wcs
1721 result = self.astrometry.run(stars, exposure, exposure_region=exposure_region)
1723 # Record astrometry stats to metadata
1724 self.metadata["astrometry_matches_count"] = len(result.matches)
1725 if exposure.wcs is not None:
1726 if initialWcs is not None:
1727 center = exposure.getBBox().getCenter()
1728 self.metadata['initial_to_final_wcs'] = (
1729 initialWcs.pixelToSky(center).separation(
1730 exposure.wcs.pixelToSky(center)
1731 ).asArcseconds()
1732 )
1733 else:
1734 self.metadata['initial_to_final_wcs'] = float("nan")
1735 self.metadata['astrom_offset_mean'] = exposure.metadata['SFM_ASTROM_OFFSET_MEAN']
1736 self.metadata['astrom_offset_std'] = exposure.metadata['SFM_ASTROM_OFFSET_STD']
1737 self.metadata['astrom_offset_median'] = exposure.metadata['SFM_ASTROM_OFFSET_MEDIAN']
1738 else:
1739 self.metadata['initial_to_final_wcs'] = float("nan")
1740 self.metadata['astrom_offset_mean'] = float("nan")
1741 self.metadata['astrom_offset_std'] = float("nan")
1742 self.metadata['astrom_offset_median'] = float("nan")
1744 return result.matches, result.matchMeta
1746 def _fit_photometry(self, exposure, stars):
1747 """Fit a photometric model to the data and return the reference
1748 matches used in the fit, and the fitted PhotoCalib.
1750 Parameters
1751 ----------
1752 exposure : `lsst.afw.image.Exposure`
1753 Exposure that is being fit, to get PSF and other metadata from.
1754 Has the fit `lsst.afw.image.PhotoCalib` attached, with pixel values
1755 unchanged.
1756 stars : `lsst.afw.table.SourceCatalog`
1757 Good stars selected for use in calibration.
1758 background : `lsst.afw.math.BackgroundList`
1759 Background that was fit to the exposure during detection of the
1760 above stars.
1762 Returns
1763 -------
1764 calibrated_stars : `lsst.afw.table.SourceCatalog`
1765 Star catalog with flux/magnitude columns computed from the fitted
1766 photoCalib (instFlux columns are retained as well).
1767 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1768 Reference/stars matches used in the fit.
1769 matchMeta : `lsst.daf.base.PropertyList`
1770 Metadata needed to unpersist matches, as returned by the matcher.
1771 photo_calib : `lsst.afw.image.PhotoCalib`
1772 Photometric calibration that was fit to the star catalog.
1773 """
1774 result = self.photometry.run(exposure, stars)
1775 calibrated_stars = result.photoCalib.calibrateCatalog(stars)
1776 exposure.setPhotoCalib(result.photoCalib)
1777 return calibrated_stars, result.matches, result.matchMeta, result.photoCalib
1779 def _apply_photometry(self, exposure, background, background_to_photometric_ratio=None):
1780 """Apply the photometric model attached to the exposure to the
1781 exposure's pixels and an associated background model.
1783 Parameters
1784 ----------
1785 exposure : `lsst.afw.image.Exposure`
1786 Exposure with the target `lsst.afw.image.PhotoCalib` attached.
1787 On return, pixel values will be calibrated and an identity
1788 photometric transform will be attached.
1789 background : `lsst.afw.math.BackgroundList`
1790 Background model to convert to nanojansky units in place.
1791 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1792 Image to convert photometric-flattened image to
1793 background-flattened image.
1794 """
1795 photo_calib = exposure.getPhotoCalib()
1796 exposure.maskedImage = photo_calib.calibrateImage(exposure.maskedImage)
1797 identity = afwImage.PhotoCalib(1.0,
1798 photo_calib.getCalibrationErr(),
1799 bbox=exposure.getBBox())
1800 exposure.setPhotoCalib(identity)
1801 exposure.metadata["BUNIT"] = "nJy"
1803 assert photo_calib._isConstant, \
1804 "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that."
1806 for bg in background:
1807 # The statsImage is a view, but we can't assign to a function call
1808 # in python.
1809 binned_image = bg[0].getStatsImage()
1810 binned_image *= photo_calib.getCalibrationMean()
1812 def _summarize(self, exposure, stars, background, summary=None):
1813 """Compute summary statistics on the exposure and update in-place the
1814 calibrations attached to it.
1816 Parameters
1817 ----------
1818 exposure : `lsst.afw.image.Exposure`
1819 Exposure that was calibrated, to get PSF and other metadata from.
1820 Should be in instrumental units with the photometric calibration
1821 attached.
1822 Modified to contain the computed summary statistics.
1823 stars : `SourceCatalog`
1824 Good stars selected used in calibration.
1825 background : `lsst.afw.math.BackgroundList`
1826 Background that was fit to the exposure during detection of the
1827 above stars. Should be in instrumental units.
1828 """
1829 summary = self.compute_summary_stats.run(exposure, stars, background, summary=summary)
1830 exposure.info.setSummaryStats(summary)
1832 def _recordMaskedPixelFractions(self, exposure):
1833 """Record the fraction of all the pixels in an exposure
1834 that are masked with a given flag. Each fraction is
1835 recorded in the task metadata. One record per flag type.
1837 Parameters
1838 ----------
1839 exposure : `lsst.afw.image.ExposureF`
1840 The target exposure to calculate masked pixel fractions for.
1841 """
1843 mask = exposure.mask
1844 maskPlanes = list(mask.getMaskPlaneDict().keys())
1845 for maskPlane in maskPlanes:
1846 self.metadata[f"{maskPlane.lower()}_mask_fraction"] = (
1847 evaluateMaskFraction(mask, maskPlane)
1848 )
1850 def _update_wcs_with_camera_model(self, exposure, cameraModel, exposure_record=None):
1851 """Replace the existing WCS with one generated using the pointing
1852 and the input camera model.
1854 If the current detector does not exist in the cameraModel, the pointing
1855 will still get updated, but the original distortion model will be used.
1857 Parameters
1858 ----------
1859 exposure : `lsst.afw.image.ExposureF`
1860 The exposure whose WCS will be updated.
1861 cameraModel : `lsst.afw.cameraGeom.Camera`
1862 Camera to be used when constructing updated WCS.
1863 exposure_record : `lsst.daf.butler.DimensionRecord`, optional
1864 Butler metadata for the ``exposure`` dimension. Used if
1865 constructing an updated WCS instead of the boresight and
1866 orientation angle from the visit info.
1867 """
1868 if cameraModel.get(exposure.detector.getId()):
1869 self.log.info("Updating WCS with the provided camera model.")
1870 detector = cameraModel[exposure.detector.getId()]
1871 else:
1872 self.log.warning(
1873 "useButlerCamera=True, but detector %s is not available in the provided camera. The "
1874 "astrometry fit will use the initial distortion model for this detector.",
1875 exposure.detector.getId())
1876 detector = exposure.detector
1877 if exposure_record is None:
1878 boresight = exposure.visitInfo.getBoresightRaDec()
1879 orientation = exposure.visitInfo.getBoresightRotAngle()
1880 else:
1881 boresight = SpherePoint(exposure_record.tracking_ra, exposure_record.tracking_dec, degrees)
1882 orientation = exposure_record.sky_angle * degrees
1883 self.log.info(
1884 "Pointing from exposure record is %.6f deg (%.3f arcsec) from initial pointing; "
1885 "orientation difference is %.6f deg (%.3f arcsec).",
1886 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asDegrees(),
1887 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asArcseconds(),
1888 (orientation - exposure.visitInfo.getBoresightRotAngle()).asDegrees(),
1889 (orientation - exposure.visitInfo.getBoresightRotAngle()).asArcseconds(),
1890 )
1891 newWcs = createInitialSkyWcsFromBoresight(boresight, orientation, detector)
1892 exposure.setWcs(newWcs)
1894 def _update_wcs_to_exposure_record(self, exposure, exposure_record):
1895 """Replace the existing WCS with the one from the butler derived
1896 exposure record pointing.
1898 Parameters
1899 ----------
1900 exposure : `lsst.afw.image.ExposureF`
1901 The exposure whose WCS will be updated.
1902 exposure_record : `lsst.daf.butler.DimensionRecord`, optional
1903 Butler metadata for the ``exposure`` dimension. Used if
1904 constructing an updated WCS instead of the boresight and
1905 orientation angle from the visit info.
1906 """
1907 detector = exposure.detector
1908 boresight = SpherePoint(exposure_record.tracking_ra, exposure_record.tracking_dec, degrees)
1909 orientation = exposure_record.sky_angle * degrees
1910 self.log.info(
1911 "Pointing from exposure record is %.6f deg (%.3f arcsec) from initial pointing; "
1912 "orientation difference is %.6f deg (%.3f arcsec).",
1913 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asDegrees(),
1914 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asArcseconds(),
1915 (orientation - exposure.visitInfo.getBoresightRotAngle()).asDegrees(),
1916 (orientation - exposure.visitInfo.getBoresightRotAngle()).asArcseconds(),
1917 )
1918 newWcs = createInitialSkyWcsFromBoresight(boresight, orientation, detector)
1919 exposure.setWcs(newWcs)
1921 def _compute_mask_fraction(self, mask, mask_planes, bad_mask_planes):
1922 """Evaluate the fraction of masked pixels in a (set of) mask plane(s).
1924 Parameters
1925 ----------
1926 mask : `lsst.afw.image.Mask`
1927 The mask on which to evaluate the fraction.
1928 mask_planes : `list`, `str`
1929 The mask planes for which to evaluate the fraction.
1930 bad_mask_planes : `list`, `str`
1931 The mask planes to exclude from the computation.
1933 Returns
1934 -------
1935 detected_fraction : `float`
1936 The calculated fraction of masked pixels
1937 """
1938 bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes)
1939 n_good_pix = np.sum(mask.array & bad_pixel_mask == 0)
1940 if n_good_pix == 0:
1941 detected_fraction = float("nan")
1942 return detected_fraction
1943 detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes)
1944 n_detected_pix = np.sum((mask.array & detected_pixel_mask != 0)
1945 & (mask.array & bad_pixel_mask == 0))
1946 detected_fraction = n_detected_pix/n_good_pix
1947 return detected_fraction
1949 def _compute_per_amp_fraction(self, exposure, detected_fraction, mask_planes, bad_mask_planes):
1950 """Evaluate the maximum per-amplifier fraction of masked pixels.
1952 Parameters
1953 ----------
1954 exposure : `lsst.afw.image.ExposureF`
1955 The exposure on which to compute the per-amp masked fraction.
1956 detected_fraction : `float`
1957 The current detected_fraction of the ``mask_planes`` for the
1958 full detector.
1959 mask_planes : `list`, `str`
1960 The mask planes for which to evaluate the fraction.
1961 bad_mask_planes : `list`, `str`
1962 The mask planes to exclude from the computation.
1964 Returns
1965 -------
1966 n_above_max_per_amp : `int`
1967 The number of amplifiers with masked fractions above a maximum
1968 value (set by the current full-detector ``detected_fraction``).
1969 highest_detected_fraction_per_amp : `float`
1970 The highest value of the per-amplifier fraction of masked pixels.
1971 no_zero_det_amps : `bool`
1972 A boolean representing whether any of the amplifiers has zero
1973 masked pixels.
1974 """
1975 highest_detected_fraction_per_amp = -9.99
1976 n_above_max_per_amp = 0
1977 n_no_zero_det_amps = 0
1978 no_zero_det_amps = True
1979 amps = exposure.detector.getAmplifiers()
1980 if amps is not None:
1981 for ia, amp in enumerate(amps):
1982 amp_bbox = amp.getBBox()
1983 exp_bbox = exposure.getBBox()
1984 if not exp_bbox.contains(amp_bbox):
1985 self.log.info("Bounding box of amplifier (%s) does not fit in exposure's "
1986 "bounding box (%s). Skipping...", amp_bbox, exp_bbox)
1987 continue
1988 sub_image = exposure.subset(amp.getBBox())
1989 detected_fraction_amp = self._compute_mask_fraction(sub_image.mask,
1990 mask_planes,
1991 bad_mask_planes)
1992 self.log.debug("Current detected fraction for amplifier %s = %.3f",
1993 amp.getName(), detected_fraction_amp)
1994 if detected_fraction_amp < 0.002:
1995 n_no_zero_det_amps += 1
1996 if n_no_zero_det_amps > 2:
1997 no_zero_det_amps = False
1998 break
1999 highest_detected_fraction_per_amp = max(detected_fraction_amp,
2000 highest_detected_fraction_per_amp)
2001 if highest_detected_fraction_per_amp > min(0.998, max(0.8, 3.0*detected_fraction)):
2002 n_above_max_per_amp += 1
2003 if n_above_max_per_amp > 2:
2004 break
2005 else:
2006 self.log.info("No amplifier object for detector %d, so skipping per-amp "
2007 "detection fraction checks.", exposure.detector.getId())
2008 return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps
2010 def _remeasure_star_background(self, result, background_to_photometric_ratio=None):
2011 """Remeasure the exposure's background with iterative adaptive
2012 threshold detection.
2014 Parameters
2015 ----------
2016 result : `lsst.pipe.base.Struct`
2017 The current state of the result Strut from the run method of
2018 CalibrateImageTask. Will be modified in place.
2019 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
2020 Image to convert photometric-flattened image to
2021 background-flattened image.
2023 Returns
2024 -------
2025 result : `lsst.pipe.base.Struct`
2026 The modified result Struct with the new background subtracted.
2027 """
2028 # Restore the previously measured background and remeasure it
2029 # using an adaptive threshold detection iteration to ensure a
2030 # "Goldilocks Zone" for the fraction of detected pixels.
2031 backgroundOrig = result.background.clone()
2032 median_background_orig = np.median(backgroundOrig.getImage().array)
2033 self.log.info("Original median_background = %.2f", median_background_orig)
2034 result.exposure.image.array += result.background.getImage().array
2035 result.background = afwMath.BackgroundList()
2037 origMask = result.exposure.mask.clone()
2038 bad_mask_planes = ["BAD", "EDGE", "NO_DATA"]
2039 detected_mask_planes = ["DETECTED", "DETECTED_NEGATIVE"]
2040 detected_fraction_orig = self._compute_mask_fraction(result.exposure.mask,
2041 detected_mask_planes,
2042 bad_mask_planes)
2043 minDetFracForFinalBg = 0.02
2044 maxDetFracForFinalBg = 0.93
2045 # Dilate the current detected mask planes and don't clear
2046 # them in the detection step.
2047 nPixToDilate = 10
2048 maxAdaptiveDetIter = 10
2049 for dilateIter in range(maxAdaptiveDetIter):
2050 dilatedMask = origMask.clone()
2051 for maskName in detected_mask_planes:
2052 # Compute the grown detection mask plane using SpanSet
2053 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)
2054 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)
2055 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)
2056 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())
2057 # Clear the detected mask plane
2058 detectedMask = dilatedMask.getMaskPlane(maskName)
2059 dilatedMask.clearMaskPlane(detectedMask)
2060 # Set the mask plane to the dilated one
2061 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)
2063 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,
2064 detected_mask_planes,
2065 bad_mask_planes)
2066 if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1:
2067 break
2068 else:
2069 nPixToDilate -= 1
2071 result.exposure.mask = dilatedMask
2072 self.log.debug("detected_fraction_orig = %.3f; detected_fraction_dilated = %.3f",
2073 detected_fraction_orig, detected_fraction_dilated)
2074 n_above_max_per_amp = -99
2075 highest_detected_fraction_per_amp = float("nan")
2077 # Check the per-amplifier detection fractions.
2078 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \
2079 self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated,
2080 detected_mask_planes, bad_mask_planes)
2081 self.log.debug("Dilated mask: n_above_max_per_amp = %d; "
2082 "highest_detected_fraction_per_amp = %.3f",
2083 n_above_max_per_amp, highest_detected_fraction_per_amp)
2085 bgIgnoreMasksToAdd = ["SAT", "SUSPECT", "SPIKE"]
2086 detected_fraction = 1.0
2087 nFootprintTemp = 1e12
2088 starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()
2089 for maskName in bgIgnoreMasksToAdd:
2090 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
2091 and maskName not in starBackgroundDetectionConfig.background.ignoredPixelMask):
2092 starBackgroundDetectionConfig.background.ignoredPixelMask += [maskName]
2093 starBackgroundDetectionConfig.doTempLocalBackground = False
2094 starBackgroundDetectionConfig.nSigmaToGrow = 70.0
2095 starBackgroundDetectionConfig.reEstimateBackground = False
2096 starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0
2097 starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background_orig)
2098 starBackgroundDetectionConfig.thresholdType = "pixel_stdev"
2100 n_above_max_per_amp = -99
2101 highest_detected_fraction_per_amp = float("nan")
2102 doCheckPerAmpDetFraction = True
2104 minFootprints = self.config.star_background_min_footprints
2105 maxIter = 40
2106 for nIter in range(maxIter):
2107 currentThresh = starBackgroundDetectionConfig.thresholdValue
2108 nZeroEncountered = 0
2109 if nFootprintTemp == 0:
2110 zeroFactor = min(0.98, 0.9 + 0.01*nZeroEncountered)
2111 starBackgroundDetectionConfig.thresholdValue = zeroFactor*currentThresh
2112 self.log.warning("No footprints detected. Decreasing threshold to %.2f and rerunning.",
2113 starBackgroundDetectionConfig.thresholdValue)
2114 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
2115 config=starBackgroundDetectionConfig)
2116 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)
2117 tempDetections = starBackgroundDetectionTask.run(
2118 table=table, exposure=result.exposure, clearMask=True)
2119 nFootprintTemp = len(tempDetections.sources)
2120 minFootprints = max(self.config.star_background_min_footprints,
2121 int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))
2122 minFootprints = min(200, minFootprints)
2123 nZeroEncountered += 1
2124 if nFootprintTemp >= minFootprints:
2125 detected_fraction = self._compute_mask_fraction(result.exposure.mask,
2126 detected_mask_planes,
2127 bad_mask_planes)
2128 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
2129 "DETECTED_NEGATIVE in star_background_detection = %.3f "
2130 "(max is %.3f; min is %.3f) nFootprint = %d (current min is %d)",
2131 nIter, starBackgroundDetectionConfig.thresholdValue,
2132 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
2133 nFootprintTemp, minFootprints)
2134 self.metadata['adaptive_threshold_value'] = starBackgroundDetectionConfig.thresholdValue
2136 break
2137 else:
2138 # Still not enough footprints, so make sure this loop is
2139 # entered again.
2140 if nFootprintTemp > 0 and nFootprintTemp < minFootprints:
2141 nFootprintTemp = 0
2142 continue
2143 if detected_fraction > maxDetFracForFinalBg or nFootprintTemp <= minFootprints:
2144 starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh
2145 if nFootprintTemp < minFootprints and detected_fraction > 0.9*maxDetFracForFinalBg:
2146 if nFootprintTemp == 1:
2147 starBackgroundDetectionConfig.thresholdValue = 1.4*currentThresh
2148 else:
2149 starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh
2151 if n_above_max_per_amp > 1:
2152 starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh
2153 if detected_fraction < minDetFracForFinalBg:
2154 starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh
2155 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(
2156 config=starBackgroundDetectionConfig)
2157 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)
2158 tempDetections = starBackgroundDetectionTask.run(
2159 table=table, exposure=result.exposure, clearMask=True)
2160 result.exposure.mask |= dilatedMask
2161 nFootprintTemp = len(tempDetections.sources)
2162 minFootprints = max(self.config.star_background_min_footprints,
2163 int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))
2164 minFootprints = min(200, minFootprints)
2165 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
2166 bad_mask_planes)
2167 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "
2168 "DETECTED_NEGATIVE in star_background_detection = %.3f "
2169 "(max is %.3f; min is %.3f); nFootprint = %d (current min is %d)",
2170 nIter, starBackgroundDetectionConfig.thresholdValue,
2171 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,
2172 nFootprintTemp, minFootprints)
2173 self.metadata['adaptive_threshold_value'] = starBackgroundDetectionConfig.thresholdValue
2175 n_amp = len(result.exposure.detector.getAmplifiers())
2176 if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg:
2177 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \
2178 self._compute_per_amp_fraction(result.exposure, detected_fraction,
2179 detected_mask_planes, bad_mask_planes)
2181 if not no_zero_det_amps:
2182 starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh
2184 if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg
2185 and n_above_max_per_amp < int(0.75*n_amp)
2186 and no_zero_det_amps
2187 and nFootprintTemp >= minFootprints):
2188 if (n_above_max_per_amp < max(1, int(0.15*n_amp))
2189 or detected_fraction < 0.85*maxDetFracForFinalBg):
2190 break
2191 else:
2192 starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh
2193 self.log.debug("Number of amplifiers with detected fraction above the maximum "
2194 "(%.2f) excedes the maximum allowed (%d >= %d). Making a small "
2195 "tweak to the threshold (from %.2f to %.2f) and rerunning.",
2196 maxDetFracForFinalBg, n_above_max_per_amp, int(0.75*n_amp),
2197 currentThresh, starBackgroundDetectionConfig.thresholdValue)
2198 self.log.debug("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp))
2200 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,
2201 bad_mask_planes)
2202 self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f "
2203 "(highest per amp section = %.5f)",
2204 detected_fraction, highest_detected_fraction_per_amp)
2206 if detected_fraction > maxDetFracForFinalBg:
2207 result.exposure.mask = dilatedMask
2208 self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE "
2209 "was too large in star_background_detection = %.3f (max = %.3f). "
2210 "Reverting to dilated mask from PSF detection.",
2211 detected_fraction, maxDetFracForFinalBg)
2212 self.metadata['adaptive_threshold_value'] = float("nan")
2213 star_background = self.star_background.run(
2214 exposure=result.exposure,
2215 backgroundToPhotometricRatio=background_to_photometric_ratio,
2216 ).background
2217 result.background.append(star_background[0])
2219 median_background = np.median(result.background.getImage().array)
2220 self.log.info("New initial median_background = %.2f", median_background)
2222 # Perform one more round of background subtraction that is just an
2223 # overall pedestal (order = 0). This is intended to account for
2224 # any potential gross oversubtraction imposed by the higher-order
2225 # subtraction.
2226 # Dilate DETECTED mask a bit more if it's below 50% detected.
2227 nPixToDilate = 2
2228 if detected_fraction < 0.5:
2229 dilatedMask = result.exposure.mask.clone()
2230 for maskName in detected_mask_planes:
2231 # Compute the grown detection mask plane using SpanSet
2232 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)
2233 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)
2234 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)
2235 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())
2236 # Clear the detected mask plane
2237 detectedMask = dilatedMask.getMaskPlane(maskName)
2238 dilatedMask.clearMaskPlane(detectedMask)
2239 # Set the mask plane to the dilated one
2240 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)
2242 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,
2243 detected_mask_planes,
2244 bad_mask_planes)
2245 result.exposure.mask = dilatedMask
2246 self.log.debug("detected_fraction_orig = %.3f; detected_fraction_dilated = %.3f",
2247 detected_fraction_orig, detected_fraction_dilated)
2249 bbox = result.exposure.getBBox()
2251 # Now measure a final background pedestal level (largely to account
2252 # for the over-subtraction often seen in the first pass, especially
2253 # in high fill-factor fields). Since the pedestal measurement can be
2254 # sensitive to the bin size, start an iteration with a small bin size,
2255 # then double it on each iteration until the relative or absolute
2256 # change criterion is met. If those are never achieved, the iteration
2257 # stops when the bin size gets bigger than the exposure's bounding box.
2259 # Initialize the parameters for the pedestal iteration.
2260 pedestalBackgroundConfig = lsst.meas.algorithms.SubtractBackgroundConfig()
2261 for maskName in bgIgnoreMasksToAdd:
2262 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()
2263 and maskName not in pedestalBackgroundConfig.ignoredPixelMask):
2264 pedestalBackgroundConfig.ignoredPixelMask += [maskName]
2265 pedestalBackgroundConfig.statisticsProperty = "MEDIAN"
2266 pedestalBackgroundConfig.approxOrderX = 0
2267 pedestalBackgroundConfig.binSize = min(32, int(math.ceil(max(bbox.width, bbox.height)/4.0)))
2268 self.log.info("Initial pedestal binSize = %d pixels", pedestalBackgroundConfig.binSize)
2269 cumulativePedestalLevel = 0.0
2270 # When the cumulative pedestal value changes by less than 5% from one
2271 # bin size to the next we assume we are converged.
2272 relativeStoppingCriterion = 0.05
2273 # When the cumulative pedestal value changes by less than 0.5 counts
2274 # (electrons or adu) from one bin size to the next, we assume we are
2275 # converged.
2276 absoluteStoppingCriterion = 0.5
2277 inPedestalIteration = True
2278 while inPedestalIteration:
2279 cumulativePedestalPrev = cumulativePedestalLevel
2280 pedestalBackgroundTask = lsst.meas.algorithms.SubtractBackgroundTask(
2281 config=pedestalBackgroundConfig)
2282 pedestalBackgroundList = pedestalBackgroundTask.run(
2283 exposure=result.exposure,
2284 background=result.background,
2285 backgroundToPhotometricRatio=background_to_photometric_ratio,
2286 ).background
2287 # Isolate the most recent pedestal background measured to log the
2288 # current computed value and add it to the cumulative value.
2289 pedestalBackground = afwMath.BackgroundList()
2290 pedestalBackground.append(pedestalBackgroundList[-1])
2291 pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0]
2292 cumulativePedestalLevel += pedestalBackgroundLevel
2293 absoluteDelta = abs(cumulativePedestalLevel - cumulativePedestalPrev)
2294 if cumulativePedestalPrev != 0.0:
2295 relativeDelta = abs(absoluteDelta/cumulativePedestalPrev)
2296 else:
2297 relativeDelta = 1.0
2298 self.log.info("Subtracted pedestalBackgroundLevel = %.4f (cumulative = %.4f; "
2299 "relativeDelta = %.4f; absoluteDelta = %.3f)",
2300 pedestalBackgroundLevel, cumulativePedestalLevel,
2301 relativeDelta, absoluteDelta)
2303 if (relativeDelta < relativeStoppingCriterion
2304 or absoluteDelta < absoluteStoppingCriterion):
2305 # We are converged.
2306 inPedestalIteration = False
2307 else:
2308 # We have not yet converged; grow the bin size if possible.
2309 if pedestalBackgroundConfig.binSize*2 < 2*max(bbox.width, bbox.height):
2310 pedestalBackgroundConfig.binSize = int(pedestalBackgroundConfig.binSize*2)
2311 self.log.info("Increasing pedestal binSize to %d pixels and remeasuring.",
2312 pedestalBackgroundConfig.binSize)
2313 else:
2314 # We have reached the maximum bin size.
2315 self.log.warning("Reached maximum bin size. Exiting pedestal loop without meeting "
2316 "the convergence criteria (relativeDelta <= %.2f or "
2317 "absoluteDelta <= %.2f).",
2318 relativeStoppingCriterion, absoluteStoppingCriterion)
2319 inPedestalIteration = False
2320 self.log.info("Final subtracted cumulativePedestalBackgroundLevel = %.4f", cumulativePedestalLevel)
2322 # Clear detected mask plane before final round of detection
2323 mask = result.exposure.mask
2324 for mp in detected_mask_planes:
2325 if mp not in mask.getMaskPlaneDict():
2326 mask.addMaskPlane(mp)
2327 mask &= ~mask.getPlaneBitMask(detected_mask_planes)
2329 return result