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