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