LSST Applications g042eb84c57+730a74494b,g04e9c324dd+8c5ae1fdc5,g134cb467dc+1f1e3e7524,g199a45376c+0ba108daf9,g1fd858c14a+fa7d31856b,g210f2d0738+f66ac109ec,g262e1987ae+83a3acc0e5,g29ae962dfc+d856a2cb1f,g2cef7863aa+aef1011c0b,g35bb328faa+8c5ae1fdc5,g3fd5ace14f+a1e0c9f713,g47891489e3+0d594cb711,g4d44eb3520+c57ec8f3ed,g4d7b6aa1c5+f66ac109ec,g53246c7159+8c5ae1fdc5,g56a1a4eaf3+fd7ad03fde,g64539dfbff+f66ac109ec,g67b6fd64d1+0d594cb711,g67fd3c3899+f66ac109ec,g6985122a63+0d594cb711,g74acd417e5+3098891321,g786e29fd12+668abc6043,g81db2e9a8d+98e2ab9f28,g87389fa792+8856018cbb,g89139ef638+0d594cb711,g8d7436a09f+80fda9ce03,g8ea07a8fe4+760ca7c3fc,g90f42f885a+033b1d468d,g97be763408+a8a29bda4b,g99822b682c+e3ec3c61f9,g9d5c6a246b+0d5dac0c3d,ga41d0fce20+9243b26dd2,gbf99507273+8c5ae1fdc5,gd7ef33dd92+0d594cb711,gdab6d2f7ff+3098891321,ge410e46f29+0d594cb711,geaed405ab2+c4bbc419c6,gf9a733ac38+8c5ae1fdc5,w.2025.38
LSST Data Management Base Package
Loading...
Searching...
No Matches
calibrateImage.py
Go to the documentation of this file.
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/>.
21
22__all__ = ["CalibrateImageTask", "CalibrateImageConfig", "NoPsfStarsToStarsMatchError",
23 "AllCentroidsFlaggedError"]
24
25import numpy as np
26import requests
27import os
28
29import lsst.afw.table as afwTable
30import lsst.afw.image as afwImage
31from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache
36import lsst.meas.base
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
45
46from . import measurePsf, repair, photoCal, computeExposureSummaryStats, snapCombine
47
48
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
63
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 }
72
73
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
85
86 @property
87 def metadata(self):
88 return {"n_psf_stars": self.n_psf_stars,
89 "n_stars": self.n_stars
90 }
91
92
93class CalibrateImageConnections(pipeBase.PipelineTaskConnections,
94 dimensions=("instrument", "visit", "detector")):
95
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 )
127
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 )
142
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 )
149
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 )
195
196 # Optional outputs
197 psf_stars_footprints = connectionTypes.Output(
198 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination; "
199 "includes source footprints.",
200 name="initial_psf_stars_footprints_detector",
201 storageClass="SourceCatalog",
202 dimensions=["instrument", "visit", "detector"],
203 )
204 psf_stars = connectionTypes.Output(
205 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination.",
206 name="initial_psf_stars_detector",
207 storageClass="ArrowAstropy",
208 dimensions=["instrument", "visit", "detector"],
209 )
210 astrometry_matches = connectionTypes.Output(
211 doc="Source to reference catalog matches from the astrometry solver.",
212 name="initial_astrometry_match_detector",
213 storageClass="Catalog",
214 dimensions=("instrument", "visit", "detector"),
215 )
216 photometry_matches = connectionTypes.Output(
217 doc="Source to reference catalog matches from the photometry solver.",
218 name="initial_photometry_match_detector",
219 storageClass="Catalog",
220 dimensions=("instrument", "visit", "detector"),
221 )
222
223 def __init__(self, *, config=None):
224 super().__init__(config=config)
225 if "psf_stars" not in config.optional_outputs:
226 del self.psf_stars
227 if "psf_stars_footprints" not in config.optional_outputs:
228 del self.psf_stars_footprints
229 if "astrometry_matches" not in config.optional_outputs:
230 del self.astrometry_matches
231 if "photometry_matches" not in config.optional_outputs:
232 del self.photometry_matches
233 if not config.do_calibrate_pixels:
234 del self.applied_photo_calib
235 if not config.do_illumination_correction:
236 del self.background_flat
239 if not config.useButlerCamera:
240 del self.camera_model
241
242
243class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CalibrateImageConnections):
244 optional_outputs = pexConfig.ListField(
245 doc="Which optional outputs to save (as their connection name)?",
246 dtype=str,
247 # TODO: note somewhere to disable this for benchmarking, but should
248 # we always have it on for production runs?
249 default=["psf_stars", "psf_stars_footprints", "astrometry_matches", "photometry_matches"],
250 optional=False
251 )
252
253 # To generate catalog ids consistently across subtasks.
254 id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()
255
256 snap_combine = pexConfig.ConfigurableField(
258 doc="Task to combine two snaps to make one exposure.",
259 )
260
261 # subtasks used during psf characterization
262 install_simple_psf = pexConfig.ConfigurableField(
264 doc="Task to install a simple PSF model into the input exposure to use "
265 "when detecting bright sources for PSF estimation.",
266 )
267 psf_repair = pexConfig.ConfigurableField(
268 target=repair.RepairTask,
269 doc="Task to repair cosmic rays on the exposure before PSF determination.",
270 )
271 psf_subtract_background = pexConfig.ConfigurableField(
273 doc="Task to perform intial background subtraction, before first detection pass.",
274 )
275 psf_detection = pexConfig.ConfigurableField(
277 doc="Task to detect sources for PSF determination."
278 )
279 psf_source_measurement = pexConfig.ConfigurableField(
281 doc="Task to measure sources to be used for psf estimation."
282 )
283 psf_measure_psf = pexConfig.ConfigurableField(
285 doc="Task to measure the psf on bright sources."
286 )
287 psf_normalized_calibration_flux = pexConfig.ConfigurableField(
289 doc="Task to normalize the calibration flux (e.g. compensated tophats) "
290 "for the bright stars used for psf estimation.",
291 )
292
293 # TODO DM-39203: we can remove aperture correction from this task once we are
294 # using the shape-based star/galaxy code.
295 measure_aperture_correction = pexConfig.ConfigurableField(
297 doc="Task to compute the aperture correction from the bright stars."
298 )
299
300 # subtasks used during star measurement
301 star_detection = pexConfig.ConfigurableField(
303 doc="Task to detect stars to return in the output catalog."
304 )
305 star_sky_sources = pexConfig.ConfigurableField(
307 doc="Task to generate sky sources ('empty' regions where there are no detections).",
308 )
309 do_downsample_footprints = pexConfig.Field(
310 dtype=bool,
311 default=False,
312 doc="Downsample footprints prior to deblending to optimize speed?",
313 )
314 downsample_max_footprints = pexConfig.Field(
315 dtype=int,
316 default=1000,
317 doc="Maximum number of non-sky-source footprints to use if do_downsample_footprints is True,",
318 )
319 star_deblend = pexConfig.ConfigurableField(
321 doc="Split blended sources into their components."
322 )
323 star_measurement = pexConfig.ConfigurableField(
325 doc="Task to measure stars to return in the output catalog."
326 )
327 star_normalized_calibration_flux = pexConfig.ConfigurableField(
329 doc="Task to apply the normalization for calibration fluxes (e.g. compensated tophats) "
330 "for the final output star catalog.",
331 )
332 star_apply_aperture_correction = pexConfig.ConfigurableField(
334 doc="Task to apply aperture corrections to the selected stars."
335 )
336 star_catalog_calculation = pexConfig.ConfigurableField(
338 doc="Task to compute extendedness values on the star catalog, "
339 "for the star selector to remove extended sources."
340 )
341 star_set_primary_flags = pexConfig.ConfigurableField(
343 doc="Task to add isPrimary to the catalog."
344 )
345 star_selector = lsst.meas.algorithms.sourceSelectorRegistry.makeField(
346 default="science",
347 doc="Task to select reliable stars to use for calibration."
348 )
349
350 # final calibrations and statistics
351 astrometry = pexConfig.ConfigurableField(
353 doc="Task to perform astrometric calibration to fit a WCS.",
354 )
355 astrometry_ref_loader = pexConfig.ConfigField(
357 doc="Configuration of reference object loader for astrometric fit.",
358 )
359 photometry = pexConfig.ConfigurableField(
361 doc="Task to perform photometric calibration to fit a PhotoCalib.",
362 )
363 photometry_ref_loader = pexConfig.ConfigField(
365 doc="Configuration of reference object loader for photometric fit.",
366 )
367
368 compute_summary_stats = pexConfig.ConfigurableField(
370 doc="Task to to compute summary statistics on the calibrated exposure."
371 )
372
373 do_illumination_correction = pexConfig.Field(
374 dtype=bool,
375 default=False,
376 doc="If True, apply the illumination correction. This assumes that the "
377 "input image has already been flat-fielded such that it is suitable "
378 "for background subtraction.",
379 )
380
381 do_calibrate_pixels = pexConfig.Field(
382 dtype=bool,
383 default=True,
384 doc=(
385 "If True, apply the photometric calibration to the image pixels "
386 "and background model, and attach an identity PhotoCalib to "
387 "the output image to reflect this. If False`, leave the image "
388 "and background uncalibrated and attach the PhotoCalib that maps "
389 "them to physical units."
390 )
391 )
392
393 do_include_astrometric_errors = pexConfig.Field(
394 dtype=bool,
395 default=True,
396 doc="If True, include astrometric errors in the output catalog.",
397 )
398
399 run_sattle = pexConfig.Field(
400 dtype=bool,
401 default=False,
402 doc="If True, the sattle service will populate a cache for later use "
403 "in ip_diffim.detectAndMeasure alert verification."
404 )
405
406 sattle_historical = pexConfig.Field(
407 dtype=bool,
408 default=False,
409 doc="If re-running a pipeline that requires sattle, this should be set "
410 "to True. This will populate sattle's cache with the historic data "
411 "closest in time to the exposure.",
412 )
413
414 useButlerCamera = pexConfig.Field(
415 dtype=bool,
416 default=False,
417 doc="If True, use a camera distortion model generated elsewhere in "
418 "the pipeline combined with the telescope boresight as a starting "
419 "point for fitting the WCS, instead of using the WCS attached to "
420 "the exposure, which is generated from the boresight and the "
421 "camera model from the obs_* package."
422 )
423
424 def setDefaults(self):
425 super().setDefaults()
426
427 # Use a very broad PSF here, to throughly reject CRs.
428 # TODO investigation: a large initial psf guess may make stars look
429 # like CRs for very good seeing images.
430 self.install_simple_psf.fwhm = 4
431
432 # S/N>=50 sources for PSF determination, but detection to S/N=10.
433 # The thresholdValue sets the minimum flux in a pixel to be included in the
434 # footprint, while peaks are only detected when they are above
435 # thresholdValue * includeThresholdMultiplier. The low thresholdValue
436 # ensures that the footprints are large enough for the noise replacer
437 # to mask out faint undetected neighbors that are not to be measured.
438 self.psf_detection.thresholdValue = 10.0
439 self.psf_detection.includeThresholdMultiplier = 5.0
440 # TODO investigation: Probably want False here, but that may require
441 # tweaking the background spatial scale, to make it small enough to
442 # prevent extra peaks in the wings of bright objects.
443 self.psf_detection.doTempLocalBackground = False
444 # NOTE: we do want reEstimateBackground=True in psf_detection, so that
445 # each measurement step is done with the best background available.
446
447 # Minimal measurement plugins for PSF determination.
448 # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use
449 # shapeHSM/moments for star/galaxy separation.
450 # TODO DM-39203: we can remove aperture correction from this task once
451 # we are using the shape-based star/galaxy code.
452 self.psf_source_measurement.plugins = ["base_PixelFlags",
453 "base_SdssCentroid",
454 "ext_shapeHSM_HsmSourceMoments",
455 "base_CircularApertureFlux",
456 "base_GaussianFlux",
457 "base_PsfFlux",
458 "base_CompensatedTophatFlux",
459 ]
460 self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
461 # Only measure apertures we need for PSF measurement.
462 self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
463 self.psf_source_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
464 # TODO DM-40843: Remove this line once this is the psfex default.
465 self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \
466 "base_CircularApertureFlux_12_0_instFlux"
467
468 # No extendeness information available: we need the aperture
469 # corrections to determine that.
470 self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False
471 self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]
472 self.measure_aperture_correction.sourceSelector["science"].flags.bad = []
473
474 # Detection for good S/N for astrometry/photometry and other
475 # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks.
476 self.star_detection.thresholdValue = 5.0
477 self.star_detection.includeThresholdMultiplier = 2.0
478 self.star_measurement.plugins = ["base_PixelFlags",
479 "base_SdssCentroid",
480 "ext_shapeHSM_HsmSourceMoments",
481 "ext_shapeHSM_HsmPsfMoments",
482 "base_GaussianFlux",
483 "base_PsfFlux",
484 "base_CircularApertureFlux",
485 "base_ClassificationSizeExtendedness",
486 "base_CompensatedTophatFlux",
487 ]
488 self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"
489 self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"
490 # Only measure the apertures we need for star selection.
491 self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]
492 self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]
493
494 # We measure and apply the normalization aperture correction with the
495 # psf_normalized_calibration_flux task, and we only apply the normalization
496 # aperture correction for the full list of stars.
497 self.star_normalized_calibration_flux.do_measure_ap_corr = False
498
499 # Select stars with reliable measurements and no bad flags.
500 self.star_selector["science"].doFlags = True
501 self.star_selector["science"].doUnresolved = True
502 self.star_selector["science"].doSignalToNoise = True
503 self.star_selector["science"].signalToNoise.minimum = 10.0
504 # Keep sky sources in the output catalog, even though they aren't
505 # wanted for calibration.
506 self.star_selector["science"].doSkySources = True
507 # Set the flux and error fields
508 self.star_selector["science"].signalToNoise.fluxField = "slot_CalibFlux_instFlux"
509 self.star_selector["science"].signalToNoise.errField = "slot_CalibFlux_instFluxErr"
510
511 # Use the affine WCS fitter (assumes we have a good camera geometry).
512 self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)
513 # phot_g_mean is the primary Gaia band for all input bands.
514 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"
515
516 # Only reject sky sources; we already selected good stars.
517 self.astrometry.sourceSelector["science"].doFlags = True
518 self.astrometry.sourceSelector["science"].flags.good = ["calib_psf_candidate"]
519 self.astrometry.sourceSelector["science"].flags.bad = []
520 self.astrometry.sourceSelector["science"].doUnresolved = False
521 self.astrometry.sourceSelector["science"].doIsolated = False
522 self.astrometry.sourceSelector["science"].doRequirePrimary = False
523 self.photometry.match.sourceSelection.doFlags = True
524 self.photometry.match.sourceSelection.flags.bad = ["sky_source"]
525 # Unset the (otherwise reasonable, but we've already made the
526 # selections we want above) selection settings in PhotoCalTask.
527 self.photometry.match.sourceSelection.doRequirePrimary = False
528 self.photometry.match.sourceSelection.doUnresolved = False
529
530 def validate(self):
531 super().validate()
532
533 # Ensure that the normalization calibration flux tasks
534 # are configured correctly.
535 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:
536 msg = ("psf_normalized_calibration_flux task must be configured with do_measure_ap_corr=True "
537 "or else the normalization and calibration flux will not be properly measured.")
538 raise pexConfig.FieldValidationError(
539 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,
540 )
541 if self.star_normalized_calibration_flux.do_measure_ap_corr:
542 msg = ("star_normalized_calibration_flux task must be configured with do_measure_ap_corr=False "
543 "to apply the previously measured normalization to the full catalog of calibration "
544 "fluxes.")
545 raise pexConfig.FieldValidationError(
546 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,
547 )
548
549 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,
550 # because they'd be running too early to pick up the fitted PhotoCalib
551 # and WCS.
552 if "base_LocalWcs" in self.psf_source_measurement.plugins.names:
553 raise pexConfig.FieldValidationError(
554 CalibrateImageConfig.psf_source_measurement,
555 self,
556 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
557 )
558 if "base_LocalWcs" in self.star_measurement.plugins.names:
559 raise pexConfig.FieldValidationError(
560 CalibrateImageConfig.star_measurement,
561 self,
562 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."
563 )
564 if "base_LocalPhotoCalib" in self.psf_source_measurement.plugins.names:
565 raise pexConfig.FieldValidationError(
566 CalibrateImageConfig.psf_source_measurement,
567 self,
568 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
569 "as it would be run before the photometry fit."
570 )
571 if "base_LocalPhotoCalib" in self.star_measurement.plugins.names:
572 raise pexConfig.FieldValidationError(
573 CalibrateImageConfig.star_measurement,
574 self,
575 "base_LocalPhotoCalib cannot run CalibrateImageTask, "
576 "as it would be run before the photometry fit."
577 )
578
579 # Check for illumination correction and background consistency.
581 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:
582 raise pexConfig.FieldValidationError(
583 CalibrateImageConfig.psf_subtract_background,
584 self,
585 "CalibrateImageTask.psf_subtract_background must be configured with "
586 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
587 )
588 if self.psf_detection.reEstimateBackground:
589 if not self.psf_detection.doApplyFlatBackgroundRatio:
590 raise pexConfig.FieldValidationError(
591 CalibrateImageConfig.psf_detection,
592 self,
593 "CalibrateImageTask.psf_detection background must be configured with "
594 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
595 )
596 if self.star_detection.reEstimateBackground:
597 if not self.star_detection.doApplyFlatBackgroundRatio:
598 raise pexConfig.FieldValidationError(
599 CalibrateImageConfig.star_detection,
600 self,
601 "CalibrateImageTask.star_detection background must be configured with "
602 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."
603 )
604
605 if self.run_sattle:
606 if not os.getenv("SATTLE_URI_BASE"):
607 raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self,
608 "Sattle requested but URI environment variable not set.")
609
610
611class CalibrateImageTask(pipeBase.PipelineTask):
612 """Compute the PSF, aperture corrections, astrometric and photometric
613 calibrations, and summary statistics for a single science exposure, and
614 produce a catalog of brighter stars that were used to calibrate it.
615
616 Parameters
617 ----------
618 initial_stars_schema : `lsst.afw.table.Schema`
619 Schema of the initial_stars output catalog.
620 """
621 _DefaultName = "calibrateImage"
622 ConfigClass = CalibrateImageConfig
623
624 def __init__(self, initial_stars_schema=None, **kwargs):
625 super().__init__(**kwargs)
626 self.makeSubtask("snap_combine")
627
628 # PSF determination subtasks
629 self.makeSubtask("install_simple_psf")
630 self.makeSubtask("psf_repair")
631 self.makeSubtask("psf_subtract_background")
632 self.psf_schema = afwTable.SourceTable.makeMinimalSchema()
633 self.psf_schema.addField(
634 'psf_max_value',
635 type=np.float32,
636 doc="PSF max value.",
637 )
638 afwTable.CoordKey.addErrorFields(self.psf_schema)
639 self.makeSubtask("psf_detection", schema=self.psf_schema)
640 self.makeSubtask("psf_source_measurement", schema=self.psf_schema)
641 self.makeSubtask("psf_measure_psf", schema=self.psf_schema)
642 self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema)
643
644 self.makeSubtask("measure_aperture_correction", schema=self.psf_schema)
645 self.makeSubtask("astrometry", schema=self.psf_schema)
646
647 # star measurement subtasks
648 if initial_stars_schema is None:
649 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()
650
651 # These fields let us track which sources were used for psf modeling,
652 # astrometric fitting, and aperture correction calculations.
653 self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved",
654 "calib_astrometry_used",
655 # TODO DM-39203: these can be removed once apcorr is gone.
656 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",
657 "apcorr_base_PsfFlux_used",)
658 for field in self.psf_fields:
659 item = self.psf_schema.find(field)
660 initial_stars_schema.addField(item.getField())
661 id_type = self.psf_schema["id"].asField().getTypeString()
662 psf_max_value_type = self.psf_schema['psf_max_value'].asField().getTypeString()
663 initial_stars_schema.addField("psf_id",
664 type=id_type,
665 doc="id of this source in psf_stars; 0 if there is no match.")
666 initial_stars_schema.addField("psf_max_value",
667 type=psf_max_value_type,
668 doc="Maximum value in the star image used to train PSF.")
669
670 afwTable.CoordKey.addErrorFields(initial_stars_schema)
671 self.makeSubtask("star_detection", schema=initial_stars_schema)
672 self.makeSubtask("star_sky_sources", schema=initial_stars_schema)
673 self.makeSubtask("star_deblend", schema=initial_stars_schema)
674 self.makeSubtask("star_measurement", schema=initial_stars_schema)
675 self.makeSubtask("star_normalized_calibration_flux", schema=initial_stars_schema)
676
677 self.makeSubtask("star_apply_aperture_correction", schema=initial_stars_schema)
678 self.makeSubtask("star_catalog_calculation", schema=initial_stars_schema)
679 self.makeSubtask("star_set_primary_flags", schema=initial_stars_schema, isSingleFrame=True)
680 self.makeSubtask("star_selector")
681 self.makeSubtask("photometry", schema=initial_stars_schema)
682 self.makeSubtask("compute_summary_stats")
683
684 # The final catalog will have calibrated flux columns, which we add to
685 # the init-output schema by calibrating our zero-length catalog with an
686 # arbitrary dummy PhotoCalib. We also use this schema to initialze
687 # the stars catalog in order to ensure it's the same even when we hit
688 # an error (and write partial outputs) before calibrating the catalog
689 # - note that calibrateCatalog will happily reuse existing output
690 # columns.
691 dummy_photo_calib = afwImage.PhotoCalib(1.0, 0, bbox=lsst.geom.Box2I())
692 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(
693 afwTable.SourceCatalog(initial_stars_schema)
694 )
695
696 def runQuantum(self, butlerQC, inputRefs, outputRefs):
697 inputs = butlerQC.get(inputRefs)
698 exposures = inputs.pop("exposures")
699
700 id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)
701
703 dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],
704 refCats=inputs.pop("astrometry_ref_cat"),
705 name=self.config.connections.astrometry_ref_cat,
706 config=self.config.astrometry_ref_loader, log=self.log)
707 self.astrometry.setRefObjLoader(astrometry_loader)
708
710 dataIds=[ref.datasetRef.dataId for ref in inputRefs.photometry_ref_cat],
711 refCats=inputs.pop("photometry_ref_cat"),
712 name=self.config.connections.photometry_ref_cat,
713 config=self.config.photometry_ref_loader, log=self.log)
714 self.photometry.match.setRefObjLoader(photometry_loader)
715
716 if self.config.do_illumination_correction:
717 background_flat = inputs.pop("background_flat")
718 illumination_correction = inputs.pop("illumination_correction")
719 else:
720 background_flat = None
721 illumination_correction = None
722
723 camera_model = None
724 if self.config.useButlerCamera:
725 if "camera_model" in inputs:
726 camera_model = inputs.pop("camera_model")
727 else:
728 self.log.warning("useButlerCamera=True, but camera is not available for filter %s. The "
729 "astrometry fit will use the WCS already attached to the exposure.",
730 exposures[0].filter.bandLabel)
731
732 # This should not happen with a properly configured execution context.
733 assert not inputs, "runQuantum got more inputs than expected"
734
735 # Specify the fields that `annotate` needs below, to ensure they
736 # exist, even as None.
737 result = pipeBase.Struct(
738 exposure=None,
739 stars_footprints=None,
740 psf_stars_footprints=None,
741 background_to_photometric_ratio=None,
742 )
743 try:
744 self.run(
745 exposures=exposures,
746 result=result,
747 id_generator=id_generator,
748 background_flat=background_flat,
749 illumination_correction=illumination_correction,
750 camera_model=camera_model,
751 )
752 except pipeBase.AlgorithmError as e:
753 error = pipeBase.AnnotatedPartialOutputsError.annotate(
754 e,
755 self,
756 result.exposure,
757 result.psf_stars_footprints,
758 result.stars_footprints,
759 log=self.log
760 )
761 butlerQC.put(result, outputRefs)
762 raise error from e
763
764 butlerQC.put(result, outputRefs)
765
766 @timeMethod
767 def run(
768 self,
769 *,
770 exposures,
771 id_generator=None,
772 result=None,
773 background_flat=None,
774 illumination_correction=None,
775 camera_model=None,
776 ):
777 """Find stars and perform psf measurement, then do a deeper detection
778 and measurement and calibrate astrometry and photometry from that.
779
780 Parameters
781 ----------
782 exposures : `lsst.afw.image.Exposure` or `list` [`lsst.afw.image.Exposure`]
783 Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter.
784 Modified in-place during processing if only one is passed.
785 If two exposures are passed, treat them as snaps and combine
786 before doing further processing.
787 id_generator : `lsst.meas.base.IdGenerator`, optional
788 Object that generates source IDs and provides random seeds.
789 result : `lsst.pipe.base.Struct`, optional
790 Result struct that is modified to allow saving of partial outputs
791 for some failure conditions. If the task completes successfully,
792 this is also returned.
793 background_flat : `lsst.afw.image.Exposure`, optional
794 Background flat-field image.
795 illumination_correction : `lsst.afw.image.Exposure`, optional
796 Illumination correction image.
797 camera_model : `lsst.afw.cameraGeom.Camera`, optional
798 Camera to be used if constructing updated WCS.
799
800 Returns
801 -------
802 result : `lsst.pipe.base.Struct`
803 Results as a struct with attributes:
804
805 ``exposure``
806 Calibrated exposure, with pixels in nJy units.
807 (`lsst.afw.image.Exposure`)
808 ``stars``
809 Stars that were used to calibrate the exposure, with
810 calibrated fluxes and magnitudes.
811 (`astropy.table.Table`)
812 ``stars_footprints``
813 Footprints of stars that were used to calibrate the exposure.
814 (`lsst.afw.table.SourceCatalog`)
815 ``psf_stars``
816 Stars that were used to determine the image PSF.
817 (`astropy.table.Table`)
818 ``psf_stars_footprints``
819 Footprints of stars that were used to determine the image PSF.
820 (`lsst.afw.table.SourceCatalog`)
821 ``background``
822 Background that was fit to the exposure when detecting
823 ``stars``. (`lsst.afw.math.BackgroundList`)
824 ``applied_photo_calib``
825 Photometric calibration that was fit to the star catalog and
826 applied to the exposure. (`lsst.afw.image.PhotoCalib`)
827 This is `None` if ``config.do_calibrate_pixels`` is `False`.
828 ``astrometry_matches``
829 Reference catalog stars matches used in the astrometric fit.
830 (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`)
831 ``photometry_matches``
832 Reference catalog stars matches used in the photometric fit.
833 (`list` [`lsst.afw.table.ReferenceMatch`] or `lsst.afw.table.BaseCatalog`)
834 """
835 if result is None:
836 result = pipeBase.Struct()
837 if id_generator is None:
838 id_generator = lsst.meas.base.IdGenerator()
839
840 result.exposure = self.snap_combine.run(exposures).exposure
841 self._recordMaskedPixelFractions(result.exposure)
842 self.log.info("Initial PhotoCalib: %s", result.exposure.getPhotoCalib())
843
844 result.exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = False
845
846 # Check input image processing.
847 if self.config.do_illumination_correction:
848 if not result.exposure.metadata.get("LSST ISR FLAT APPLIED", False):
849 raise pipeBase.InvalidQuantumError(
850 "Cannot use do_illumination_correction with an image that has not had a flat applied",
851 )
852
853 if camera_model:
854 if camera_model.get(result.exposure.detector.getId()):
855 self.log.info("Updating WCS with the provided camera model.")
856 self._update_wcs_with_camera_model(result.exposure, camera_model)
857 else:
858 self.log.warning(
859 "useButlerCamera=True, but detector %s is not available in the provided camera. The "
860 "astrometry fit will use the WCS already attached to the exposure.",
861 result.exposure.detector.getId())
862
863 result.background = None
864 summary_stat_catalog = None
865 # Some exposure components are set to initial placeholder objects
866 # while we try to bootstrap them. If we fail before we fit for them,
867 # we want to reset those components to None so the placeholders don't
868 # masquerade as the real thing.
869 have_fit_psf = False
870 have_fit_astrometry = False
871 have_fit_photometry = False
872 try:
873 result.background_to_photometric_ratio = self._apply_illumination_correction(
874 result.exposure,
875 background_flat,
876 illumination_correction,
877 )
878
879 result.psf_stars_footprints, result.background, _ = self._compute_psf(
880 result.exposure,
881 id_generator,
882 background_to_photometric_ratio=result.background_to_photometric_ratio,
883 )
884 have_fit_psf = True
885
886 # Check if all centroids have been flagged. This should happen
887 # after the PSF is fit and recorded because even a junky PSF
888 # model is informative.
889 if result.psf_stars_footprints["slot_Centroid_flag"].all():
890 psf_shape = result.exposure.psf.computeShape(result.exposure.psf.getAveragePosition())
892 n_sources=len(result.psf_stars_footprints),
893 psf_shape_ixx=psf_shape.getIxx(),
894 psf_shape_iyy=psf_shape.getIyy(),
895 psf_shape_ixy=psf_shape.getIxy(),
896 psf_size=psf_shape.getDeterminantRadius(),
897 )
898
899 self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)
900 result.psf_stars = result.psf_stars_footprints.asAstropy()
901 # Run astrometry using PSF candidate stars
902 astrometry_matches, astrometry_meta = self._fit_astrometry(
903 result.exposure, result.psf_stars_footprints
904 )
905 self.metadata["astrometry_matches_count"] = len(astrometry_matches)
906 if "astrometry_matches" in self.config.optional_outputs:
907 result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches,
908 astrometry_meta)
909 result.psf_stars = result.psf_stars_footprints.asAstropy()
910
911 # Run the stars_detection subtask for the photometric calibration.
912 result.stars_footprints = self._find_stars(
913 result.exposure,
914 result.background,
915 id_generator,
916 background_to_photometric_ratio=result.background_to_photometric_ratio,
917 )
918 self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints)
919
920 # Update the source cooordinates with the current wcs.
922 result.exposure.wcs,
923 sourceList=result.stars_footprints,
924 include_covariance=self.config.do_include_astrometric_errors
925 )
926
927 summary_stat_catalog = result.stars_footprints
928 result.stars = result.stars_footprints.asAstropy()
929 self.metadata["star_count"] = np.sum(~result.stars["sky_source"])
930
931 # Validate the astrometric fit. Send in the stars_footprints
932 # catalog so that its coords get set to NaN if the fit is deemed
933 # a failure.
934 self.astrometry.check(result.exposure, result.stars_footprints, len(astrometry_matches))
935 result.stars = result.stars_footprints.asAstropy()
936 have_fit_astrometry = True
937
938 result.stars_footprints, photometry_matches, \
939 photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints)
940 have_fit_photometry = True
941 self.metadata["photometry_matches_count"] = len(photometry_matches)
942 # fit_photometry returns a new catalog, so we need a new astropy table view.
943 result.stars = result.stars_footprints.asAstropy()
944 # summary stats don't make use of the calibrated fluxes, but we
945 # might as well use the best catalog we've got in case that
946 # changes, and help the old one get garbage-collected.
947 summary_stat_catalog = result.stars_footprints
948 if "photometry_matches" in self.config.optional_outputs:
949 result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,
950 photometry_meta)
951 except pipeBase.AlgorithmError:
952 if not have_fit_psf:
953 result.exposure.setPsf(None)
954 if not have_fit_astrometry:
955 result.exposure.setWcs(None)
956 if not have_fit_photometry:
957 result.exposure.setPhotoCalib(None)
958 # Summary stat calculations can handle missing components gracefully,
959 # but we want to run them as late as possible (but still before we
960 # calibrate pixels, if we do that at all).
961 # So we run them after we succeed or if we get an AlgorithmError. We
962 # intentionally don't use 'finally' here because we don't want to run
963 # them if we get some other kind of error.
964 self._summarize(result.exposure, summary_stat_catalog, result.background)
965 raise
966 else:
967 self._summarize(result.exposure, summary_stat_catalog, result.background)
968
969 if self.config.do_calibrate_pixels:
971 result.exposure,
972 result.background,
973 background_to_photometric_ratio=result.background_to_photometric_ratio,
974 )
975 result.applied_photo_calib = photo_calib
976 else:
977 result.applied_photo_calib = None
978
979 if self.config.run_sattle:
980 # send boresight and timing information to sattle so the cache
981 # is populated by the time we reach ip_diffim detectAndMeasure.
982 try:
983 populate_sattle_visit_cache(result.exposure.getInfo().getVisitInfo(),
984 historical=self.config.sattle_historical)
985 self.log.info('Successfully triggered load of sattle visit cache')
986 except requests.exceptions.HTTPError:
987 self.log.exception("Sattle visit cache update failed; continuing with image processing")
988
989 return result
990
991 def _apply_illumination_correction(self, exposure, background_flat, illumination_correction):
992 """Apply the illumination correction to a background-flattened image.
993
994 Parameters
995 ----------
996 exposure : `lsst.afw.image.Exposure`
997 Exposure to convert to a photometric-flattened image.
998 background_flat : `lsst.afw.image.Exposure`
999 Flat image that had previously been applied to exposure.
1000 illumination_correction : `lsst.afw.image.Exposure`
1001 Illumination correction image to convert to photometric-flattened image.
1002
1003 Returns
1004 -------
1005 background_to_photometric_ratio : `lsst.afw.image.Image`
1006 Ratio image to convert a photometric-flattened image to/from
1007 a background-flattened image. Will be None if task not
1008 configured to use the illumination correction.
1009 """
1010 if not self.config.do_illumination_correction:
1011 return None
1012
1013 # From a raw image to a background-flattened image, we have:
1014 # bfi = image / background_flat
1015 # From a raw image to a photometric-flattened image, we have:
1016 # pfi = image / reference_flux_flat
1017 # pfi = image / (dome_flat * illumination_correction),
1018 # where the illumination correction contains the jacobian
1019 # of the wcs, converting to fluence units.
1020 # Currently background_flat == dome_flat, so we have for the
1021 # "background_to_photometric_ratio", the ratio of the background-
1022 # flattened image to the photometric-flattened image:
1023 # bfi / pfi = illumination_correction.
1024
1025 background_to_photometric_ratio = illumination_correction.image.clone()
1026
1027 # Dividing the ratio will convert a background-flattened image to
1028 # a photometric-flattened image.
1029 exposure.maskedImage /= background_to_photometric_ratio
1030
1031 exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = True
1032
1033 return background_to_photometric_ratio
1034
1035 def _compute_psf(self, exposure, id_generator, background_to_photometric_ratio=None):
1036 """Find bright sources detected on an exposure and fit a PSF model to
1037 them, repairing likely cosmic rays before detection.
1038
1039 Repair, detect, measure, and compute PSF twice, to ensure the PSF
1040 model does not include contributions from cosmic rays.
1041
1042 Parameters
1043 ----------
1044 exposure : `lsst.afw.image.Exposure`
1045 Exposure to detect and measure bright stars on.
1046 id_generator : `lsst.meas.base.IdGenerator`
1047 Object that generates source IDs and provides random seeds.
1048 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1049 Image to convert photometric-flattened image to
1050 background-flattened image.
1051
1052 Returns
1053 -------
1054 sources : `lsst.afw.table.SourceCatalog`
1055 Catalog of detected bright sources.
1056 background : `lsst.afw.math.BackgroundList`
1057 Background that was fit to the exposure during detection.
1058 cell_set : `lsst.afw.math.SpatialCellSet`
1059 PSF candidates returned by the psf determiner.
1060 """
1061 def log_psf(msg, addToMetadata=False):
1062 """Log the parameters of the psf and background, with a prepended
1063 message. There is also the option to add the PSF sigma to the task
1064 metadata.
1065
1066 Parameters
1067 ----------
1068 msg : `str`
1069 Message to prepend the log info with.
1070 addToMetadata : `bool`, optional
1071 Whether to add the final psf sigma value to the task metadata
1072 (the default is False).
1073 """
1074 position = exposure.psf.getAveragePosition()
1075 sigma = exposure.psf.computeShape(position).getDeterminantRadius()
1076 dimensions = exposure.psf.computeImage(position).getDimensions()
1077 median_background = np.median(background.getImage().array)
1078 self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",
1079 msg, sigma, dimensions, median_background)
1080 if addToMetadata:
1081 self.metadata["final_psf_sigma"] = sigma
1082
1083 self.log.info("First pass detection with Guassian PSF FWHM=%s pixels",
1084 self.config.install_simple_psf.fwhm)
1085 self.install_simple_psf.run(exposure=exposure)
1086
1087 background = self.psf_subtract_background.run(
1088 exposure=exposure,
1089 backgroundToPhotometricRatio=background_to_photometric_ratio,
1090 ).background
1091 log_psf("Initial PSF:")
1092 self.psf_repair.run(exposure=exposure, keepCRs=True)
1093
1094 table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())
1095 # Re-estimate the background during this detection step, so that
1096 # measurement uses the most accurate background-subtraction.
1097 detections = self.psf_detection.run(
1098 table=table,
1099 exposure=exposure,
1100 background=background,
1101 backgroundToPhotometricRatio=background_to_photometric_ratio,
1102 )
1103 self.metadata["initial_psf_positive_footprint_count"] = detections.numPos
1104 self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg
1105 self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks
1106 self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks
1107 self.psf_source_measurement.run(detections.sources, exposure)
1108 psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
1109 # Replace the initial PSF with something simpler for the second
1110 # repair/detect/measure/measure_psf step: this can help it converge.
1111 self.install_simple_psf.run(exposure=exposure)
1112
1113 log_psf("Rerunning with simple PSF:")
1114 # TODO investigation: Should we only re-run repair here, to use the
1115 # new PSF? Maybe we *do* need to re-run measurement with PsfFlux, to
1116 # use the fitted PSF?
1117 # TODO investigation: do we need a separate measurement task here
1118 # for the post-psf_measure_psf step, since we only want to do PsfFlux
1119 # and GaussianFlux *after* we have a PSF? Maybe that's not relevant
1120 # once DM-39203 is merged?
1121 self.psf_repair.run(exposure=exposure, keepCRs=True)
1122 # Re-estimate the background during this detection step, so that
1123 # measurement uses the most accurate background-subtraction.
1124 detections = self.psf_detection.run(
1125 table=table,
1126 exposure=exposure,
1127 background=background,
1128 backgroundToPhotometricRatio=background_to_photometric_ratio,
1129 )
1130 self.metadata["simple_psf_positive_footprint_count"] = detections.numPos
1131 self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg
1132 self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks
1133 self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks
1134 self.psf_source_measurement.run(detections.sources, exposure)
1135 psf_result = self.psf_measure_psf.run(exposure=exposure, sources=detections.sources)
1136
1137 log_psf("Final PSF:", addToMetadata=True)
1138
1139 # Final repair with final PSF, removing cosmic rays this time.
1140 self.psf_repair.run(exposure=exposure)
1141 # Final measurement with the CRs removed.
1142 self.psf_source_measurement.run(detections.sources, exposure)
1143
1144 # PSF is set on exposure; candidates are returned to use for
1145 # calibration flux normalization and aperture corrections.
1146 return detections.sources, background, psf_result.cellSet
1147
1148 def _measure_aperture_correction(self, exposure, bright_sources):
1149 """Measure and set the ApCorrMap on the Exposure, using
1150 previously-measured bright sources.
1151
1152 This function first normalizes the calibration flux and then
1153 the full set of aperture corrections are measured relative
1154 to this normalized calibration flux.
1155
1156 Parameters
1157 ----------
1158 exposure : `lsst.afw.image.Exposure`
1159 Exposure to set the ApCorrMap on.
1160 bright_sources : `lsst.afw.table.SourceCatalog`
1161 Catalog of detected bright sources; modified to include columns
1162 necessary for point source determination for the aperture correction
1163 calculation.
1164 """
1165 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(
1166 exposure=exposure,
1167 catalog=bright_sources,
1168 ).ap_corr_map
1169
1170 ap_corr_map = self.measure_aperture_correction.run(exposure, bright_sources).apCorrMap
1171
1172 # Need to merge the aperture correction map from the normalization.
1173 for key in norm_ap_corr_map:
1174 ap_corr_map[key] = norm_ap_corr_map[key]
1175
1176 exposure.info.setApCorrMap(ap_corr_map)
1177
1178 def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None):
1179 """Detect stars on an exposure that has a PSF model, and measure their
1180 PSF, circular aperture, compensated gaussian fluxes.
1181
1182 Parameters
1183 ----------
1184 exposure : `lsst.afw.image.Exposure`
1185 Exposure to detect and measure stars on.
1186 background : `lsst.afw.math.BackgroundList`
1187 Background that was fit to the exposure during detection;
1188 modified in-place during subsequent detection.
1189 id_generator : `lsst.meas.base.IdGenerator`
1190 Object that generates source IDs and provides random seeds.
1191 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1192 Image to convert photometric-flattened image to
1193 background-flattened image.
1194
1195 Returns
1196 -------
1197 stars : `SourceCatalog`
1198 Sources that are very likely to be stars, with a limited set of
1199 measurements performed on them.
1200 """
1201 table = afwTable.SourceTable.make(self.initial_stars_schema.schema,
1202 id_generator.make_table_id_factory())
1203 # Re-estimate the background during this detection step, so that
1204 # measurement uses the most accurate background-subtraction.
1205 detections = self.star_detection.run(
1206 table=table,
1207 exposure=exposure,
1208 background=background,
1209 backgroundToPhotometricRatio=background_to_photometric_ratio,
1210 )
1211 sources = detections.sources
1212 self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources)
1213
1214 n_sky_sources = np.sum(sources["sky_source"])
1215 if (self.config.do_downsample_footprints
1216 and (len(sources) - n_sky_sources) > self.config.downsample_max_footprints):
1217 if exposure.info.id is None:
1218 self.log.warning("Exposure does not have a proper id; using 0 seed for downsample.")
1219 seed = 0
1220 else:
1221 seed = exposure.info.id & 0xFFFFFFFF
1222
1223 gen = np.random.RandomState(seed)
1224
1225 # Isolate the sky sources before downsampling.
1226 indices = np.arange(len(sources))[~sources["sky_source"]]
1227 indices = gen.choice(
1228 indices,
1229 size=self.config.downsample_max_footprints,
1230 replace=False,
1231 )
1232 skyIndices, = np.where(sources["sky_source"])
1233 indices = np.concatenate((indices, skyIndices))
1234
1235 self.log.info("Downsampling from %d to %d non-sky-source footprints.", len(sources), len(indices))
1236
1237 sel = np.zeros(len(sources), dtype=bool)
1238 sel[indices] = True
1239 sources = sources[sel]
1240
1241 # TODO investigation: Could this deblender throw away blends of non-PSF sources?
1242 self.star_deblend.run(exposure=exposure, sources=sources)
1243 # The deblender may not produce a contiguous catalog; ensure
1244 # contiguity for subsequent tasks.
1245 if not sources.isContiguous():
1246 sources = sources.copy(deep=True)
1247
1248 # Measure everything, and use those results to select only stars.
1249 self.star_measurement.run(sources, exposure)
1250 self.metadata["post_deblend_source_count"] = np.sum(~sources["sky_source"])
1251 self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"])
1252 self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"])
1253
1254 # Run the normalization calibration flux task to apply the
1255 # normalization correction to create normalized
1256 # calibration fluxes.
1257 self.star_normalized_calibration_flux.run(exposure=exposure, catalog=sources)
1258 self.star_apply_aperture_correction.run(sources, exposure.apCorrMap)
1259 self.star_catalog_calculation.run(sources)
1260 self.star_set_primary_flags.run(sources)
1261
1262 result = self.star_selector.run(sources)
1263 # The star selector may not produce a contiguous catalog.
1264 if not result.sourceCat.isContiguous():
1265 return result.sourceCat.copy(deep=True)
1266 else:
1267 return result.sourceCat
1268
1269 def _match_psf_stars(self, psf_stars, stars):
1270 """Match calibration stars to psf stars, to identify which were psf
1271 candidates, and which were used or reserved during psf measurement
1272 and the astrometric fit.
1273
1274 Parameters
1275 ----------
1276 psf_stars : `lsst.afw.table.SourceCatalog`
1277 PSF candidate stars that were sent to the psf determiner and
1278 used in the astrometric fit. Used to populate psf and astrometry
1279 related flag fields.
1280 stars : `lsst.afw.table.SourceCatalog`
1281 Stars that will be used for calibration; psf-related fields will
1282 be updated in-place.
1283
1284 Notes
1285 -----
1286 This code was adapted from CalibrateTask.copyIcSourceFields().
1287 """
1288 control = afwTable.MatchControl()
1289 # Return all matched objects, to separate blends.
1290 control.findOnlyClosest = False
1291 matches = afwTable.matchXy(psf_stars, stars, 3.0, control)
1292 deblend_key = stars.schema["deblend_nChild"].asKey()
1293 matches = [m for m in matches if m[1].get(deblend_key) == 0]
1294
1295 # Because we had to allow multiple matches to handle parents, we now
1296 # need to prune to the best (closest) matches.
1297 # Closest matches is a dict of psf_stars source ID to Match record
1298 # (psf_stars source, sourceCat source, distance in pixels).
1299 best = {}
1300 for match_psf, match_stars, d in matches:
1301 match = best.get(match_psf.getId())
1302 if match is None or d <= match[2]:
1303 best[match_psf.getId()] = (match_psf, match_stars, d)
1304 matches = list(best.values())
1305 # We'll use this to construct index arrays into each catalog.
1306 ids = np.array([(match_psf.getId(), match_stars.getId()) for match_psf, match_stars, d in matches]).T
1307
1308 if (n_matches := len(matches)) == 0:
1309 raise NoPsfStarsToStarsMatchError(n_psf_stars=len(psf_stars), n_stars=len(stars))
1310
1311 self.log.info("%d psf/astrometry stars out of %d matched %d calib stars",
1312 n_matches, len(psf_stars), len(stars))
1313 self.metadata["matched_psf_star_count"] = n_matches
1314
1315 # Check that no stars sources are listed twice; we already know
1316 # that each match has a unique psf_stars id, due to using as the key
1317 # in best above.
1318 n_unique = len(set(m[1].getId() for m in matches))
1319 if n_unique != n_matches:
1320 self.log.warning("%d psf_stars matched only %d stars", n_matches, n_unique)
1321
1322 # The indices of the IDs, so we can update the flag fields as arrays.
1323 idx_psf_stars = np.searchsorted(psf_stars["id"], ids[0])
1324 idx_stars = np.searchsorted(stars["id"], ids[1])
1325 for field in self.psf_fields:
1326 result = np.zeros(len(stars), dtype=bool)
1327 result[idx_stars] = psf_stars[field][idx_psf_stars]
1328 stars[field] = result
1329 stars['psf_id'][idx_stars] = psf_stars['id'][idx_psf_stars]
1330 stars['psf_max_value'][idx_stars] = psf_stars['psf_max_value'][idx_psf_stars]
1331
1332 def _fit_astrometry(self, exposure, stars):
1333 """Fit an astrometric model to the data and return the reference
1334 matches used in the fit, and the fitted WCS.
1335
1336 Parameters
1337 ----------
1338 exposure : `lsst.afw.image.Exposure`
1339 Exposure that is being fit, to get PSF and other metadata from.
1340 Modified to add the fitted skyWcs.
1341 stars : `SourceCatalog`
1342 Good stars selected for use in calibration, with RA/Dec coordinates
1343 computed from the pixel positions and fitted WCS.
1344
1345 Returns
1346 -------
1347 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1348 Reference/stars matches used in the fit.
1349 """
1350 result = self.astrometry.run(stars, exposure)
1351 return result.matches, result.matchMeta
1352
1353 def _fit_photometry(self, exposure, stars):
1354 """Fit a photometric model to the data and return the reference
1355 matches used in the fit, and the fitted PhotoCalib.
1356
1357 Parameters
1358 ----------
1359 exposure : `lsst.afw.image.Exposure`
1360 Exposure that is being fit, to get PSF and other metadata from.
1361 Has the fit `lsst.afw.image.PhotoCalib` attached, with pixel values
1362 unchanged.
1363 stars : `lsst.afw.table.SourceCatalog`
1364 Good stars selected for use in calibration.
1365 background : `lsst.afw.math.BackgroundList`
1366 Background that was fit to the exposure during detection of the
1367 above stars.
1368
1369 Returns
1370 -------
1371 calibrated_stars : `lsst.afw.table.SourceCatalog`
1372 Star catalog with flux/magnitude columns computed from the fitted
1373 photoCalib (instFlux columns are retained as well).
1374 matches : `list` [`lsst.afw.table.ReferenceMatch`]
1375 Reference/stars matches used in the fit.
1376 matchMeta : `lsst.daf.base.PropertyList`
1377 Metadata needed to unpersist matches, as returned by the matcher.
1378 photo_calib : `lsst.afw.image.PhotoCalib`
1379 Photometric calibration that was fit to the star catalog.
1380 """
1381 result = self.photometry.run(exposure, stars)
1382 calibrated_stars = result.photoCalib.calibrateCatalog(stars)
1383 exposure.setPhotoCalib(result.photoCalib)
1384 return calibrated_stars, result.matches, result.matchMeta, result.photoCalib
1385
1386 def _apply_photometry(self, exposure, background, background_to_photometric_ratio=None):
1387 """Apply the photometric model attached to the exposure to the
1388 exposure's pixels and an associated background model.
1389
1390 Parameters
1391 ----------
1392 exposure : `lsst.afw.image.Exposure`
1393 Exposure with the target `lsst.afw.image.PhotoCalib` attached.
1394 On return, pixel values will be calibrated and an identity
1395 photometric transform will be attached.
1396 background : `lsst.afw.math.BackgroundList`
1397 Background model to convert to nanojansky units in place.
1398 background_to_photometric_ratio : `lsst.afw.image.Image`, optional
1399 Image to convert photometric-flattened image to
1400 background-flattened image.
1401 """
1402 photo_calib = exposure.getPhotoCalib()
1403 exposure.maskedImage = photo_calib.calibrateImage(exposure.maskedImage)
1404 identity = afwImage.PhotoCalib(1.0,
1405 photo_calib.getCalibrationErr(),
1406 bbox=exposure.getBBox())
1407 exposure.setPhotoCalib(identity)
1408 exposure.metadata["BUNIT"] = "nJy"
1409
1410 assert photo_calib._isConstant, \
1411 "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that."
1412
1413 for bg in background:
1414 # The statsImage is a view, but we can't assign to a function call in python.
1415 binned_image = bg[0].getStatsImage()
1416 binned_image *= photo_calib.getCalibrationMean()
1417
1418 def _summarize(self, exposure, stars, background):
1419 """Compute summary statistics on the exposure and update in-place the
1420 calibrations attached to it.
1421
1422 Parameters
1423 ----------
1424 exposure : `lsst.afw.image.Exposure`
1425 Exposure that was calibrated, to get PSF and other metadata from.
1426 Should be in instrumental units with the photometric calibration
1427 attached.
1428 Modified to contain the computed summary statistics.
1429 stars : `SourceCatalog`
1430 Good stars selected used in calibration.
1431 background : `lsst.afw.math.BackgroundList`
1432 Background that was fit to the exposure during detection of the
1433 above stars. Should be in instrumental units.
1434 """
1435 summary = self.compute_summary_stats.run(exposure, stars, background)
1436 exposure.info.setSummaryStats(summary)
1437
1438 def _recordMaskedPixelFractions(self, exposure):
1439 """Record the fraction of all the pixels in an exposure
1440 that are masked with a given flag. Each fraction is
1441 recorded in the task metadata. One record per flag type.
1442
1443 Parameters
1444 ----------
1445 exposure : `lsst.afw.image.ExposureF`
1446 The target exposure to calculate masked pixel fractions for.
1447 """
1448
1449 mask = exposure.mask
1450 maskPlanes = list(mask.getMaskPlaneDict().keys())
1451 for maskPlane in maskPlanes:
1452 self.metadata[f"{maskPlane.lower()}_mask_fraction"] = (
1453 evaluateMaskFraction(mask, maskPlane)
1454 )
1455
1456 def _update_wcs_with_camera_model(self, exposure, cameraModel):
1457 """Replace the existing WCS with one generated using the pointing
1458 and the input camera model.
1459
1460 Parameters
1461 ----------
1462 exposure : `lsst.afw.image.ExposureF`
1463 The exposure whose WCS will be updated.
1464 cameraModel : `lsst.afw.cameraGeom.Camera`
1465 Camera to be used when constructing updated WCS.
1466 """
1467 detector = cameraModel[exposure.detector.getId()]
1468 newWcs = createInitialSkyWcs(exposure.visitInfo, detector)
1469 exposure.setWcs(newWcs)
The photometric calibration of an exposure.
Definition PhotoCalib.h:114
Pass parameters to algorithms that match list of sources.
Definition Match.h:45
An integer coordinate rectangle.
Definition Box.h:55
__init__(self, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size)
runQuantum(self, butlerQC, inputRefs, outputRefs)
_update_wcs_with_camera_model(self, exposure, cameraModel)
_find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None)
_apply_illumination_correction(self, exposure, background_flat, illumination_correction)
_measure_aperture_correction(self, exposure, bright_sources)
_apply_photometry(self, exposure, background, background_to_photometric_ratio=None)
run(self, *, exposures, id_generator=None, result=None, background_flat=None, illumination_correction=None, camera_model=None)
__init__(self, initial_stars_schema=None, **kwargs)
_compute_psf(self, exposure, id_generator, background_to_photometric_ratio=None)
_summarize(self, exposure, stars, background)
SourceMatchVector matchXy(SourceCatalog const &cat1, SourceCatalog const &cat2, double radius, MatchControl const &mc=MatchControl())
Compute all tuples (s1,s2,d) where s1 belings to cat1, s2 belongs to cat2 and d, the distance between...
Definition Match.cc:305
void updateSourceCoords(geom::SkyWcs const &wcs, SourceCollection &sourceList, bool include_covariance=true)
Update sky coordinates in a collection of source objects.
Definition wcsUtils.cc:125