Coverage for python/lsst/pipe/tasks/calibrateImage.py: 18%
479 statements
« prev ^ index » next coverage.py v7.11.0, created at 2025-11-06 02:11 -0800
1# This file is part of pipe_tasks.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22__all__ = ["CalibrateImageTask", "CalibrateImageConfig", "NoPsfStarsToStarsMatchError",
23 "AllCentroidsFlaggedError"]
25import numpy as np
26import requests
27import os
29import lsst.afw.table as afwTable
30import lsst.afw.image as afwImage
31from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache
32import lsst.meas.algorithms
33import lsst.meas.algorithms.installGaussianPsf
34import lsst.meas.algorithms.measureApCorr
35import lsst.meas.algorithms.setPrimaryFlags
36import lsst.meas.base
37import lsst.meas.astrom
38import lsst.meas.deblender
39import lsst.meas.extensions.shapeHSM
40from lsst.obs.base import createInitialSkyWcs
41import lsst.pex.config as pexConfig
42import lsst.pipe.base as pipeBase
43from lsst.pipe.base import connectionTypes
44from lsst.utils.timer import timeMethod
46from . import measurePsf, repair, photoCal, computeExposureSummaryStats, snapCombine, diffractionSpikeMask
49class AllCentroidsFlaggedError(pipeBase.AlgorithmError):
50 """Raised when there are no valid centroids during psf fitting.
51 """
52 def __init__(self, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size):
53 msg = (f"All source centroids (out of {n_sources}) flagged during PSF fitting. "
54 "Original image PSF is likely unuseable; best-fit PSF shape parameters: "
55 f"Ixx={psf_shape_ixx}, Iyy={psf_shape_iyy}, Ixy={psf_shape_ixy}, size={psf_size}"
56 )
57 super().__init__(msg)
58 self.n_sources = n_sources
59 self.psf_shape_ixx = psf_shape_ixx
60 self.psf_shape_iyy = psf_shape_iyy
61 self.psf_shape_ixy = psf_shape_ixy
62 self.psf_size = psf_size
64 @property
65 def metadata(self):
66 return {"n_sources": self.n_sources,
67 "psf_shape_ixx": self.psf_shape_ixx,
68 "psf_shape_iyy": self.psf_shape_iyy,
69 "psf_shape_ixy": self.psf_shape_ixy,
70 "psf_size": self.psf_size
71 }
74class NoPsfStarsToStarsMatchError(pipeBase.AlgorithmError):
75 """Raised when there are no matches between the psf_stars and stars
76 catalogs.
77 """
78 def __init__(self, *, n_psf_stars, n_stars):
79 msg = (f"No psf stars out of {n_psf_stars} matched {n_stars} calib stars."
80 " Downstream processes probably won't have useful {calib_type} stars in this case."
81 " Is `star_selector` too strict or is this a bad image?")
82 super().__init__(msg)
83 self.n_psf_stars = n_psf_stars
84 self.n_stars = n_stars
86 @property
87 def metadata(self):
88 return {"n_psf_stars": self.n_psf_stars,
89 "n_stars": self.n_stars
90 }
93class CalibrateImageConnections(pipeBase.PipelineTaskConnections,
94 dimensions=("instrument", "visit", "detector")):
96 astrometry_ref_cat = connectionTypes.PrerequisiteInput(
97 doc="Reference catalog to use for astrometric calibration.",
98 name="gaia_dr3_20230707",
99 storageClass="SimpleCatalog",
100 dimensions=("skypix",),
101 deferLoad=True,
102 multiple=True,
103 )
104 photometry_ref_cat = connectionTypes.PrerequisiteInput(
105 doc="Reference catalog to use for photometric calibration.",
106 name="ps1_pv3_3pi_20170110",
107 storageClass="SimpleCatalog",
108 dimensions=("skypix",),
109 deferLoad=True,
110 multiple=True
111 )
112 camera_model = connectionTypes.PrerequisiteInput(
113 doc="Camera distortion model to use for astrometric calibration.",
114 name="astrometry_camera",
115 storageClass="Camera",
116 dimensions=("instrument", "physical_filter"),
117 isCalibration=True,
118 minimum=0, # Can fall back on WCS already attached to exposure.
119 )
120 exposures = connectionTypes.Input(
121 doc="Exposure (or two snaps) to be calibrated, and detected and measured on.",
122 name="postISRCCD",
123 storageClass="Exposure",
124 multiple=True, # to handle 1 exposure or 2 snaps
125 dimensions=["instrument", "exposure", "detector"],
126 )
128 background_flat = connectionTypes.PrerequisiteInput(
129 name="flat",
130 doc="Flat calibration frame used for background correction.",
131 storageClass="ExposureF",
132 dimensions=["instrument", "detector", "physical_filter"],
133 isCalibration=True,
134 )
135 illumination_correction = connectionTypes.PrerequisiteInput(
136 name="illuminationCorrection",
137 doc="Illumination correction frame.",
138 storageClass="Exposure",
139 dimensions=["instrument", "detector", "physical_filter"],
140 isCalibration=True,
141 )
143 # outputs
144 initial_stars_schema = connectionTypes.InitOutput(
145 doc="Schema of the output initial stars catalog.",
146 name="initial_stars_schema",
147 storageClass="SourceCatalog",
148 )
150 # TODO DM-38732: We want some kind of flag on Exposures/Catalogs to make
151 # it obvious which components had failed to be computed/persisted.
152 exposure = connectionTypes.Output(
153 doc="Photometrically calibrated, background-subtracted exposure with fitted calibrations and "
154 "summary statistics. To recover the original exposure, first add the background "
155 "(`initial_pvi_background`), and then uncalibrate (divide by `initial_photoCalib_detector`).",
156 name="initial_pvi",
157 storageClass="ExposureF",
158 dimensions=("instrument", "visit", "detector"),
159 )
160 stars = connectionTypes.Output(
161 doc="Catalog of unresolved sources detected on the calibrated exposure.",
162 name="initial_stars_detector",
163 storageClass="ArrowAstropy",
164 dimensions=["instrument", "visit", "detector"],
165 )
166 stars_footprints = connectionTypes.Output(
167 doc="Catalog of unresolved sources detected on the calibrated exposure; "
168 "includes source footprints.",
169 name="initial_stars_footprints_detector",
170 storageClass="SourceCatalog",
171 dimensions=["instrument", "visit", "detector"],
172 )
173 applied_photo_calib = connectionTypes.Output(
174 doc=(
175 "Photometric calibration that was applied to exposure's pixels. "
176 "This connection is disabled when do_calibrate_pixels=False."
177 ),
178 name="initial_photoCalib_detector",
179 storageClass="PhotoCalib",
180 dimensions=("instrument", "visit", "detector"),
181 )
182 background = connectionTypes.Output(
183 doc="Background models estimated during calibration task; calibrated to be in nJy units.",
184 name="initial_pvi_background",
185 storageClass="Background",
186 dimensions=("instrument", "visit", "detector"),
187 )
188 background_to_photometric_ratio = connectionTypes.Output(
189 doc="Ratio of a background-flattened image to a photometric-flattened image. Only persisted "
190 "if do_illumination_correction is True.",
191 name="background_to_photometric_ratio",
192 storageClass="Image",
193 dimensions=("instrument", "visit", "detector"),
194 )
196 # Optional outputs
197 mask = connectionTypes.Output(
198 doc="The mask plane of the calibrated exposure, written as a separate "
199 "output so that it can be passed to downstream tasks even if the "
200 "exposure is removed to save storage.",
201 name="preliminary_visit_mask",
202 storageClass="Mask",
203 dimensions=("instrument", "visit", "detector"),
204 )
205 psf_stars_footprints = connectionTypes.Output(
206 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination; "
207 "includes source footprints.",
208 name="initial_psf_stars_footprints_detector",
209 storageClass="SourceCatalog",
210 dimensions=["instrument", "visit", "detector"],
211 )
212 psf_stars = connectionTypes.Output(
213 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination.",
214 name="initial_psf_stars_detector",
215 storageClass="ArrowAstropy",
216 dimensions=["instrument", "visit", "detector"],
217 )
218 astrometry_matches = connectionTypes.Output(
219 doc="Source to reference catalog matches from the astrometry solver.",
220 name="initial_astrometry_match_detector",
221 storageClass="Catalog",
222 dimensions=("instrument", "visit", "detector"),
223 )
224 photometry_matches = connectionTypes.Output(
225 doc="Source to reference catalog matches from the photometry solver.",
226 name="initial_photometry_match_detector",
227 storageClass="Catalog",
228 dimensions=("instrument", "visit", "detector"),
229 )
231 def __init__(self, *, config=None):
232 super().__init__(config=config)
233 if "psf_stars" not in config.optional_outputs:
234 del self.psf_stars
235 if "psf_stars_footprints" not in config.optional_outputs:
236 del self.psf_stars_footprints
237 if "astrometry_matches" not in config.optional_outputs:
238 del self.astrometry_matches
239 if "photometry_matches" not in config.optional_outputs:
240 del self.photometry_matches
241 if "mask" not in config.optional_outputs:
242 del self.mask
243 if not config.do_calibrate_pixels:
244 del self.applied_photo_calib
245 if not config.do_illumination_correction:
246 del self.background_flat
247 del self.illumination_correction
248 del self.background_to_photometric_ratio
249 if not config.useButlerCamera:
250 del self.camera_model
253class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CalibrateImageConnections):
254 optional_outputs = pexConfig.ListField(
255 doc="Which optional outputs to save (as their connection name)?",
256 dtype=str,
257 # TODO: note somewhere to disable this for benchmarking, but should
258 # we always have it on for production runs?
259 default=["psf_stars", "psf_stars_footprints", "astrometry_matches", "photometry_matches", "mask"],
260 optional=False
261 )
263 # To generate catalog ids consistently across subtasks.
264 id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()
266 snap_combine = pexConfig.ConfigurableField(
267 target=snapCombine.SnapCombineTask,
268 doc="Task to combine two snaps to make one exposure.",
269 )
271 # subtasks used during psf characterization
272 install_simple_psf = pexConfig.ConfigurableField(
273 target=lsst.meas.algorithms.installGaussianPsf.InstallGaussianPsfTask,
274 doc="Task to install a simple PSF model into the input exposure to use "
275 "when detecting bright sources for PSF estimation.",
276 )
277 psf_repair = pexConfig.ConfigurableField(
278 target=repair.RepairTask,
279 doc="Task to repair cosmic rays on the exposure before PSF determination.",
280 )
281 psf_subtract_background = pexConfig.ConfigurableField(
282 target=lsst.meas.algorithms.SubtractBackgroundTask,
283 doc="Task to perform intial background subtraction, before first detection pass.",
284 )
285 psf_detection = pexConfig.ConfigurableField(
286 target=lsst.meas.algorithms.SourceDetectionTask,
287 doc="Task to detect sources for PSF determination."
288 )
289 psf_source_measurement = pexConfig.ConfigurableField(
290 target=lsst.meas.base.SingleFrameMeasurementTask,
291 doc="Task to measure sources to be used for psf estimation."
292 )
293 psf_measure_psf = pexConfig.ConfigurableField(
294 target=measurePsf.MeasurePsfTask,
295 doc="Task to measure the psf on bright sources."
296 )
297 psf_normalized_calibration_flux = pexConfig.ConfigurableField(
298 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
299 doc="Task to normalize the calibration flux (e.g. compensated tophats) "
300 "for the bright stars used for psf estimation.",
301 )
303 # TODO DM-39203: we can remove aperture correction from this task once we are
304 # using the shape-based star/galaxy code.
305 measure_aperture_correction = pexConfig.ConfigurableField(
306 target=lsst.meas.algorithms.measureApCorr.MeasureApCorrTask,
307 doc="Task to compute the aperture correction from the bright stars."
308 )
310 # subtasks used during star measurement
311 star_detection = pexConfig.ConfigurableField(
312 target=lsst.meas.algorithms.SourceDetectionTask,
313 doc="Task to detect stars to return in the output catalog."
314 )
315 star_sky_sources = pexConfig.ConfigurableField(
316 target=lsst.meas.algorithms.SkyObjectsTask,
317 doc="Task to generate sky sources ('empty' regions where there are no detections).",
318 )
319 do_downsample_footprints = pexConfig.Field(
320 dtype=bool,
321 default=False,
322 doc="Downsample footprints prior to deblending to optimize speed?",
323 )
324 downsample_max_footprints = pexConfig.Field(
325 dtype=int,
326 default=1000,
327 doc="Maximum number of non-sky-source footprints to use if do_downsample_footprints is True,",
328 )
329 star_deblend = pexConfig.ConfigurableField(
330 target=lsst.meas.deblender.SourceDeblendTask,
331 doc="Split blended sources into their components."
332 )
333 star_measurement = pexConfig.ConfigurableField(
334 target=lsst.meas.base.SingleFrameMeasurementTask,
335 doc="Task to measure stars to return in the output catalog."
336 )
337 star_normalized_calibration_flux = pexConfig.ConfigurableField(
338 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,
339 doc="Task to apply the normalization for calibration fluxes (e.g. compensated tophats) "
340 "for the final output star catalog.",
341 )
342 star_apply_aperture_correction = pexConfig.ConfigurableField(
343 target=lsst.meas.base.ApplyApCorrTask,
344 doc="Task to apply aperture corrections to the selected stars."
345 )
346 star_catalog_calculation = pexConfig.ConfigurableField(
347 target=lsst.meas.base.CatalogCalculationTask,
348 doc="Task to compute extendedness values on the star catalog, "
349 "for the star selector to remove extended sources."
350 )
351 star_set_primary_flags = pexConfig.ConfigurableField(
352 target=lsst.meas.algorithms.setPrimaryFlags.SetPrimaryFlagsTask,
353 doc="Task to add isPrimary to the catalog."
354 )
355 star_selector = lsst.meas.algorithms.sourceSelectorRegistry.makeField(
356 default="science",
357 doc="Task to select reliable stars to use for calibration."
358 )
360 # final calibrations and statistics
361 astrometry = pexConfig.ConfigurableField(
362 target=lsst.meas.astrom.AstrometryTask,
363 doc="Task to perform astrometric calibration to fit a WCS.",
364 )
365 astrometry_ref_loader = pexConfig.ConfigField(
366 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
367 doc="Configuration of reference object loader for astrometric fit.",
368 )
369 photometry = pexConfig.ConfigurableField(
370 target=photoCal.PhotoCalTask,
371 doc="Task to perform photometric calibration to fit a PhotoCalib.",
372 )
373 photometry_ref_loader = pexConfig.ConfigField(
374 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,
375 doc="Configuration of reference object loader for photometric fit.",
376 )
377 doMaskDiffractionSpikes = pexConfig.Field(
378 dtype=bool,
379 default=False,
380 doc="If True, load the photometric reference catalog again but select "
381 "only bright stars. Use the bright star catalog to set the SPIKE "
382 "mask for regions likely contaminated by diffraction spikes.",
383 )
384 diffractionSpikeMask = pexConfig.ConfigurableField(
385 target=diffractionSpikeMask.DiffractionSpikeMaskTask,
386 doc="Task to identify and mask the diffraction spikes of bright stars.",
387 )
389 compute_summary_stats = pexConfig.ConfigurableField(
390 target=computeExposureSummaryStats.ComputeExposureSummaryStatsTask,
391 doc="Task to to compute summary statistics on the calibrated exposure."
392 )
394 do_illumination_correction = pexConfig.Field(
395 dtype=bool,
396 default=False,
397 doc="If True, apply the illumination correction. This assumes that the "
398 "input image has already been flat-fielded such that it is suitable "
399 "for background subtraction.",
400 )
402 do_calibrate_pixels = pexConfig.Field(
403 dtype=bool,
404 default=True,
405 doc=(
406 "If True, apply the photometric calibration to the image pixels "
407 "and background model, and attach an identity PhotoCalib to "
408 "the output image to reflect this. If False`, leave the image "
409 "and background uncalibrated and attach the PhotoCalib that maps "
410 "them to physical units."
411 )
412 )
414 do_include_astrometric_errors = pexConfig.Field(
415 dtype=bool,
416 default=True,
417 doc="If True, include astrometric errors in the output catalog.",
418 )
420 run_sattle = pexConfig.Field(
421 dtype=bool,
422 default=False,
423 doc="If True, the sattle service will populate a cache for later use "
424 "in ip_diffim.detectAndMeasure alert verification."
425 )
427 sattle_historical = pexConfig.Field(
428 dtype=bool,
429 default=False,
430 doc="If re-running a pipeline that requires sattle, this should be set "
431 "to True. This will populate sattle's cache with the historic data "
432 "closest in time to the exposure.",
433 )
435 useButlerCamera = pexConfig.Field(
436 dtype=bool,
437 default=False,
438 doc="If True, use a camera distortion model generated elsewhere in "
439 "the pipeline combined with the telescope boresight as a starting "
440 "point for fitting the WCS, instead of using the WCS attached to "
441 "the exposure, which is generated from the boresight and the "
442 "camera model from the obs_* package."
443 )
445 def setDefaults(self):
446 super().setDefaults()
448 # Use a very broad PSF here, to throughly reject CRs.
449 # TODO investigation: a large initial psf guess may make stars look
450 # like CRs for very good seeing images.
451 self.install_simple_psf.fwhm = 4
453 # S/N>=50 sources for PSF determination, but detection to S/N=10.
454 # The thresholdValue sets the minimum flux in a pixel to be included in the
455 # footprint, while peaks are only detected when they are above
456 # thresholdValue * includeThresholdMultiplier. The low thresholdValue
457 # ensures that the footprints are large enough for the noise replacer
458 # to mask out faint undetected neighbors that are not to be measured.
459 self.psf_detection.thresholdValue = 10.0
460 self.psf_detection.includeThresholdMultiplier = 5.0
461 # TODO investigation: Probably want False here, but that may require
462 # tweaking the background spatial scale, to make it small enough to
463 # prevent extra peaks in the wings of bright objects.
464 self.psf_detection.doTempLocalBackground = False
465 # NOTE: we do want reEstimateBackground=True in psf_detection, so that
466 # each measurement step is done with the best background available.
468 # Minimal measurement plugins for PSF determination.
469 # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use
470 # shapeHSM/moments for star/galaxy separation.
471 # TODO DM-39203: we can remove aperture correction from this task once
472 # we are using the shape-based star/galaxy code.
473 self.psf_source_measurement.plugins = ["base_PixelFlags",
474 "base_SdssCentroid",
475 "ext_shapeHSM_HsmSourceMoments",
476 "base_CircularApertureFlux",
477 "base_GaussianFlux",
478 "base_PsfFlux",
479 "base_CompensatedTophatFlux",
480 ]
481 self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
482 # Only measure apertures we need for PSF measurement.
483 self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
484 self.psf_source_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
485 # TODO DM-40843: Remove this line once this is the psfex default.
486 self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \
487 "base_CircularApertureFlux_12_0_instFlux"
489 # No extendeness information available: we need the aperture
490 # corrections to determine that.
491 self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False
492 self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]
493 self.measure_aperture_correction.sourceSelector["science"].flags.bad = []
495 # Detection for good S/N for astrometry/photometry and other
496 # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks.
497 self.star_detection.thresholdValue = 5.0
498 self.star_detection.includeThresholdMultiplier = 2.0
499 self.star_measurement.plugins = ["base_PixelFlags",
500 "base_SdssCentroid",
501 "ext_shapeHSM_HsmSourceMoments",
502 "ext_shapeHSM_HsmPsfMoments",
503 "base_GaussianFlux",
504 "base_PsfFlux",
505 "base_CircularApertureFlux",
506 "base_ClassificationSizeExtendedness",
507 "base_CompensatedTophatFlux",
508 ]
509 self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
510 self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
511 # Only measure the apertures we need for star selection.
512 self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
513 self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
515 # We measure and apply the normalization aperture correction with the
516 # psf_normalized_calibration_flux task, and we only apply the normalization
517 # aperture correction for the full list of stars.
518 self.star_normalized_calibration_flux.do_measure_ap_corr = False
520 # Select stars with reliable measurements and no bad flags.
521 self.star_selector["science"].doFlags = True
522 self.star_selector["science"].doUnresolved = True
523 self.star_selector["science"].doSignalToNoise = True
524 self.star_selector["science"].signalToNoise.minimum = 10.0
525 # Keep sky sources in the output catalog, even though they aren't
526 # wanted for calibration.
527 self.star_selector["science"].doSkySources = True
528 # Set the flux and error fields
529 self.star_selector["science"].signalToNoise.fluxField = "slot_CalibFlux_instFlux"
530 self.star_selector["science"].signalToNoise.errField = "slot_CalibFlux_instFluxErr"
532 # Use the affine WCS fitter (assumes we have a good camera geometry).
533 self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)
534 # phot_g_mean is the primary Gaia band for all input bands.
535 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"
537 # Only reject sky sources; we already selected good stars.
538 self.astrometry.sourceSelector["science"].doFlags = True
539 self.astrometry.sourceSelector["science"].flags.good = ["calib_psf_candidate"]
540 self.astrometry.sourceSelector["science"].flags.bad = []
541 self.astrometry.sourceSelector["science"].doUnresolved = False
542 self.astrometry.sourceSelector["science"].doIsolated = False
543 self.astrometry.sourceSelector["science"].doRequirePrimary = False
544 self.photometry.match.sourceSelection.doFlags = True
545 self.photometry.match.sourceSelection.flags.bad = ["sky_source"]
546 # Unset the (otherwise reasonable, but we've already made the
547 # selections we want above) selection settings in PhotoCalTask.
548 self.photometry.match.sourceSelection.doRequirePrimary = False
549 self.photometry.match.sourceSelection.doUnresolved = False
551 def validate(self):
552 super().validate()
554 # Ensure that the normalization calibration flux tasks
555 # are configured correctly.
556 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:
557 msg = ("psf_normalized_calibration_flux task must be configured with do_measure_ap_corr=True "
558 "or else the normalization and calibration flux will not be properly measured.")
559 raise pexConfig.FieldValidationError(
560 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,
561 )
562 if self.star_normalized_calibration_flux.do_measure_ap_corr:
563 msg = ("star_normalized_calibration_flux task must be configured with do_measure_ap_corr=False "
564 "to apply the previously measured normalization to the full catalog of calibration "
565 "fluxes.")
566 raise pexConfig.FieldValidationError(
567 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,
568 )
570 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,
571 # because they'd be running too early to pick up the fitted PhotoCalib
572 # and WCS.
573 if "base_LocalWcs" in self.psf_source_measurement.plugins.names:
574 raise pexConfig.FieldValidationError(
575 CalibrateImageConfig.psf_source_measurement,
576 self,
577 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
578 )
579 if "base_LocalWcs" in self.star_measurement.plugins.names:
580 raise pexConfig.FieldValidationError(
581 CalibrateImageConfig.star_measurement,
582 self,
583 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
584 )
585 if "base_LocalPhotoCalib" in self.psf_source_measurement.plugins.names:
586 raise pexConfig.FieldValidationError(
587 CalibrateImageConfig.psf_source_measurement,
588 self,
589 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
590 "as it would be run before the photometry fit."
591 )
592 if "base_LocalPhotoCalib" in self.star_measurement.plugins.names:
593 raise pexConfig.FieldValidationError(
594 CalibrateImageConfig.star_measurement,
595 self,
596 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
597 "as it would be run before the photometry fit."
598 )
600 # Check for illumination correction and background consistency.
601 if self.do_illumination_correction:
602 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:
603 raise pexConfig.FieldValidationError(
604 CalibrateImageConfig.psf_subtract_background,
605 self,
606 "CalibrateImageTask.psf_subtract_background must be configured with "
607 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
608 )
609 if self.psf_detection.reEstimateBackground:
610 if not self.psf_detection.doApplyFlatBackgroundRatio:
611 raise pexConfig.FieldValidationError(
612 CalibrateImageConfig.psf_detection,
613 self,
614 "CalibrateImageTask.psf_detection background must be configured with "
615 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
616 )
617 if self.star_detection.reEstimateBackground:
618 if not self.star_detection.doApplyFlatBackgroundRatio:
619 raise pexConfig.FieldValidationError(
620 CalibrateImageConfig.star_detection,
621 self,
622 "CalibrateImageTask.star_detection background must be configured with "
623 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
624 )
626 if self.run_sattle:
627 if not os.getenv("SATTLE_URI_BASE"):
628 raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self,
629 "Sattle requested but URI environment variable not set.")
632class CalibrateImageTask(pipeBase.PipelineTask):
633 """Compute the PSF, aperture corrections, astrometric and photometric
634 calibrations, and summary statistics for a single science exposure, and
635 produce a catalog of brighter stars that were used to calibrate it.
637 Parameters
638 ----------
639 initial_stars_schema : `lsst.afw.table.Schema`
640 Schema of the initial_stars output catalog.
641 """
642 _DefaultName = "calibrateImage"
643 ConfigClass = CalibrateImageConfig
645 def __init__(self, initial_stars_schema=None, **kwargs):
646 super().__init__(**kwargs)
647 self.makeSubtask("snap_combine")
649 # PSF determination subtasks
650 self.makeSubtask("install_simple_psf")
651 self.makeSubtask("psf_repair")
652 self.makeSubtask("psf_subtract_background")
653 self.psf_schema = afwTable.SourceTable.makeMinimalSchema()
654 self.psf_schema.addField(
655 'psf_max_value',
656 type=np.float32,
657 doc="PSF max value.",
658 )
659 afwTable.CoordKey.addErrorFields(self.psf_schema)
660 self.makeSubtask("psf_detection", schema=self.psf_schema)
661 self.makeSubtask("psf_source_measurement", schema=self.psf_schema)
662 self.makeSubtask("psf_measure_psf", schema=self.psf_schema)
663 self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema)
665 self.makeSubtask("measure_aperture_correction", schema=self.psf_schema)
666 self.makeSubtask("astrometry", schema=self.psf_schema)
668 # star measurement subtasks
669 if initial_stars_schema is None:
670 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()
672 # These fields let us track which sources were used for psf modeling,
673 # astrometric fitting, and aperture correction calculations.
674 self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved",
675 "calib_astrometry_used",
676 # TODO DM-39203: these can be removed once apcorr is gone.
677 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",
678 "apcorr_base_PsfFlux_used",)
679 for field in self.psf_fields:
680 item = self.psf_schema.find(field)
681 initial_stars_schema.addField(item.getField())
682 id_type = self.psf_schema["id"].asField().getTypeString()
683 psf_max_value_type = self.psf_schema['psf_max_value'].asField().getTypeString()
684 initial_stars_schema.addField("psf_id",
685 type=id_type,
686 doc="id of this source in psf_stars; 0 if there is no match.")
687 initial_stars_schema.addField("psf_max_value",
688 type=psf_max_value_type,
689 doc="Maximum value in the star image used to train PSF.")
691 afwTable.CoordKey.addErrorFields(initial_stars_schema)
692 self.makeSubtask("star_detection", schema=initial_stars_schema)
693 self.makeSubtask("star_sky_sources", schema=initial_stars_schema)
694 self.makeSubtask("star_deblend", schema=initial_stars_schema)
695 self.makeSubtask("star_measurement", schema=initial_stars_schema)
696 self.makeSubtask("star_normalized_calibration_flux", schema=initial_stars_schema)
698 self.makeSubtask("star_apply_aperture_correction", schema=initial_stars_schema)
699 self.makeSubtask("star_catalog_calculation", schema=initial_stars_schema)
700 self.makeSubtask("star_set_primary_flags", schema=initial_stars_schema, isSingleFrame=True)
701 self.makeSubtask("star_selector")
702 self.makeSubtask("photometry", schema=initial_stars_schema)
703 if self.config.doMaskDiffractionSpikes:
704 self.makeSubtask("diffractionSpikeMask")
705 self.makeSubtask("compute_summary_stats")
707 # The final catalog will have calibrated flux columns, which we add to
708 # the init-output schema by calibrating our zero-length catalog with an
709 # arbitrary dummy PhotoCalib. We also use this schema to initialze
710 # the stars catalog in order to ensure it's the same even when we hit
711 # an error (and write partial outputs) before calibrating the catalog
712 # - note that calibrateCatalog will happily reuse existing output
713 # columns.
714 dummy_photo_calib = afwImage.PhotoCalib(1.0, 0, bbox=lsst.geom.Box2I())
715 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(
716 afwTable.SourceCatalog(initial_stars_schema)
717 )
719 def runQuantum(self, butlerQC, inputRefs, outputRefs):
720 inputs = butlerQC.get(inputRefs)
721 exposures = inputs.pop("exposures")
723 id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)
725 astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
726 dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],
727 refCats=inputs.pop("astrometry_ref_cat"),
728 name=self.config.connections.astrometry_ref_cat,
729 config=self.config.astrometry_ref_loader, log=self.log)
730 self.astrometry.setRefObjLoader(astrometry_loader)
732 photometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(
733 dataIds=[ref.datasetRef.dataId for ref in inputRefs.photometry_ref_cat],
734 refCats=inputs.pop("photometry_ref_cat"),
735 name=self.config.connections.photometry_ref_cat,
736 config=self.config.photometry_ref_loader, log=self.log)
737 self.photometry.match.setRefObjLoader(photometry_loader)
739 if self.config.doMaskDiffractionSpikes:
740 # Use the same photometry reference catalog for the bright star mask
741 self.diffractionSpikeMask.setRefObjLoader(photometry_loader)
743 if self.config.do_illumination_correction:
744 background_flat = inputs.pop("background_flat")
745 illumination_correction = inputs.pop("illumination_correction")
746 else:
747 background_flat = None
748 illumination_correction = None
750 camera_model = None
751 if self.config.useButlerCamera:
752 if "camera_model" in inputs:
753 camera_model = inputs.pop("camera_model")
754 else:
755 self.log.warning("useButlerCamera=True, but camera is not available for filter %s. The "
756 "astrometry fit will use the WCS already attached to the exposure.",
757 exposures[0].filter.bandLabel)
759 # This should not happen with a properly configured execution context.
760 assert not inputs, "runQuantum got more inputs than expected"
762 # Specify the fields that `annotate` needs below, to ensure they
763 # exist, even as None.
764 result = pipeBase.Struct(
765 exposure=None,
766 stars_footprints=None,
767 psf_stars_footprints=None,
768 background_to_photometric_ratio=None,
769 )
770 try:
771 self.run(
772 exposures=exposures,
773 result=result,
774 id_generator=id_generator,
775 background_flat=background_flat,
776 illumination_correction=illumination_correction,
777 camera_model=camera_model,
778 )
779 except pipeBase.AlgorithmError as e:
780 error = pipeBase.AnnotatedPartialOutputsError.annotate(
781 e,
782 self,
783 result.exposure,
784 result.psf_stars_footprints,
785 result.stars_footprints,
786 log=self.log
787 )
788 butlerQC.put(result, outputRefs)
789 raise error from e
791 butlerQC.put(result, outputRefs)
793 @timeMethod
794 def run(
795 self,
796 *,
797 exposures,
798 id_generator=None,
799 result=None,
800 background_flat=None,
801 illumination_correction=None,
802 camera_model=None,
803 ):
804 """Find stars and perform psf measurement, then do a deeper detection
805 and measurement and calibrate astrometry and photometry from that.
807 Parameters
808 ----------
809 exposures : `lsst.afw.image.Exposure` or `list` [`lsst.afw.image.Exposure`]
810 Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter.
811 Modified in-place during processing if only one is passed.
812 If two exposures are passed, treat them as snaps and combine
813 before doing further processing.
814 id_generator : `lsst.meas.base.IdGenerator`, optional
815 Object that generates source IDs and provides random seeds.
816 result : `lsst.pipe.base.Struct`, optional
817 Result struct that is modified to allow saving of partial outputs
818 for some failure conditions. If the task completes successfully,
819 this is also returned.
820 background_flat : `lsst.afw.image.Exposure`, optional
821 Background flat-field image.
822 illumination_correction : `lsst.afw.image.Exposure`, optional
823 Illumination correction image.
824 camera_model : `lsst.afw.cameraGeom.Camera`, optional
825 Camera to be used if constructing updated WCS.
827 Returns
828 -------
829 result : `lsst.pipe.base.Struct`
830 Results as a struct with attributes:
832 ``exposure``
833 Calibrated exposure, with pixels in nJy units.
834 (`lsst.afw.image.Exposure`)
835 ``stars``
836 Stars that were used to calibrate the exposure, with
837 calibrated fluxes and magnitudes.
838 (`astropy.table.Table`)
839 ``stars_footprints``
840 Footprints of stars that were used to calibrate the exposure.
841 (`lsst.afw.table.SourceCatalog`)
842 ``psf_stars``
843 Stars that were used to determine the image PSF.
844 (`astropy.table.Table`)
845 ``psf_stars_footprints``
846 Footprints of stars that were used to determine the image PSF.
847 (`lsst.afw.table.SourceCatalog`)
848 ``background``
849 Background that was fit to the exposure when detecting
850 ``stars``. (`lsst.afw.math.BackgroundList`)
851 ``applied_photo_calib``
852 Photometric calibration that was fit to the star catalog and
853 applied to the exposure. (`lsst.afw.image.PhotoCalib`)
854 This is `None` if ``config.do_calibrate_pixels`` is `False`.
855 ``astrometry_matches``
856 Reference catalog stars matches used in the astrometric fit.
857 (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`)
858 ``photometry_matches``
859 Reference catalog stars matches used in the photometric fit.
860 (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`)
861 ``mask``
862 Copy of the mask plane of `exposure`.
863 (`lsst.afw.image.Mask`)
864 """
865 if result is None:
866 result = pipeBase.Struct()
867 if id_generator is None:
868 id_generator = lsst.meas.base.IdGenerator()
870 result.exposure = self.snap_combine.run(exposures).exposure
871 self._recordMaskedPixelFractions(result.exposure)
872 self.log.info("Initial PhotoCalib: %s", result.exposure.getPhotoCalib())
874 result.exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = False
876 # Check input image processing.
877 if self.config.do_illumination_correction:
878 if not result.exposure.metadata.get("LSST ISR FLAT APPLIED", False):
879 raise pipeBase.InvalidQuantumError(
880 "Cannot use do_illumination_correction with an image that has not had a flat applied",
881 )
883 if camera_model:
884 if camera_model.get(result.exposure.detector.getId()):
885 self.log.info("Updating WCS with the provided camera model.")
886 self._update_wcs_with_camera_model(result.exposure, camera_model)
887 else:
888 self.log.warning(
889 "useButlerCamera=True, but detector %s is not available in the provided camera. The "
890 "astrometry fit will use the WCS already attached to the exposure.",
891 result.exposure.detector.getId())
893 result.background = None
894 summary_stat_catalog = None
895 # Some exposure components are set to initial placeholder objects
896 # while we try to bootstrap them. If we fail before we fit for them,
897 # we want to reset those components to None so the placeholders don't
898 # masquerade as the real thing.
899 have_fit_psf = False
900 have_fit_astrometry = False
901 have_fit_photometry = False
902 try:
903 result.background_to_photometric_ratio = self._apply_illumination_correction(
904 result.exposure,
905 background_flat,
906 illumination_correction,
907 )
909 result.psf_stars_footprints, result.background, _ = self._compute_psf(
910 result.exposure,
911 id_generator,
912 background_to_photometric_ratio=result.background_to_photometric_ratio,
913 )
914 have_fit_psf = True
916 # Check if all centroids have been flagged. This should happen
917 # after the PSF is fit and recorded because even a junky PSF
918 # model is informative.
919 if result.psf_stars_footprints["slot_Centroid_flag"].all():
920 psf_shape = result.exposure.psf.computeShape(result.exposure.psf.getAveragePosition())
921 raise AllCentroidsFlaggedError(
922 n_sources=len(result.psf_stars_footprints),
923 psf_shape_ixx=psf_shape.getIxx(),
924 psf_shape_iyy=psf_shape.getIyy(),
925 psf_shape_ixy=psf_shape.getIxy(),
926 psf_size=psf_shape.getDeterminantRadius(),
927 )
929 self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)
930 result.psf_stars = result.psf_stars_footprints.asAstropy()
931 # Run astrometry using PSF candidate stars
932 astrometry_matches, astrometry_meta = self._fit_astrometry(
933 result.exposure, result.psf_stars_footprints
934 )
935 self.metadata["astrometry_matches_count"] = len(astrometry_matches)
936 if "astrometry_matches" in self.config.optional_outputs:
937 result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches,
938 astrometry_meta)
939 result.psf_stars = result.psf_stars_footprints.asAstropy()
941 # Run the stars_detection subtask for the photometric calibration.
942 result.stars_footprints = self._find_stars(
943 result.exposure,
944 result.background,
945 id_generator,
946 background_to_photometric_ratio=result.background_to_photometric_ratio,
947 )
948 self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints)
950 # Update the source cooordinates with the current wcs.
951 afwTable.updateSourceCoords(
952 result.exposure.wcs,
953 sourceList=result.stars_footprints,
954 include_covariance=self.config.do_include_astrometric_errors
955 )
957 summary_stat_catalog = result.stars_footprints
958 result.stars = result.stars_footprints.asAstropy()
959 self.metadata["star_count"] = np.sum(~result.stars["sky_source"])
961 # Validate the astrometric fit. Send in the stars_footprints
962 # catalog so that its coords get set to NaN if the fit is deemed
963 # a failure.
964 self.astrometry.check(result.exposure, result.stars_footprints, len(astrometry_matches))
965 result.stars = result.stars_footprints.asAstropy()
966 have_fit_astrometry = True
968 result.stars_footprints, photometry_matches, \
969 photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints)
970 have_fit_photometry = True
971 self.metadata["photometry_matches_count"] = len(photometry_matches)
972 # fit_photometry returns a new catalog, so we need a new astropy table view.
973 result.stars = result.stars_footprints.asAstropy()
974 # summary stats don't make use of the calibrated fluxes, but we
975 # might as well use the best catalog we've got in case that
976 # changes, and help the old one get garbage-collected.
977 summary_stat_catalog = result.stars_footprints
978 if "photometry_matches" in self.config.optional_outputs:
979 result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,
980 photometry_meta)
981 if self.config.doMaskDiffractionSpikes:
982 self.diffractionSpikeMask.run(result.exposure)
983 if "mask" in self.config.optional_outputs:
984 result.mask = result.exposure.mask.clone()
985 except pipeBase.AlgorithmError:
986 if not have_fit_psf:
987 result.exposure.setPsf(None)
988 if not have_fit_astrometry:
989 result.exposure.setWcs(None)
990 if not have_fit_photometry:
991 result.exposure.setPhotoCalib(None)
992 # Summary stat calculations can handle missing components gracefully,
993 # but we want to run them as late as possible (but still before we
994 # calibrate pixels, if we do that at all).
995 # So we run them after we succeed or if we get an AlgorithmError. We
996 # intentionally don't use 'finally' here because we don't want to run
997 # them if we get some other kind of error.
998 self._summarize(result.exposure, summary_stat_catalog, result.background)
999 raise
1000 else:
1001 self._summarize(result.exposure, summary_stat_catalog, result.background)
1003 if self.config.do_calibrate_pixels:
1004 self._apply_photometry(
1005 result.exposure,
1006 result.background,
1007 background_to_photometric_ratio=result.background_to_photometric_ratio,
1008 )
1009 result.applied_photo_calib = photo_calib
1010 else:
1011 result.applied_photo_calib = None
1013 if self.config.run_sattle:
1014 # send boresight and timing information to sattle so the cache
1015 # is populated by the time we reach ip_diffim detectAndMeasure.
1016 try:
1017 populate_sattle_visit_cache(result.exposure.getInfo().getVisitInfo(),
1018 historical=self.config.sattle_historical)
1019 self.log.info('Successfully triggered load of sattle visit cache')
1020 except requests.exceptions.HTTPError:
1021 self.log.exception("Sattle visit cache update failed; continuing with image processing")
1023 return result
1025 def _apply_illumination_correction(self, exposure, background_flat, illumination_correction):
1026 """Apply the illumination correction to a background-flattened image.
1028 Parameters
1029 ----------
1030 exposure : `lsst.afw.image.Exposure`
1031 Exposure to convert to a photometric-flattened image.
1032 background_flat : `lsst.afw.image.Exposure`
1033 Flat image that had previously been applied to exposure.
1034 illumination_correction : `lsst.afw.image.Exposure`
1035 Illumination correction image to convert to photometric-flattened image.
1037 Returns
1038 -------
1039 background_to_photometric_ratio : `lsst.afw.image.Image`
1040 Ratio image to convert a photometric-flattened image to/from
1041 a background-flattened image. Will be None if task not
1042 configured to use the illumination correction.
1043 """
1044 if not self.config.do_illumination_correction:
1045 return None
1047 # From a raw image to a background-flattened image, we have:
1048 # bfi = image / background_flat
1049 # From a raw image to a photometric-flattened image, we have:
1050 # pfi = image / reference_flux_flat
1051 # pfi = image / (dome_flat * illumination_correction),
1052 # where the illumination correction contains the jacobian
1053 # of the wcs, converting to fluence units.
1054 # Currently background_flat == dome_flat, so we have for the
1055 # "background_to_photometric_ratio", the ratio of the background-
1056 # flattened image to the photometric-flattened image:
1057 # bfi / pfi = illumination_correction.
1059 background_to_photometric_ratio = illumination_correction.image.clone()
1061 # Dividing the ratio will convert a background-flattened image to
1062 # a photometric-flattened image.
1063 exposure.maskedImage /= background_to_photometric_ratio
1065 exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = True
1067 return background_to_photometric_ratio
1069 def _compute_psf(self, exposure, id_generator, background_to_photometric_ratio=None):
1070 """Find bright sources detected on an exposure and fit a PSF model to
1071 them, repairing likely cosmic rays before detection.
1073 Repair, detect, measure, and compute PSF twice, to ensure the PSF
1074 model does not include contributions from cosmic rays.
1076 Parameters
1077 ----------
1078 exposure : `lsst.afw.image.Exposure`
1079 Exposure to detect and measure bright stars on.
1080 id_generator : `lsst.meas.base.IdGenerator`
1081 Object that generates source IDs and provides random seeds.
1082 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1083 Image to convert photometric-flattened image to
1084 background-flattened image.
1086 Returns
1087 -------
1088 sources : `lsst.afw.table.SourceCatalog`
1089 Catalog of detected bright sources.
1090 background : `lsst.afw.math.BackgroundList`
1091 Background that was fit to the exposure during detection.
1092 cell_set : `lsst.afw.math.SpatialCellSet`
1093 PSF candidates returned by the psf determiner.
1094 """
1095 def log_psf(msg, addToMetadata=False):
1096 """Log the parameters of the psf and background, with a prepended
1097 message. There is also the option to add the PSF sigma to the task
1098 metadata.
1100 Parameters
1101 ----------
1102 msg : `str`
1103 Message to prepend the log info with.
1104 addToMetadata : `bool`, optional
1105 Whether to add the final psf sigma value to the task metadata
1106 (the default is False).
1107 """
1108 position = exposure.psf.getAveragePosition()
1109 sigma = exposure.psf.computeShape(position).getDeterminantRadius()
1110 dimensions = exposure.psf.computeImage(position).getDimensions()
1111 median_background = np.median(background.getImage().array)
1112 self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",
1113 msg, sigma, dimensions, median_background)
1114 if addToMetadata:
1115 self.metadata["final_psf_sigma"] = sigma
1117 self.log.info("First pass detection with Guassian PSF FWHM=%s pixels",
1118 self.config.install_simple_psf.fwhm)
1119 self.install_simple_psf.run(exposure=exposure)
1121 background = self.psf_subtract_background.run(
1122 exposure=exposure,
1123 backgroundToPhotometricRatio=background_to_photometric_ratio,
1124 ).background
1125 log_psf("Initial PSF:")
1126 self.psf_repair.run(exposure=exposure, keepCRs=True)
1128 table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())
1129 # Re-estimate the background during this detection step, so that
1130 # measurement uses the most accurate background-subtraction.
1131 detections = self.psf_detection.run(
1132 table=table,
1133 exposure=exposure,
1134 background=background,
1135 backgroundToPhotometricRatio=background_to_photometric_ratio,
1136 )
1137 self.metadata["initial_psf_positive_footprint_count"] = detections.numPos
1138 self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg
1139 self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks
1140 self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks
1141 self.psf_source_measurement.run(detections.sources, exposure)
1142 psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
1143 # Replace the initial PSF with something simpler for the second
1144 # repair/detect/measure/measure_psf step: this can help it converge.
1145 self.install_simple_psf.run(exposure=exposure)
1147 log_psf("Rerunning with simple PSF:")
1148 # TODO investigation: Should we only re-run repair here, to use the
1149 # new PSF? Maybe we *do* need to re-run measurement with PsfFlux, to
1150 # use the fitted PSF?
1151 # TODO investigation: do we need a separate measurement task here
1152 # for the post-psf_measure_psf step, since we only want to do PsfFlux
1153 # and GaussianFlux *after* we have a PSF? Maybe that's not relevant
1154 # once DM-39203 is merged?
1155 self.psf_repair.run(exposure=exposure, keepCRs=True)
1156 # Re-estimate the background during this detection step, so that
1157 # measurement uses the most accurate background-subtraction.
1158 detections = self.psf_detection.run(
1159 table=table,
1160 exposure=exposure,
1161 background=background,
1162 backgroundToPhotometricRatio=background_to_photometric_ratio,
1163 )
1164 self.metadata["simple_psf_positive_footprint_count"] = detections.numPos
1165 self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg
1166 self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks
1167 self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks
1168 self.psf_source_measurement.run(detections.sources, exposure)
1169 psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
1171 log_psf("Final PSF:", addToMetadata=True)
1173 # Final repair with final PSF, removing cosmic rays this time.
1174 self.psf_repair.run(exposure=exposure)
1175 # Final measurement with the CRs removed.
1176 self.psf_source_measurement.run(detections.sources, exposure)
1178 # PSF is set on exposure; candidates are returned to use for
1179 # calibration flux normalization and aperture corrections.
1180 return detections.sources, background, psf_result.cellSet
1182 def _measure_aperture_correction(self, exposure, bright_sources):
1183 """Measure and set the ApCorrMap on the Exposure, using
1184 previously-measured bright sources.
1186 This function first normalizes the calibration flux and then
1187 the full set of aperture corrections are measured relative
1188 to this normalized calibration flux.
1190 Parameters
1191 ----------
1192 exposure : `lsst.afw.image.Exposure`
1193 Exposure to set the ApCorrMap on.
1194 bright_sources : `lsst.afw.table.SourceCatalog`
1195 Catalog of detected bright sources; modified to include columns
1196 necessary for point source determination for the aperture correction
1197 calculation.
1198 """
1199 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(
1200 exposure=exposure,
1201 catalog=bright_sources,
1202 ).ap_corr_map
1204 ap_corr_map = self.measure_aperture_correction.run(exposure, bright_sources).apCorrMap
1206 # Need to merge the aperture correction map from the normalization.
1207 for key in norm_ap_corr_map:
1208 ap_corr_map[key] = norm_ap_corr_map[key]
1210 exposure.info.setApCorrMap(ap_corr_map)
1212 def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None):
1213 """Detect stars on an exposure that has a PSF model, and measure their
1214 PSF, circular aperture, compensated gaussian fluxes.
1216 Parameters
1217 ----------
1218 exposure : `lsst.afw.image.Exposure`
1219 Exposure to detect and measure stars on.
1220 background : `lsst.afw.math.BackgroundList`
1221 Background that was fit to the exposure during detection;
1222 modified in-place during subsequent detection.
1223 id_generator : `lsst.meas.base.IdGenerator`
1224 Object that generates source IDs and provides random seeds.
1225 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1226 Image to convert photometric-flattened image to
1227 background-flattened image.
1229 Returns
1230 -------
1231 stars : `SourceCatalog`
1232 Sources that are very likely to be stars, with a limited set of
1233 measurements performed on them.
1234 """
1235 table = afwTable.SourceTable.make(self.initial_stars_schema.schema,
1236 id_generator.make_table_id_factory())
1237 # Re-estimate the background during this detection step, so that
1238 # measurement uses the most accurate background-subtraction.
1239 detections = self.star_detection.run(
1240 table=table,
1241 exposure=exposure,
1242 background=background,
1243 backgroundToPhotometricRatio=background_to_photometric_ratio,
1244 )
1245 sources = detections.sources
1246 self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources)
1248 n_sky_sources = np.sum(sources["sky_source"])
1249 if (self.config.do_downsample_footprints
1250 and (len(sources) - n_sky_sources) > self.config.downsample_max_footprints):
1251 if exposure.info.id is None:
1252 self.log.warning("Exposure does not have a proper id; using 0 seed for downsample.")
1253 seed = 0
1254 else:
1255 seed = exposure.info.id & 0xFFFFFFFF
1257 gen = np.random.RandomState(seed)
1259 # Isolate the sky sources before downsampling.
1260 indices = np.arange(len(sources))[~sources["sky_source"]]
1261 indices = gen.choice(
1262 indices,
1263 size=self.config.downsample_max_footprints,
1264 replace=False,
1265 )
1266 skyIndices, = np.where(sources["sky_source"])
1267 indices = np.concatenate((indices, skyIndices))
1269 self.log.info("Downsampling from %d to %d non-sky-source footprints.", len(sources), len(indices))
1271 sel = np.zeros(len(sources), dtype=bool)
1272 sel[indices] = True
1273 sources = sources[sel]
1275 # TODO investigation: Could this deblender throw away blends of non-PSF sources?
1276 self.star_deblend.run(exposure=exposure, sources=sources)
1277 # The deblender may not produce a contiguous catalog; ensure
1278 # contiguity for subsequent tasks.
1279 if not sources.isContiguous():
1280 sources = sources.copy(deep=True)
1282 # Measure everything, and use those results to select only stars.
1283 self.star_measurement.run(sources, exposure)
1284 self.metadata["post_deblend_source_count"] = np.sum(~sources["sky_source"])
1285 self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"])
1286 self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"])
1288 # Run the normalization calibration flux task to apply the
1289 # normalization correction to create normalized
1290 # calibration fluxes.
1291 self.star_normalized_calibration_flux.run(exposure=exposure, catalog=sources)
1292 self.star_apply_aperture_correction.run(sources, exposure.apCorrMap)
1293 self.star_catalog_calculation.run(sources)
1294 self.star_set_primary_flags.run(sources)
1296 result = self.star_selector.run(sources)
1297 # The star selector may not produce a contiguous catalog.
1298 if not result.sourceCat.isContiguous():
1299 return result.sourceCat.copy(deep=True)
1300 else:
1301 return result.sourceCat
1303 def _match_psf_stars(self, psf_stars, stars):
1304 """Match calibration stars to psf stars, to identify which were psf
1305 candidates, and which were used or reserved during psf measurement
1306 and the astrometric fit.
1308 Parameters
1309 ----------
1310 psf_stars : `lsst.afw.table.SourceCatalog`
1311 PSF candidate stars that were sent to the psf determiner and
1312 used in the astrometric fit. Used to populate psf and astrometry
1313 related flag fields.
1314 stars : `lsst.afw.table.SourceCatalog`
1315 Stars that will be used for calibration; psf-related fields will
1316 be updated in-place.
1318 Notes
1319 -----
1320 This code was adapted from CalibrateTask.copyIcSourceFields().
1321 """
1322 control = afwTable.MatchControl()
1323 # Return all matched objects, to separate blends.
1324 control.findOnlyClosest = False
1325 matches = afwTable.matchXy(psf_stars, stars, 3.0, control)
1326 deblend_key = stars.schema["deblend_nChild"].asKey()
1327 matches = [m for m in matches if m[1].get(deblend_key) == 0]
1329 # Because we had to allow multiple matches to handle parents, we now
1330 # need to prune to the best (closest) matches.
1331 # Closest matches is a dict of psf_stars source ID to Match record
1332 # (psf_stars source, sourceCat source, distance in pixels).
1333 best = {}
1334 for match_psf, match_stars, d in matches:
1335 match = best.get(match_psf.getId())
1336 if match is None or d <= match[2]:
1337 best[match_psf.getId()] = (match_psf, match_stars, d)
1338 matches = list(best.values())
1339 # We'll use this to construct index arrays into each catalog.
1340 ids = np.array([(match_psf.getId(), match_stars.getId()) for match_psf, match_stars, d in matches]).T
1342 if (n_matches := len(matches)) == 0:
1343 raise NoPsfStarsToStarsMatchError(n_psf_stars=len(psf_stars), n_stars=len(stars))
1345 self.log.info("%d psf/astrometry stars out of %d matched %d calib stars",
1346 n_matches, len(psf_stars), len(stars))
1347 self.metadata["matched_psf_star_count"] = n_matches
1349 # Check that no stars sources are listed twice; we already know
1350 # that each match has a unique psf_stars id, due to using as the key
1351 # in best above.
1352 n_unique = len(set(m[1].getId() for m in matches))
1353 if n_unique != n_matches:
1354 self.log.warning("%d psf_stars matched only %d stars", n_matches, n_unique)
1356 # The indices of the IDs, so we can update the flag fields as arrays.
1357 idx_psf_stars = np.searchsorted(psf_stars["id"], ids[0])
1358 idx_stars = np.searchsorted(stars["id"], ids[1])
1359 for field in self.psf_fields:
1360 result = np.zeros(len(stars), dtype=bool)
1361 result[idx_stars] = psf_stars[field][idx_psf_stars]
1362 stars[field] = result
1363 stars['psf_id'][idx_stars] = psf_stars['id'][idx_psf_stars]
1364 stars['psf_max_value'][idx_stars] = psf_stars['psf_max_value'][idx_psf_stars]
1366 def _fit_astrometry(self, exposure, stars):
1367 """Fit an astrometric model to the data and return the reference
1368 matches used in the fit, and the fitted WCS.
1370 Parameters
1371 ----------
1372 exposure : `lsst.afw.image.Exposure`
1373 Exposure that is being fit, to get PSF and other metadata from.
1374 Modified to add the fitted skyWcs.
1375 stars : `SourceCatalog`
1376 Good stars selected for use in calibration, with RA/Dec coordinates
1377 computed from the pixel positions and fitted WCS.
1379 Returns
1380 -------
1381 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1382 Reference/stars matches used in the fit.
1383 """
1384 result = self.astrometry.run(stars, exposure)
1385 return result.matches, result.matchMeta
1387 def _fit_photometry(self, exposure, stars):
1388 """Fit a photometric model to the data and return the reference
1389 matches used in the fit, and the fitted PhotoCalib.
1391 Parameters
1392 ----------
1393 exposure : `lsst.afw.image.Exposure`
1394 Exposure that is being fit, to get PSF and other metadata from.
1395 Has the fit `lsst.afw.image.PhotoCalib` attached, with pixel values
1396 unchanged.
1397 stars : `lsst.afw.table.SourceCatalog`
1398 Good stars selected for use in calibration.
1399 background : `lsst.afw.math.BackgroundList`
1400 Background that was fit to the exposure during detection of the
1401 above stars.
1403 Returns
1404 -------
1405 calibrated_stars : `lsst.afw.table.SourceCatalog`
1406 Star catalog with flux/magnitude columns computed from the fitted
1407 photoCalib (instFlux columns are retained as well).
1408 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1409 Reference/stars matches used in the fit.
1410 matchMeta : `lsst.daf.base.PropertyList`
1411 Metadata needed to unpersist matches, as returned by the matcher.
1412 photo_calib : `lsst.afw.image.PhotoCalib`
1413 Photometric calibration that was fit to the star catalog.
1414 """
1415 result = self.photometry.run(exposure, stars)
1416 calibrated_stars = result.photoCalib.calibrateCatalog(stars)
1417 exposure.setPhotoCalib(result.photoCalib)
1418 return calibrated_stars, result.matches, result.matchMeta, result.photoCalib
1420 def _apply_photometry(self, exposure, background, background_to_photometric_ratio=None):
1421 """Apply the photometric model attached to the exposure to the
1422 exposure's pixels and an associated background model.
1424 Parameters
1425 ----------
1426 exposure : `lsst.afw.image.Exposure`
1427 Exposure with the target `lsst.afw.image.PhotoCalib` attached.
1428 On return, pixel values will be calibrated and an identity
1429 photometric transform will be attached.
1430 background : `lsst.afw.math.BackgroundList`
1431 Background model to convert to nanojansky units in place.
1432 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1433 Image to convert photometric-flattened image to
1434 background-flattened image.
1435 """
1436 photo_calib = exposure.getPhotoCalib()
1437 exposure.maskedImage = photo_calib.calibrateImage(exposure.maskedImage)
1438 identity = afwImage.PhotoCalib(1.0,
1439 photo_calib.getCalibrationErr(),
1440 bbox=exposure.getBBox())
1441 exposure.setPhotoCalib(identity)
1442 exposure.metadata["BUNIT"] = "nJy"
1444 assert photo_calib._isConstant, \
1445 "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that."
1447 for bg in background:
1448 # The statsImage is a view, but we can't assign to a function call in python.
1449 binned_image = bg[0].getStatsImage()
1450 binned_image *= photo_calib.getCalibrationMean()
1452 def _summarize(self, exposure, stars, background):
1453 """Compute summary statistics on the exposure and update in-place the
1454 calibrations attached to it.
1456 Parameters
1457 ----------
1458 exposure : `lsst.afw.image.Exposure`
1459 Exposure that was calibrated, to get PSF and other metadata from.
1460 Should be in instrumental units with the photometric calibration
1461 attached.
1462 Modified to contain the computed summary statistics.
1463 stars : `SourceCatalog`
1464 Good stars selected used in calibration.
1465 background : `lsst.afw.math.BackgroundList`
1466 Background that was fit to the exposure during detection of the
1467 above stars. Should be in instrumental units.
1468 """
1469 summary = self.compute_summary_stats.run(exposure, stars, background)
1470 exposure.info.setSummaryStats(summary)
1472 def _recordMaskedPixelFractions(self, exposure):
1473 """Record the fraction of all the pixels in an exposure
1474 that are masked with a given flag. Each fraction is
1475 recorded in the task metadata. One record per flag type.
1477 Parameters
1478 ----------
1479 exposure : `lsst.afw.image.ExposureF`
1480 The target exposure to calculate masked pixel fractions for.
1481 """
1483 mask = exposure.mask
1484 maskPlanes = list(mask.getMaskPlaneDict().keys())
1485 for maskPlane in maskPlanes:
1486 self.metadata[f"{maskPlane.lower()}_mask_fraction"] = (
1487 evaluateMaskFraction(mask, maskPlane)
1488 )
1490 def _update_wcs_with_camera_model(self, exposure, cameraModel):
1491 """Replace the existing WCS with one generated using the pointing
1492 and the input camera model.
1494 Parameters
1495 ----------
1496 exposure : `lsst.afw.image.ExposureF`
1497 The exposure whose WCS will be updated.
1498 cameraModel : `lsst.afw.cameraGeom.Camera`
1499 Camera to be used when constructing updated WCS.
1500 """
1501 detector = cameraModel[exposure.detector.getId()]
1502 newWcs = createInitialSkyWcs(exposure.visitInfo, detector)
1503 exposure.setWcs(newWcs)