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