Coverage for python/lsst/pipe/tasks/calibrateImage.py: 12%

1# This file is part of pipe_tasks.

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.

9# This program is free software: you can redistribute it and/or modify

10# it under the terms of the GNU General Public License as published by

11# the Free Software Foundation, either version 3 of the License, or

12# (at your option) any later version.

13#

14# This program is distributed in the hope that it will be useful,

15# but WITHOUT ANY WARRANTY; without even the implied warranty of

16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the

17# GNU General Public License for more details.

18#

19# You should have received a copy of the GNU General Public License

20# along with this program. If not, see <https://www.gnu.org/licenses/>.

22__all__ = ["CalibrateImageTask", "CalibrateImageConfig", "NoPsfStarsToStarsMatchError",

23 "AllCentroidsFlaggedError"]

25import math

26import numpy as np

27import requests

28import os

30from lsst.geom import SpherePoint, degrees

31from lsst.afw.geom import SpanSet

32import lsst.afw.table as afwTable

33import lsst.afw.image as afwImage

34import lsst.afw.math as afwMath

35from lsst.ip.diffim.utils import evaluateMaskFraction, populate_sattle_visit_cache

36import lsst.meas.algorithms

37import lsst.meas.algorithms.installGaussianPsf

38import lsst.meas.algorithms.measureApCorr

39import lsst.meas.algorithms.setPrimaryFlags

40from lsst.meas.algorithms.adaptive_thresholds import AdaptiveThresholdDetectionTask

41import lsst.meas.base

42import lsst.meas.astrom

43import lsst.meas.deblender

44import lsst.meas.extensions.shapeHSM

45from lsst.obs.base.utils import createInitialSkyWcsFromBoresight

46import lsst.pex.config as pexConfig

47import lsst.pipe.base as pipeBase

48from lsst.pipe.base import connectionTypes

49from lsst.utils.timer import timeMethod

51from . import measurePsf, repair, photoCal, computeExposureSummaryStats, snapCombine, diffractionSpikeMask

54class AllCentroidsFlaggedError(pipeBase.AlgorithmError):

55 """Raised when there are no valid centroids during psf fitting.

56 """

57 def __init__(self, n_sources, psf_shape_ixx, psf_shape_iyy, psf_shape_ixy, psf_size):

58 msg = (f"All source centroids (out of {n_sources}) flagged during PSF fitting. "

59 "Original image PSF is likely unusable; best-fit PSF shape parameters: "

60 f"Ixx={psf_shape_ixx}, Iyy={psf_shape_iyy}, Ixy={psf_shape_ixy}, size={psf_size}"

61 )

62 super().__init__(msg)

63 self.n_sources = n_sources

64 self.psf_shape_ixx = psf_shape_ixx

65 self.psf_shape_iyy = psf_shape_iyy

66 self.psf_shape_ixy = psf_shape_ixy

67 self.psf_size = psf_size

69 @property

70 def metadata(self):

71 return {"n_sources": self.n_sources,

72 "psf_shape_ixx": self.psf_shape_ixx,

73 "psf_shape_iyy": self.psf_shape_iyy,

74 "psf_shape_ixy": self.psf_shape_ixy,

75 "psf_size": self.psf_size

76 }

79class NoPsfStarsToStarsMatchError(pipeBase.AlgorithmError):

80 """Raised when there are no matches between the psf_stars and stars

81 catalogs.

82 """

83 def __init__(self, *, n_psf_stars, n_stars):

84 msg = (f"No psf stars out of {n_psf_stars} matched {n_stars} calib stars."

85 " Downstream processes probably won't have useful {calib_type} stars in this case."

86 " Is `star_selector` too strict or is this a bad image?")

87 super().__init__(msg)

88 self.n_psf_stars = n_psf_stars

89 self.n_stars = n_stars

91 @property

92 def metadata(self):

93 return {"n_psf_stars": self.n_psf_stars,

94 "n_stars": self.n_stars

95 }

98class CalibrateImageConnections(pipeBase.PipelineTaskConnections,

99 dimensions=("instrument", "visit", "detector")):

100

101 astrometry_ref_cat = connectionTypes.PrerequisiteInput(

102 doc="Reference catalog to use for astrometric calibration.",

103 name="gaia_dr3_20230707",

104 storageClass="SimpleCatalog",

105 dimensions=("skypix",),

106 deferLoad=True,

107 multiple=True,

108 )

109 photometry_ref_cat = connectionTypes.PrerequisiteInput(

110 doc="Reference catalog to use for photometric calibration.",

111 name="ps1_pv3_3pi_20170110",

112 storageClass="SimpleCatalog",

113 dimensions=("skypix",),

114 deferLoad=True,

115 multiple=True

116 )

117 camera_model = connectionTypes.PrerequisiteInput(

118 doc="Camera distortion model to use for astrometric calibration.",

119 name="astrometry_camera",

120 storageClass="Camera",

121 dimensions=("instrument", "physical_filter"),

122 isCalibration=True,

123 minimum=0, # Can fall back on WCS already attached to exposure.

124 )

125 exposures = connectionTypes.Input(

126 doc="Exposure (or two snaps) to be calibrated, and detected and measured on.",

127 name="postISRCCD",

128 storageClass="Exposure",

129 multiple=True, # to handle 1 exposure or 2 snaps

130 dimensions=["instrument", "exposure", "detector"],

131 )

132

133 background_flat = connectionTypes.PrerequisiteInput(

134 name="flat",

135 doc="Flat calibration frame used for background correction.",

136 storageClass="ExposureF",

137 dimensions=["instrument", "detector", "physical_filter"],

138 isCalibration=True,

139 )

140 illumination_correction = connectionTypes.PrerequisiteInput(

141 name="illuminationCorrection",

142 doc="Illumination correction frame.",

143 storageClass="Exposure",

144 dimensions=["instrument", "detector", "physical_filter"],

145 isCalibration=True,

146 )

147

148 # outputs

149 initial_stars_schema = connectionTypes.InitOutput(

150 doc="Schema of the output initial stars catalog.",

151 name="initial_stars_schema",

152 storageClass="SourceCatalog",

153 )

154

155 # TODO DM-38732: We want some kind of flag on Exposures/Catalogs to make

156 # it obvious which components had failed to be computed/persisted.

157 exposure = connectionTypes.Output(

158 doc="Photometrically calibrated, background-subtracted exposure with fitted calibrations and "

159 "summary statistics. To recover the original exposure, first add the background "

160 "(`initial_pvi_background`), and then uncalibrate (divide by `initial_photoCalib_detector`).",

161 name="initial_pvi",

162 storageClass="ExposureF",

163 dimensions=("instrument", "visit", "detector"),

164 )

165 stars = connectionTypes.Output(

166 doc="Catalog of unresolved sources detected on the calibrated exposure.",

167 name="initial_stars_detector",

168 storageClass="ArrowAstropy",

169 dimensions=["instrument", "visit", "detector"],

170 )

171 stars_footprints = connectionTypes.Output(

172 doc="Catalog of unresolved sources detected on the calibrated exposure; "

173 "includes source footprints.",

174 name="initial_stars_footprints_detector",

175 storageClass="SourceCatalog",

176 dimensions=["instrument", "visit", "detector"],

177 )

178 applied_photo_calib = connectionTypes.Output(

179 doc=(

180 "Photometric calibration that was applied to exposure's pixels. "

181 "This connection is disabled when do_calibrate_pixels=False."

182 ),

183 name="initial_photoCalib_detector",

184 storageClass="PhotoCalib",

185 dimensions=("instrument", "visit", "detector"),

186 )

187 background = connectionTypes.Output(

188 doc="Background models estimated during calibration task; calibrated to be in nJy units.",

189 name="initial_pvi_background",

190 storageClass="Background",

191 dimensions=("instrument", "visit", "detector"),

192 )

193 background_to_photometric_ratio = connectionTypes.Output(

194 doc="Ratio of a background-flattened image to a photometric-flattened image. Only persisted "

195 "if do_illumination_correction is True.",

196 name="background_to_photometric_ratio",

197 storageClass="Image",

198 dimensions=("instrument", "visit", "detector"),

199 )

200

201 # Optional outputs

202 mask = connectionTypes.Output(

203 doc="The mask plane of the calibrated exposure, written as a separate "

204 "output so that it can be passed to downstream tasks even if the "

205 "exposure is removed to save storage.",

206 name="preliminary_visit_mask",

207 storageClass="Mask",

208 dimensions=("instrument", "visit", "detector"),

209 )

210 psf_stars_footprints = connectionTypes.Output(

211 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination; "

212 "includes source footprints.",

213 name="initial_psf_stars_footprints_detector",

214 storageClass="SourceCatalog",

215 dimensions=["instrument", "visit", "detector"],

216 )

217 psf_stars = connectionTypes.Output(

218 doc="Catalog of bright unresolved sources detected on the exposure used for PSF determination.",

219 name="initial_psf_stars_detector",

220 storageClass="ArrowAstropy",

221 dimensions=["instrument", "visit", "detector"],

222 )

223 astrometry_matches = connectionTypes.Output(

224 doc="Source to reference catalog matches from the astrometry solver.",

225 name="initial_astrometry_match_detector",

226 storageClass="Catalog",

227 dimensions=("instrument", "visit", "detector"),

228 )

229 photometry_matches = connectionTypes.Output(

230 doc="Source to reference catalog matches from the photometry solver.",

231 name="initial_photometry_match_detector",

232 storageClass="Catalog",

233 dimensions=("instrument", "visit", "detector"),

234 )

235

236 def __init__(self, *, config=None):

237 super().__init__(config=config)

238 if "psf_stars" not in config.optional_outputs:

239 del self.psf_stars

240 if "psf_stars_footprints" not in config.optional_outputs:

241 del self.psf_stars_footprints

242 if "astrometry_matches" not in config.optional_outputs:

243 del self.astrometry_matches

244 if "photometry_matches" not in config.optional_outputs:

245 del self.photometry_matches

246 if "mask" not in config.optional_outputs:

247 del self.mask

248 if not config.do_calibrate_pixels:

249 del self.applied_photo_calib

250 if not config.do_illumination_correction:

251 del self.background_flat

252 del self.illumination_correction

253 del self.background_to_photometric_ratio

254 if not config.useButlerCamera:

255 del self.camera_model

256

257

258class CalibrateImageConfig(pipeBase.PipelineTaskConfig, pipelineConnections=CalibrateImageConnections):

259 optional_outputs = pexConfig.ListField(

260 doc="Which optional outputs to save (as their connection name)?",

261 dtype=str,

262 # TODO: note somewhere to disable this for benchmarking, but should

263 # we always have it on for production runs?

264 default=["psf_stars", "psf_stars_footprints", "astrometry_matches", "photometry_matches", "mask"],

265 optional=False

266 )

267

268 # To generate catalog ids consistently across subtasks.

269 id_generator = lsst.meas.base.DetectorVisitIdGeneratorConfig.make_field()

270

271 snap_combine = pexConfig.ConfigurableField(

272 target=snapCombine.SnapCombineTask,

273 doc="Task to combine two snaps to make one exposure.",

274 )

275

276 # subtasks used during psf characterization

277 install_simple_psf = pexConfig.ConfigurableField(

278 target=lsst.meas.algorithms.installGaussianPsf.InstallGaussianPsfTask,

279 doc="Task to install a simple PSF model into the input exposure to use "

280 "when detecting bright sources for PSF estimation.",

281 )

282 psf_repair = pexConfig.ConfigurableField(

283 target=repair.RepairTask,

284 doc="Task to repair cosmic rays on the exposure before PSF determination.",

285 )

286 psf_subtract_background = pexConfig.ConfigurableField(

287 target=lsst.meas.algorithms.SubtractBackgroundTask,

288 doc="Task to perform initial background subtraction, before first detection pass.",

289 )

290 psf_detection = pexConfig.ConfigurableField(

291 target=lsst.meas.algorithms.SourceDetectionTask,

292 doc="Task to detect sources for PSF determination."

293 )

294 do_adaptive_threshold_detection = pexConfig.Field(

295 dtype=bool,

296 default=True,

297 doc="Implement the adaptive detection thresholding approach?",

298 )

299 psf_adaptive_threshold_detection = pexConfig.ConfigurableField(

300 target=AdaptiveThresholdDetectionTask,

301 doc="Task to adaptively detect sources for PSF determination.",

302 )

303 psf_source_measurement = pexConfig.ConfigurableField(

304 target=lsst.meas.base.SingleFrameMeasurementTask,

305 doc="Task to measure sources to be used for psf estimation."

306 )

307 psf_measure_psf = pexConfig.ConfigurableField(

308 target=measurePsf.MeasurePsfTask,

309 doc="Task to measure the psf on bright sources."

310 )

311 psf_normalized_calibration_flux = pexConfig.ConfigurableField(

312 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,

313 doc="Task to normalize the calibration flux (e.g. compensated tophats) "

314 "for the bright stars used for psf estimation.",

315 )

316

317 # TODO DM-39203: we can remove aperture correction from this task once we

318 # are using the shape-based star/galaxy code.

319 measure_aperture_correction = pexConfig.ConfigurableField(

320 target=lsst.meas.algorithms.measureApCorr.MeasureApCorrTask,

321 doc="Task to compute the aperture correction from the bright stars."

322 )

323

324 # subtasks used during star measurement

325 star_background = pexConfig.ConfigurableField(

326 target=lsst.meas.algorithms.SubtractBackgroundTask,

327 doc="Task to perform final background subtraction, just before photoCal.",

328 )

329 star_background_min_footprints = pexConfig.Field(

330 dtype=int,

331 default=3,

332 doc="Minimum number of footprints in the detection mask for star_background measurement. "

333 "This number will get adjusted to the fraction config.star_background_peak_fraction "

334 "of the detected peaks if that number is larger. If the number of footprints is less "

335 "than the minimum, the detection threshold is iteratively increased until the "

336 "threshold is met.",

337 )

338 star_background_peak_fraction = pexConfig.Field(

339 dtype=float,

340 default=0.01,

341 doc="The minimum number of footprints in the detection mask for star_background measurement. "

342 "gets set to the maximum of this fraction of the detected peaks and the value set in "

343 "config.star_background_min_footprints. If the number of footprints is less than the "

344 "current minimum set, the detection threshold is iteratively increased until the "

345 "threshold is met.",

346 )

347 star_detection = pexConfig.ConfigurableField(

348 target=lsst.meas.algorithms.SourceDetectionTask,

349 doc="Task to detect stars to return in the output catalog."

350 )

351 star_sky_sources = pexConfig.ConfigurableField(

352 target=lsst.meas.algorithms.SkyObjectsTask,

353 doc="Task to generate sky sources ('empty' regions where there are no detections).",

354 )

355 do_downsample_footprints = pexConfig.Field(

356 dtype=bool,

357 default=False,

358 doc="Downsample footprints prior to deblending to optimize speed?",

359 )

360 downsample_max_footprints = pexConfig.Field(

361 dtype=int,

362 default=1000,

363 doc="Maximum number of non-sky-source footprints to use if do_downsample_footprints is True,",

364 )

365 star_deblend = pexConfig.ConfigurableField(

366 target=lsst.meas.deblender.SourceDeblendTask,

367 doc="Split blended sources into their components."

368 )

369 star_measurement = pexConfig.ConfigurableField(

370 target=lsst.meas.base.SingleFrameMeasurementTask,

371 doc="Task to measure stars to return in the output catalog."

372 )

373 star_normalized_calibration_flux = pexConfig.ConfigurableField(

374 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,

375 doc="Task to apply the normalization for calibration fluxes (e.g. compensated tophats) "

376 "for the final output star catalog.",

377 )

378 star_apply_aperture_correction = pexConfig.ConfigurableField(

379 target=lsst.meas.base.ApplyApCorrTask,

380 doc="Task to apply aperture corrections to the selected stars."

381 )

382 star_catalog_calculation = pexConfig.ConfigurableField(

383 target=lsst.meas.base.CatalogCalculationTask,

384 doc="Task to compute extendedness values on the star catalog, "

385 "for the star selector to remove extended sources."

386 )

387 star_set_primary_flags = pexConfig.ConfigurableField(

388 target=lsst.meas.algorithms.setPrimaryFlags.SetPrimaryFlagsTask,

389 doc="Task to add isPrimary to the catalog."

390 )

391 star_selector = lsst.meas.algorithms.sourceSelectorRegistry.makeField(

392 default="science",

393 doc="Task to select reliable stars to use for calibration."

394 )

395

396 # final calibrations and statistics

397 astrometry = pexConfig.ConfigurableField(

398 target=lsst.meas.astrom.AstrometryTask,

399 doc="Task to perform astrometric calibration to fit a WCS.",

400 )

401 astrometry_ref_loader = pexConfig.ConfigField(

402 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,

403 doc="Configuration of reference object loader for astrometric fit.",

404 )

405 photometry = pexConfig.ConfigurableField(

406 target=photoCal.PhotoCalTask,

407 doc="Task to perform photometric calibration to fit a PhotoCalib.",

408 )

409 photometry_ref_loader = pexConfig.ConfigField(

410 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,

411 doc="Configuration of reference object loader for photometric fit.",

412 )

413 doMaskDiffractionSpikes = pexConfig.Field(

414 dtype=bool,

415 default=False,

416 doc="If True, load the photometric reference catalog again but select "

417 "only bright stars. Use the bright star catalog to set the SPIKE "

418 "mask for regions likely contaminated by diffraction spikes.",

419 )

420 diffractionSpikeMask = pexConfig.ConfigurableField(

421 target=diffractionSpikeMask.DiffractionSpikeMaskTask,

422 doc="Task to identify and mask the diffraction spikes of bright stars.",

423 )

424

425 compute_summary_stats = pexConfig.ConfigurableField(

426 target=computeExposureSummaryStats.ComputeExposureSummaryStatsTask,

427 doc="Task to to compute summary statistics on the calibrated exposure."

428 )

429

430 do_illumination_correction = pexConfig.Field(

431 dtype=bool,

432 default=False,

433 doc="If True, apply the illumination correction. This assumes that the "

434 "input image has already been flat-fielded such that it is suitable "

435 "for background subtraction.",

436 )

437

438 do_calibrate_pixels = pexConfig.Field(

439 dtype=bool,

440 default=True,

441 doc=(

442 "If True, apply the photometric calibration to the image pixels "

443 "and background model, and attach an identity PhotoCalib to "

444 "the output image to reflect this. If False`, leave the image "

445 "and background uncalibrated and attach the PhotoCalib that maps "

446 "them to physical units."

447 )

448 )

449

450 do_include_astrometric_errors = pexConfig.Field(

451 dtype=bool,

452 default=True,

453 doc="If True, include astrometric errors in the output catalog.",

454 )

455

456 background_stats_ignored_pixel_masks = pexConfig.ListField(

457 dtype=str,

458 default=["SAT", "SUSPECT", "SPIKE"],

459 doc="Pixel mask flags to ignore when calculating post-sky-subtraction "

460 "background statistics. These are added to those ignored by the "

461 "meas.algorithms.SubtractBackgroundConfig algorithm."

462 )

463

464 run_sattle = pexConfig.Field(

465 dtype=bool,

466 default=False,

467 doc="If True, the sattle service will populate a cache for later use "

468 "in ip_diffim.detectAndMeasure alert verification."

469 )

470

471 sattle_historical = pexConfig.Field(

472 dtype=bool,

473 default=False,

474 doc="If re-running a pipeline that requires sattle, this should be set "

475 "to True. This will populate sattle's cache with the historic data "

476 "closest in time to the exposure.",

477 )

478

479 useButlerExposureRegion = pexConfig.Field(

480 dtype=bool,

481 default=False,

482 doc="If True, use the exposure region stored in the butler as the "

483 "region for reference catalog trimming. If False, the reference loader "

484 "trimming will be set by the loader (typically from a padded version of "

485 "the detector's bounding box + current WCS)."

486 )

487

488 useButlerCamera = pexConfig.Field(

489 dtype=bool,

490 default=False,

491 doc="If True, use a camera distortion model generated elsewhere in "

492 "the pipeline combined with the telescope boresight as a starting "

493 "point for fitting the WCS, instead of using the WCS attached to "

494 "the exposure, which is generated from the boresight and the "

495 "camera model from the obs_* package."

496 )

497

498 background_stats_flux_column = pexConfig.Field(

499 dtype=str,

500 default="base_CircularApertureFlux_12_0_flux",

501 doc="Column used to generate post-subtracted background stats."

502 )

503

504 def setDefaults(self):

505 super().setDefaults()

506

507 # Use a very broad PSF here, to thoroughly reject CRs.

508 # TODO investigation: a large initial psf guess may make stars look

509 # like CRs for very good seeing images.

510 self.install_simple_psf.fwhm = 4

511

512 # S/N>=50 sources for PSF determination, but detection to S/N=10.

513 # The thresholdValue sets the minimum flux in a pixel to be included

514 # in the footprint, while peaks are only detected when they are above

515 # thresholdValue * includeThresholdMultiplier. The low thresholdValue

516 # ensures that the footprints are large enough for the noise replacer

517 # to mask out faint undetected neighbors that are not to be measured.

518 self.psf_detection.thresholdValue = 10.0

519 self.psf_detection.includeThresholdMultiplier = 5.0

520 # TODO investigation: Probably want False here, but that may require

521 # tweaking the background spatial scale, to make it small enough to

522 # prevent extra peaks in the wings of bright objects.

523 self.psf_detection.doTempLocalBackground = False

524 self.psf_detection.reEstimateBackground = False

525

526 # Minimal measurement plugins for PSF determination.

527 # TODO DM-39203: We can drop GaussianFlux and PsfFlux, if we use

528 # shapeHSM/moments for star/galaxy separation.

529 # TODO DM-39203: we can remove aperture correction from this task

530 # once we are using the shape-based star/galaxy code.

531 self.psf_source_measurement.plugins = ["base_PixelFlags",

532 "base_SdssCentroid",

533 "ext_shapeHSM_HsmSourceMoments",

534 "base_CircularApertureFlux",

535 "base_GaussianFlux",

536 "base_PsfFlux",

537 "base_CompensatedTophatFlux",

538 ]

539 self.psf_source_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"

540 # Only measure apertures we need for PSF measurement.

541 self.psf_source_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]

542 self.psf_source_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]

543 # TODO DM-40843: Remove this line once this is the psfex default.

544 self.psf_measure_psf.psfDeterminer["psfex"].photometricFluxField = \

545 "base_CircularApertureFlux_12_0_instFlux"

546

547 # No extendedness information available: we need the aperture

548 # corrections to determine that.

549 self.measure_aperture_correction.sourceSelector["science"].doUnresolved = False

550 self.measure_aperture_correction.sourceSelector["science"].flags.good = ["calib_psf_used"]

551 self.measure_aperture_correction.sourceSelector["science"].flags.bad = []

552

553 # Detection for good S/N for astrometry/photometry and other

554 # downstream tasks; detection mask to S/N>=5, but S/N>=10 peaks.

555 self.star_detection.reEstimateBackground = False

556 self.star_detection.thresholdValue = 5.0

557 self.star_detection.includeThresholdMultiplier = 2.0

558 self.star_measurement.plugins = ["base_PixelFlags",

559 "base_SdssCentroid",

560 "ext_shapeHSM_HsmSourceMoments",

561 "ext_shapeHSM_HsmPsfMoments",

562 "base_GaussianFlux",

563 "base_PsfFlux",

564 "base_CircularApertureFlux",

565 "base_ClassificationSizeExtendedness",

566 "base_CompensatedTophatFlux",

567 ]

568 self.star_measurement.slots.psfShape = "ext_shapeHSM_HsmPsfMoments"

569 self.star_measurement.slots.shape = "ext_shapeHSM_HsmSourceMoments"

570 # Only measure the apertures we need for star selection.

571 self.star_measurement.plugins["base_CircularApertureFlux"].radii = [12.0]

572 self.star_measurement.plugins["base_CompensatedTophatFlux"].apertures = [12]

573

574 # We measure and apply the normalization aperture correction with

575 # the psf_normalized_calibration_flux task, and we only apply the

576 # normalization aperture correction for the full list of stars.

577 self.star_normalized_calibration_flux.do_measure_ap_corr = False

578

579 # Select stars with reliable measurements and no bad flags.

580 self.star_selector["science"].doFlags = True

581 self.star_selector["science"].doUnresolved = True

582 self.star_selector["science"].doSignalToNoise = True

583 self.star_selector["science"].signalToNoise.minimum = 10.0

584 # Keep sky sources in the output catalog, even though they aren't

585 # wanted for calibration.

586 self.star_selector["science"].doSkySources = True

587 # Set the flux and error fields

588 self.star_selector["science"].signalToNoise.fluxField = "slot_CalibFlux_instFlux"

589 self.star_selector["science"].signalToNoise.errField = "slot_CalibFlux_instFluxErr"

590

591 # Use the affine WCS fitter (assumes we have a good camera geometry).

592 self.astrometry.wcsFitter.retarget(lsst.meas.astrom.FitAffineWcsTask)

593 # phot_g_mean is the primary Gaia band for all input bands.

594 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"

595

596 # Only reject sky sources; we already selected good stars.

597 self.astrometry.sourceSelector["science"].doFlags = True

598 self.astrometry.sourceSelector["science"].flags.good = [] # ["calib_psf_candidate"]

599 # self.astrometry.sourceSelector["science"].flags.bad = []

600 self.astrometry.sourceSelector["science"].doUnresolved = False

601 self.astrometry.sourceSelector["science"].doIsolated = False

602 self.astrometry.sourceSelector["science"].doRequirePrimary = False

603 self.photometry.match.sourceSelection.doFlags = True

604 self.photometry.match.sourceSelection.flags.bad = ["sky_source"]

605 # Unset the (otherwise reasonable, but we've already made the

606 # selections we want above) selection settings in PhotoCalTask.

607 self.photometry.match.sourceSelection.doRequirePrimary = False

608 self.photometry.match.sourceSelection.doUnresolved = False

609

610 def validate(self):

611 super().validate()

612

613 # Ensure that the normalization calibration flux tasks

614 # are configured correctly.

615 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:

616 msg = ("psf_normalized_calibration_flux task must be configured with do_measure_ap_corr=True "

617 "or else the normalization and calibration flux will not be properly measured.")

618 raise pexConfig.FieldValidationError(

619 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,

620 )

621 if self.star_normalized_calibration_flux.do_measure_ap_corr:

622 msg = ("star_normalized_calibration_flux task must be configured with do_measure_ap_corr=False "

623 "to apply the previously measured normalization to the full catalog of calibration "

624 "fluxes.")

625 raise pexConfig.FieldValidationError(

626 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,

627 )

628

629 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,

630 # because they'd be running too early to pick up the fitted PhotoCalib

631 # and WCS.

632 if "base_LocalWcs" in self.psf_source_measurement.plugins.names:

633 raise pexConfig.FieldValidationError(

634 CalibrateImageConfig.psf_source_measurement,

635 self,

636 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."

637 )

638 if "base_LocalWcs" in self.star_measurement.plugins.names:

639 raise pexConfig.FieldValidationError(

640 CalibrateImageConfig.star_measurement,

641 self,

642 "base_LocalWcs cannot run CalibrateImageTask, as it would be run before the astrometry fit."

643 )

644 if "base_LocalPhotoCalib" in self.psf_source_measurement.plugins.names:

645 raise pexConfig.FieldValidationError(

646 CalibrateImageConfig.psf_source_measurement,

647 self,

648 "base_LocalPhotoCalib cannot run CalibrateImageTask, "

649 "as it would be run before the photometry fit."

650 )

651 if "base_LocalPhotoCalib" in self.star_measurement.plugins.names:

652 raise pexConfig.FieldValidationError(

653 CalibrateImageConfig.star_measurement,

654 self,

655 "base_LocalPhotoCalib cannot run CalibrateImageTask, "

656 "as it would be run before the photometry fit."

657 )

658

659 # Check for illumination correction and background consistency.

660 if self.do_illumination_correction:

661 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:

662 raise pexConfig.FieldValidationError(

663 CalibrateImageConfig.psf_subtract_background,

664 self,

665 "CalibrateImageTask.psf_subtract_background must be configured with "

666 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

667 )

668 if self.psf_detection.reEstimateBackground:

669 if not self.psf_detection.doApplyFlatBackgroundRatio:

670 raise pexConfig.FieldValidationError(

671 CalibrateImageConfig.psf_detection,

672 self,

673 "CalibrateImageTask.psf_detection background must be configured with "

674 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

675 )

676 if self.star_detection.reEstimateBackground:

677 if not self.star_detection.doApplyFlatBackgroundRatio:

678 raise pexConfig.FieldValidationError(

679 CalibrateImageConfig.star_detection,

680 self,

681 "CalibrateImageTask.star_detection background must be configured with "

682 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

683 )

684 if not self.star_background.doApplyFlatBackgroundRatio:

685 raise pexConfig.FieldValidationError(

686 CalibrateImageConfig.star_background,

687 self,

688 "CalibrateImageTask.star_background background must be configured with "

689 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

690 )

691

692 if self.run_sattle:

693 if not os.getenv("SATTLE_URI_BASE"):

694 raise pexConfig.FieldValidationError(CalibrateImageConfig.run_sattle, self,

695 "Sattle requested but URI environment variable not set.")

696

697 if not self.do_adaptive_threshold_detection:

698 if not self.psf_detection.reEstimateBackground:

699 raise pexConfig.FieldValidationError(

700 CalibrateImageConfig.psf_detection,

701 self,

702 "If not using the adaptive threshold detection implementation (i.e. "

703 "config.do_adaptive_threshold_detection = False), CalibrateImageTask.psf_detection "

704 "background must be configured with "

705 "reEstimateBackground = True to maintain original behavior."

706 )

707 if not self.star_detection.reEstimateBackground:

708 raise pexConfig.FieldValidationError(

709 CalibrateImageConfig.psf_detection,

710 self,

711 "If not using the adaptive threshold detection implementation "

712 "(i.e. config.do_adaptive_threshold_detection = False), "

713 "CalibrateImageTask.star_detection background must be configured "

714 "with reEstimateBackground = True to maintain original behavior."

715 )

716

717

718class CalibrateImageTask(pipeBase.PipelineTask):

719 """Compute the PSF, aperture corrections, astrometric and photometric

720 calibrations, and summary statistics for a single science exposure, and

721 produce a catalog of brighter stars that were used to calibrate it.

722

723 Parameters

724 ----------

725 initial_stars_schema : `lsst.afw.table.Schema`

726 Schema of the initial_stars output catalog.

727 """

728 _DefaultName = "calibrateImage"

729 ConfigClass = CalibrateImageConfig

730

731 def __init__(self, initial_stars_schema=None, **kwargs):

732 super().__init__(**kwargs)

733 self.makeSubtask("snap_combine")

734

735 # PSF determination subtasks

736 self.makeSubtask("install_simple_psf")

737 self.makeSubtask("psf_repair")

738 self.makeSubtask("psf_subtract_background")

739 self.psf_schema = afwTable.SourceTable.makeMinimalSchema()

740 self.psf_schema.addField(

741 'psf_max_value',

742 type=np.float32,

743 doc="PSF max value.",

744 )

745 afwTable.CoordKey.addErrorFields(self.psf_schema)

746 self.makeSubtask("psf_detection", schema=self.psf_schema)

747 self.makeSubtask("psf_source_measurement", schema=self.psf_schema)

748 self.makeSubtask("psf_adaptive_threshold_detection")

749 self.makeSubtask("psf_measure_psf", schema=self.psf_schema)

750 self.makeSubtask("psf_normalized_calibration_flux", schema=self.psf_schema)

751

752 self.makeSubtask("measure_aperture_correction", schema=self.psf_schema)

753 self.makeSubtask("astrometry", schema=self.psf_schema)

754

755 # star measurement subtasks

756 if initial_stars_schema is None:

757 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()

758

759 # These fields let us track which sources were used for psf modeling,

760 # astrometric fitting, and aperture correction calculations.

761 self.psf_fields = ("calib_psf_candidate", "calib_psf_used", "calib_psf_reserved",

762 "calib_astrometry_used",

763 # TODO DM-39203:

764 # these can be removed once apcorr is gone.

765 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",

766 "apcorr_base_PsfFlux_used",)

767 for field in self.psf_fields:

768 item = self.psf_schema.find(field)

769 initial_stars_schema.addField(item.getField())

770 id_type = self.psf_schema["id"].asField().getTypeString()

771 psf_max_value_type = self.psf_schema['psf_max_value'].asField().getTypeString()

772 initial_stars_schema.addField("psf_id",

773 type=id_type,

774 doc="id of this source in psf_stars; 0 if there is no match.")

775 initial_stars_schema.addField("psf_max_value",

776 type=psf_max_value_type,

777 doc="Maximum value in the star image used to train PSF.")

778

779 afwTable.CoordKey.addErrorFields(initial_stars_schema)

780 self.makeSubtask("star_background")

781 self.makeSubtask("star_detection", schema=initial_stars_schema)

782 self.makeSubtask("star_sky_sources", schema=initial_stars_schema)

783 self.makeSubtask("star_deblend", schema=initial_stars_schema)

784 self.makeSubtask("star_measurement", schema=initial_stars_schema)

785 self.makeSubtask("star_normalized_calibration_flux", schema=initial_stars_schema)

786

787 self.makeSubtask("star_apply_aperture_correction", schema=initial_stars_schema)

788 self.makeSubtask("star_catalog_calculation", schema=initial_stars_schema)

789 self.makeSubtask("star_set_primary_flags", schema=initial_stars_schema, isSingleFrame=True)

790 self.makeSubtask("star_selector")

791 self.makeSubtask("photometry", schema=initial_stars_schema)

792 if self.config.doMaskDiffractionSpikes:

793 self.makeSubtask("diffractionSpikeMask")

794 self.makeSubtask("compute_summary_stats")

795

796 # The final catalog will have calibrated flux columns, which we add to

797 # the init-output schema by calibrating our zero-length catalog with an

798 # arbitrary dummy PhotoCalib. We also use this schema to initialize

799 # the stars catalog in order to ensure it's the same even when we hit

800 # an error (and write partial outputs) before calibrating the catalog

801 # - note that calibrateCatalog will happily reuse existing output

802 # columns.

803 dummy_photo_calib = afwImage.PhotoCalib(1.0, 0, bbox=lsst.geom.Box2I())

804 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(

805 afwTable.SourceCatalog(initial_stars_schema)

806 )

807

808 def runQuantum(self, butlerQC, inputRefs, outputRefs):

809 inputs = butlerQC.get(inputRefs)

810 exposures = inputs.pop("exposures")

811

812 exposure_record = inputRefs.exposures[0].dataId.records["exposure"]

813 if self.config.useButlerExposureRegion:

814 exposure_region = butlerQC.quantum.dataId.region

815 else:

816 exposure_region = None

817

818 id_generator = self.config.id_generator.apply(butlerQC.quantum.dataId)

819

820 astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(

821 dataIds=[ref.datasetRef.dataId for ref in inputRefs.astrometry_ref_cat],

822 refCats=inputs.pop("astrometry_ref_cat"),

823 name=self.config.connections.astrometry_ref_cat,

824 config=self.config.astrometry_ref_loader, log=self.log)

825 self.astrometry.setRefObjLoader(astrometry_loader)

826

827 photometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(

828 dataIds=[ref.datasetRef.dataId for ref in inputRefs.photometry_ref_cat],

829 refCats=inputs.pop("photometry_ref_cat"),

830 name=self.config.connections.photometry_ref_cat,

831 config=self.config.photometry_ref_loader, log=self.log)

832 self.photometry.match.setRefObjLoader(photometry_loader)

833

834 if self.config.doMaskDiffractionSpikes:

835 # Use the same photometry reference catalog for the bright star

836 # mask.

837 self.diffractionSpikeMask.setRefObjLoader(photometry_loader)

838

839 if self.config.do_illumination_correction:

840 background_flat = inputs.pop("background_flat")

841 illumination_correction = inputs.pop("illumination_correction")

842 else:

843 background_flat = None

844 illumination_correction = None

845

846 camera_model = None

847 if self.config.useButlerCamera:

848 if "camera_model" in inputs:

849 camera_model = inputs.pop("camera_model")

850 else:

851 self.log.warning("useButlerCamera=True, but camera is not available for filter %s. The "

852 "astrometry fit will use the WCS already attached to the exposure.",

853 exposures[0].filter.bandLabel)

854

855 # This should not happen with a properly configured execution context.

856 assert not inputs, "runQuantum got more inputs than expected"

857

858 # Specify the fields that `annotate` needs below, to ensure they

859 # exist, even as None.

860 result = pipeBase.Struct(

861 exposure=None,

862 stars_footprints=None,

863 psf_stars_footprints=None,

864 background_to_photometric_ratio=None,

865 )

866 try:

867 self.run(

868 exposures=exposures,

869 result=result,

870 id_generator=id_generator,

871 background_flat=background_flat,

872 illumination_correction=illumination_correction,

873 camera_model=camera_model,

874 exposure_record=exposure_record,

875 exposure_region=exposure_region,

876 )

877 except pipeBase.AlgorithmError as e:

878 error = pipeBase.AnnotatedPartialOutputsError.annotate(

879 e,

880 self,

881 result.exposure,

882 result.psf_stars_footprints,

883 result.stars_footprints,

884 log=self.log

885 )

886 butlerQC.put(result, outputRefs)

887 raise error from e

888

889 butlerQC.put(result, outputRefs)

890

891 @timeMethod

892 def run(

893 self,

894 *,

895 exposures,

896 id_generator=None,

897 result=None,

898 background_flat=None,

899 illumination_correction=None,

900 camera_model=None,

901 exposure_record=None,

902 exposure_region=None,

903 ):

904 """Find stars and perform psf measurement, then do a deeper detection

905 and measurement and calibrate astrometry and photometry from that.

906

907 Parameters

908 ----------

909 exposures : `lsst.afw.image.Exposure` or \

910 `list` [`lsst.afw.image.Exposure`]

911 Post-ISR exposure(s), with an initial WCS, VisitInfo, and Filter.

912 Modified in-place during processing if only one is passed.

913 If two exposures are passed, treat them as snaps and combine

914 before doing further processing.

915 id_generator : `lsst.meas.base.IdGenerator`, optional

916 Object that generates source IDs and provides random seeds.

917 result : `lsst.pipe.base.Struct`, optional

918 Result struct that is modified to allow saving of partial outputs

919 for some failure conditions. If the task completes successfully,

920 this is also returned.

921 background_flat : `lsst.afw.image.Exposure`, optional

922 Background flat-field image.

923 illumination_correction : `lsst.afw.image.Exposure`, optional

924 Illumination correction image.

925 camera_model : `lsst.afw.cameraGeom.Camera`, optional

926 Camera to be used if constructing updated WCS.

927 exposure_record : `lsst.daf.butler.DimensionRecord`, optional

928 Butler metadata for the ``exposure`` dimension. Used if

929 constructing an updated WCS instead of the boresight and

930 orientation angle from the visit info.

931 exposure_region : `lsst.sphgeom.Region`, optional

932 The exposure region to use for the reference catalog filtering.

933 If `None`, this region will be set as a padded bbox + current WCS

934 of the exposure.

935

936 Returns

937 -------

938 result : `lsst.pipe.base.Struct`

939 Results as a struct with attributes:

940

941 ``exposure``

942 Calibrated exposure, with pixels in nJy units.

943 (`lsst.afw.image.Exposure`)

944 ``stars``

945 Stars that were used to calibrate the exposure, with

946 calibrated fluxes and magnitudes.

947 (`astropy.table.Table`)

948 ``stars_footprints``

949 Footprints of stars that were used to calibrate the exposure.

950 (`lsst.afw.table.SourceCatalog`)

951 ``psf_stars``

952 Stars that were used to determine the image PSF.

953 (`astropy.table.Table`)

954 ``psf_stars_footprints``

955 Footprints of stars that were used to determine the image PSF.

956 (`lsst.afw.table.SourceCatalog`)

957 ``background``

958 Background that was fit to the exposure when detecting

959 ``stars``. (`lsst.afw.math.BackgroundList`)

960 ``applied_photo_calib``

961 Photometric calibration that was fit to the star catalog and

962 applied to the exposure. (`lsst.afw.image.PhotoCalib`)

963 This is `None` if ``config.do_calibrate_pixels`` is `False`.

964 ``astrometry_matches``

965 Reference catalog stars matches used in the astrometric fit.

966 (`list` [`lsst.afw.table.ReferenceMatch`] or

967 `lsst.afw.table.BaseCatalog`).

968 ``photometry_matches``

969 Reference catalog stars matches used in the photometric fit.

970 (`list` [`lsst.afw.table.ReferenceMatch`] or

971 `lsst.afw.table.BaseCatalog`).

972 ``mask``

973 Copy of the mask plane of `exposure`.

974 (`lsst.afw.image.Mask`)

975 """

976 if result is None:

977 result = pipeBase.Struct()

978 if id_generator is None:

979 id_generator = lsst.meas.base.IdGenerator()

980

981 result.exposure = self.snap_combine.run(exposures).exposure

982 self.log.info("Initial PhotoCalib: %s", result.exposure.getPhotoCalib())

983

984 result.exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = False

985

986 # Check input image processing.

987 if self.config.do_illumination_correction:

988 if not result.exposure.metadata.get("LSST ISR FLAT APPLIED", False):

989 raise pipeBase.InvalidQuantumError(

990 "Cannot use do_illumination_correction with an image that has not had a flat applied",

991 )

992

993 # Update the exposure pointing to that stored in the butler record

994 # (regardless of a camera model update). This is to pick up any

995 # updates to the pointing stored in the butler record that may not

996 # be reflected in the exposure metadata (headers).

997 if camera_model:

998 self._update_wcs_with_camera_model(

999 result.exposure, camera_model, exposure_record=exposure_record

1000 )

1001 elif exposure_record is not None:

1002 self._update_wcs_to_exposure_record(result.exposure, exposure_record)

1003

1004 result.background = None

1005 result.background_to_photometric_ratio = None

1006 summary_stat_catalog = None

1007 # Some exposure components are set to initial placeholder objects

1008 # while we try to bootstrap them. If we fail before we fit for them,

1009 # we want to reset those components to None so the placeholders don't

1010 # masquerade as the real thing.

1011 have_fit_psf = False

1012 have_fit_astrometry = False

1013 have_fit_photometry = False

1014 try:

1015 result.background_to_photometric_ratio = self._apply_illumination_correction(

1016 result.exposure,

1017 background_flat,

1018 illumination_correction,

1019 )

1020

1021 result.psf_stars_footprints, _, adaptive_det_res_struct = self._compute_psf(

1022 result,

1023 id_generator,

1024 )

1025 have_fit_psf = True

1026

1027 # Check if all centroids have been flagged. This should happen

1028 # after the PSF is fit and recorded because even a junky PSF

1029 # model is informative.

1030 if result.psf_stars_footprints["slot_Centroid_flag"].all():

1031 psf_shape = result.exposure.psf.computeShape(result.exposure.psf.getAveragePosition())

1032 raise AllCentroidsFlaggedError(

1033 n_sources=len(result.psf_stars_footprints),

1034 psf_shape_ixx=psf_shape.getIxx(),

1035 psf_shape_iyy=psf_shape.getIyy(),

1036 psf_shape_ixy=psf_shape.getIxy(),

1037 psf_size=psf_shape.getDeterminantRadius(),

1038 )

1039

1040 if result.background is None:

1041 result.background = afwMath.BackgroundList()

1042

1043 self._measure_aperture_correction(result.exposure, result.psf_stars_footprints)

1044 result.psf_stars = result.psf_stars_footprints.asAstropy()

1045 # Run astrometry using PSF candidate stars.

1046 # Update "the psf_stars" source coordinates with the current wcs.

1047 afwTable.updateSourceCoords(

1048 result.exposure.wcs,

1049 sourceList=result.psf_stars_footprints,

1050 include_covariance=self.config.do_include_astrometric_errors,

1051 )

1052 astrometry_matches, astrometry_meta = self._fit_astrometry(

1053 result.exposure, result.psf_stars_footprints, exposure_region=exposure_region

1054 )

1055 num_astrometry_matches = len(astrometry_matches)

1056 if "astrometry_matches" in self.config.optional_outputs:

1057 result.astrometry_matches = lsst.meas.astrom.denormalizeMatches(astrometry_matches,

1058 astrometry_meta)

1059 result.psf_stars = result.psf_stars_footprints.asAstropy()

1060

1061 # Add diffraction spike mask here so it can be used (i.e. avoided)

1062 # in the adaptive background estimation.

1063 if self.config.doMaskDiffractionSpikes:

1064 self.diffractionSpikeMask.run(result.exposure)

1065

1066 self.metadata['adaptive_threshold_value'] = float("nan")

1067 if self.config.do_adaptive_threshold_detection:

1068 self._remeasure_star_background(

1069 result,

1070 background_to_photometric_ratio=result.background_to_photometric_ratio,

1071 )

1072

1073 # Run the stars_detection subtask for the photometric calibration.

1074 result.stars_footprints = self._find_stars(

1075 result.exposure,

1076 result.background,

1077 id_generator,

1078 background_to_photometric_ratio=result.background_to_photometric_ratio,

1079 adaptive_det_res_struct=adaptive_det_res_struct,

1080 num_astrometry_matches=num_astrometry_matches,

1081 )

1082 psf = result.exposure.getPsf()

1083 psfSigma = psf.computeShape(result.exposure.getBBox().getCenter()).getDeterminantRadius()

1084 self._match_psf_stars(result.psf_stars_footprints, result.stars_footprints,

1085 psfSigma=psfSigma)

1086

1087 # Update the "stars" source coordinates with the current wcs.

1088 afwTable.updateSourceCoords(

1089 result.exposure.wcs,

1090 sourceList=result.stars_footprints,

1091 include_covariance=self.config.do_include_astrometric_errors,

1092 )

1093

1094 summary_stat_catalog = result.stars_footprints

1095 result.stars = result.stars_footprints.asAstropy()

1096 self.metadata["star_count"] = np.sum(~result.stars["sky_source"])

1097

1098 # Validate the astrometric fit. Send in the stars_footprints

1099 # catalog so that its coords get set to NaN if the fit is deemed

1100 # a failure.

1101 self.astrometry.check(result.exposure, result.stars_footprints, len(astrometry_matches))

1102 result.stars = result.stars_footprints.asAstropy()

1103 have_fit_astrometry = True

1104

1105 result.stars_footprints, photometry_matches, \

1106 photometry_meta, photo_calib = self._fit_photometry(result.exposure, result.stars_footprints)

1107 have_fit_photometry = True

1108 self.metadata["photometry_matches_count"] = len(photometry_matches)

1109 # fit_photometry returns a new catalog, so we need a new astropy

1110 # table view.

1111 result.stars = result.stars_footprints.asAstropy()

1112 # summary stats don't make use of the calibrated fluxes, but we

1113 # might as well use the best catalog we've got in case that

1114 # changes, and help the old one get garbage-collected.

1115 summary_stat_catalog = result.stars_footprints

1116 if "photometry_matches" in self.config.optional_outputs:

1117 result.photometry_matches = lsst.meas.astrom.denormalizeMatches(photometry_matches,

1118 photometry_meta)

1119 if "mask" in self.config.optional_outputs:

1120 result.mask = result.exposure.mask.clone()

1121 except pipeBase.AlgorithmError:

1122 if not have_fit_psf:

1123 result.exposure.setPsf(None)

1124 if not have_fit_astrometry:

1125 result.exposure.setWcs(None)

1126 if not have_fit_photometry:

1127 result.exposure.setPhotoCalib(None)

1128 # Summary stat calculations can handle missing components

1129 # gracefully, but we want to run them as late as possible (but

1130 # still before we calibrate pixels, if we do that at all).

1131 # So we run them after we succeed or if we get an AlgorithmError.

1132 # We intentionally don't use 'finally' here because we don't

1133 # want to run them if we get some other kind of error.

1134 self._summarize(result.exposure, summary_stat_catalog, result.background)

1135 raise

1136 else:

1137 self._summarize(result.exposure, summary_stat_catalog, result.background)

1138

1139 # Output post-subtraction background stats to task metadata:

1140 # Specify the pixels flags to ignore, starting with those ignored

1141 # by the subtraction.

1142 bgIgnoredPixelMasks = lsst.meas.algorithms.SubtractBackgroundConfig().ignoredPixelMask.list()

1143 for maskName in self.config.background_stats_ignored_pixel_masks:

1144 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()

1145 and maskName not in bgIgnoredPixelMasks):

1146 bgIgnoredPixelMasks += [maskName]

1147

1148 statsCtrl = afwMath.StatisticsControl()

1149 statsCtrl.setAndMask(afwImage.Mask.getPlaneBitMask(bgIgnoredPixelMasks))

1150

1151 statObj = afwMath.makeStatistics(result.exposure.getMaskedImage(), afwMath.MEDIAN, statsCtrl)

1152 median_bg, _ = statObj.getResult(afwMath.MEDIAN)

1153 self.metadata["bg_subtracted_skyPixel_instFlux_median"] = median_bg

1154

1155 statObj = afwMath.makeStatistics(result.exposure.getMaskedImage(), afwMath.STDEV, statsCtrl)

1156 stdev_bg, _ = statObj.getResult(afwMath.STDEV)

1157 self.metadata["bg_subtracted_skyPixel_instFlux_stdev"] = stdev_bg

1158

1159 self.metadata["bg_subtracted_skySource_flux_median"] = (

1160 np.median(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column])

1161 )

1162 self.metadata["bg_subtracted_skySource_flux_stdev"] = (

1163 np.std(result.stars[result.stars['sky_source']][self.config.background_stats_flux_column])

1164 )

1165

1166 if self.config.do_calibrate_pixels:

1167 self._apply_photometry(

1168 result.exposure,

1169 result.background,

1170 background_to_photometric_ratio=result.background_to_photometric_ratio,

1171 )

1172 result.applied_photo_calib = photo_calib

1173 else:

1174 result.applied_photo_calib = None

1175

1176 self._recordMaskedPixelFractions(result.exposure)

1177

1178 if self.config.run_sattle:

1179 # send boresight and timing information to sattle so the cache

1180 # is populated by the time we reach ip_diffim detectAndMeasure.

1181 try:

1182 populate_sattle_visit_cache(result.exposure.getInfo().getVisitInfo(),

1183 historical=self.config.sattle_historical)

1184 self.log.info('Successfully triggered load of sattle visit cache')

1185 except requests.exceptions.HTTPError:

1186 self.log.exception("Sattle visit cache update failed; continuing with image processing")

1187

1188 return result

1189

1190 def _apply_illumination_correction(self, exposure, background_flat, illumination_correction):

1191 """Apply the illumination correction to a background-flattened image.

1192

1193 Parameters

1194 ----------

1195 exposure : `lsst.afw.image.Exposure`

1196 Exposure to convert to a photometric-flattened image.

1197 background_flat : `lsst.afw.image.Exposure`

1198 Flat image that had previously been applied to exposure.

1199 illumination_correction : `lsst.afw.image.Exposure`

1200 Illumination correction image to convert to photometric-flattened

1201 image.

1202

1203 Returns

1204 -------

1205 background_to_photometric_ratio : `lsst.afw.image.Image`

1206 Ratio image to convert a photometric-flattened image to/from

1207 a background-flattened image. Will be None if task not

1208 configured to use the illumination correction.

1209 """

1210 if not self.config.do_illumination_correction:

1211 return None

1212

1213 # From a raw image to a background-flattened image, we have:

1214 # bfi = image / background_flat

1215 # From a raw image to a photometric-flattened image, we have:

1216 # pfi = image / reference_flux_flat

1217 # pfi = image / (dome_flat * illumination_correction),

1218 # where the illumination correction contains the jacobian

1219 # of the wcs, converting to fluence units.

1220 # Currently background_flat == dome_flat, so we have for the

1221 # "background_to_photometric_ratio", the ratio of the background-

1222 # flattened image to the photometric-flattened image:

1223 # bfi / pfi = illumination_correction.

1224

1225 background_to_photometric_ratio = illumination_correction.image.clone()

1226

1227 # Dividing the ratio will convert a background-flattened image to

1228 # a photometric-flattened image.

1229 exposure.maskedImage /= background_to_photometric_ratio

1230

1231 exposure.metadata["LSST CALIB ILLUMCORR APPLIED"] = True

1232

1233 return background_to_photometric_ratio

1234

1235 def _compute_psf(self, result, id_generator):

1236 """Find bright sources detected on an exposure and fit a PSF model to

1237 them, repairing likely cosmic rays before detection.

1238

1239 Repair, detect, measure, and compute PSF twice, to ensure the PSF

1240 model does not include contributions from cosmic rays.

1241

1242 Parameters

1243 ----------

1244 result : `lsst.pipe.base.Struct`

1245 Result struct that is modified to allow saving of partial outputs

1246 for some failure conditions. Should contain at least the following

1247 attributes:

1248

1249 - exposure : `lsst.afw.image.Exposure`

1250 Exposure to detect and measure bright stars on.

1251 - background : `lsst.afw.math.BackgroundList` | `None`

1252 Background that was fit to the exposure during detection.

1253 - background_to_photometric_ratio : `lsst.afw.image.Image` | `None`

1254 Image to convert photometric-flattened image to

1255 background-flattened image.

1256 id_generator : `lsst.meas.base.IdGenerator`

1257 Object that generates source IDs and provides random seeds.

1258

1259 Returns

1260 -------

1261 sources : `lsst.afw.table.SourceCatalog`

1262 Catalog of detected bright sources.

1263 cell_set : `lsst.afw.math.SpatialCellSet`

1264 PSF candidates returned by the psf determiner.

1265 adaptive_det_res_struct : `lsst.pipe.base.Struct`

1266 Result struct from the adaptive threshold detection.

1267

1268 Notes

1269 -----

1270 This method modifies the exposure, background and

1271 background_to_photometric_ratio attributes of the result struct

1272 in-place.

1273 """

1274 def log_psf(msg, addToMetadata=False):

1275 """Log the parameters of the psf and background, with a prepended

1276 message. There is also the option to add the PSF sigma to the task

1277 metadata.

1278

1279 Parameters

1280 ----------

1281 msg : `str`

1282 Message to prepend the log info with.

1283 addToMetadata : `bool`, optional

1284 Whether to add the final psf sigma value to the task

1285 metadata (the default is False).

1286 """

1287 position = result.exposure.psf.getAveragePosition()

1288 sigma = result.exposure.psf.computeShape(position).getDeterminantRadius()

1289 dimensions = result.exposure.psf.computeImage(position).getDimensions()

1290 if result.background is not None:

1291 median_background = np.median(result.background.getImage().array)

1292 else:

1293 median_background = 0.0

1294 self.log.info("%s sigma=%0.4f, dimensions=%s; median background=%0.2f",

1295 msg, sigma, dimensions, median_background)

1296 if addToMetadata:

1297 self.metadata["final_psf_sigma"] = sigma

1298

1299 self.log.info("First pass detection with Gaussian PSF FWHM=%s pixels",

1300 self.config.install_simple_psf.fwhm)

1301 self.install_simple_psf.run(exposure=result.exposure)

1302

1303 result.background = self.psf_subtract_background.run(

1304 exposure=result.exposure,

1305 backgroundToPhotometricRatio=result.background_to_photometric_ratio,

1306 ).background

1307 log_psf("Initial PSF:")

1308 self.psf_repair.run(exposure=result.exposure, keepCRs=True)

1309

1310 table = afwTable.SourceTable.make(self.psf_schema, id_generator.make_table_id_factory())

1311 if not self.config.do_adaptive_threshold_detection:

1312 adaptive_det_res_struct = None

1313 # Re-estimate the background during this detection step, so that

1314 # measurement uses the most accurate background-subtraction.

1315 detections = self.psf_detection.run(

1316 table=table,

1317 exposure=result.exposure,

1318 background=result.background,

1319 backgroundToPhotometricRatio=result.background_to_photometric_ratio,

1320 )

1321 else:

1322 initialThreshold = self.config.psf_detection.thresholdValue

1323 initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier

1324 adaptive_det_res_struct = self.psf_adaptive_threshold_detection.run(

1325 table, result.exposure,

1326 initialThreshold=initialThreshold,

1327 initialThresholdMultiplier=initialThresholdMultiplier,

1328 )

1329 detections = adaptive_det_res_struct.detections

1330

1331 self.metadata["initial_psf_positive_footprint_count"] = detections.numPos

1332 self.metadata["initial_psf_negative_footprint_count"] = detections.numNeg

1333 self.metadata["initial_psf_positive_peak_count"] = detections.numPosPeaks

1334 self.metadata["initial_psf_negative_peak_count"] = detections.numNegPeaks

1335 self.psf_source_measurement.run(detections.sources, result.exposure)

1336 psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)

1337

1338 # This 2nd round of PSF fitting has been deemed unnecessary (and

1339 # sometimes even causing harm), so it is being skipped for the

1340 # adaptive threshold logic, but is left here for now to maintain

1341 # consistency with the previous logic for any users wanting to

1342 # revert to it.

1343 if not self.config.do_adaptive_threshold_detection:

1344 # Replace the initial PSF with something simpler for the second

1345 # repair/detect/measure/measure_psf step: this can help it

1346 # converge.

1347 self.install_simple_psf.run(exposure=result.exposure)

1348

1349 log_psf("Rerunning with simple PSF:")

1350 # TODO investigation: Should we only re-run repair here, to use the

1351 # new PSF? Maybe we *do* need to re-run measurement with PsfFlux,

1352 # to use the fitted PSF?

1353 # TODO investigation: do we need a separate measurement task here

1354 # for the post-psf_measure_psf step, since we only want to do

1355 # PsfFlux and GaussianFlux *after* we have a PSF? Maybe that's not

1356 # relevant once DM-39203 is merged?

1357 self.psf_repair.run(exposure=result.exposure, keepCRs=True)

1358 # Re-estimate the background during this detection step, so that

1359 # measurement uses the most accurate background-subtraction.

1360 detections = self.psf_detection.run(

1361 table=table,

1362 exposure=result.exposure,

1363 background=result.background,

1364 backgroundToPhotometricRatio=result.background_to_photometric_ratio,

1365 )

1366 self.psf_source_measurement.run(detections.sources, result.exposure)

1367 psf_result = self.psf_measure_psf.run(exposure=result.exposure, sources=detections.sources)

1368 self.metadata["simple_psf_positive_footprint_count"] = detections.numPos

1369 self.metadata["simple_psf_negative_footprint_count"] = detections.numNeg

1370 self.metadata["simple_psf_positive_peak_count"] = detections.numPosPeaks

1371 self.metadata["simple_psf_negative_peak_count"] = detections.numNegPeaks

1372 log_psf("Final PSF:", addToMetadata=True)

1373

1374 # Final repair with final PSF, removing cosmic rays this time.

1375 self.psf_repair.run(exposure=result.exposure)

1376 # Final measurement with the CRs removed.

1377 self.psf_source_measurement.run(detections.sources, result.exposure)

1378

1379 # PSF is set on exposure; candidates are returned to use for

1380 # calibration flux normalization and aperture corrections.

1381 return detections.sources, psf_result.cellSet, adaptive_det_res_struct

1382

1383 def _measure_aperture_correction(self, exposure, bright_sources):

1384 """Measure and set the ApCorrMap on the Exposure, using

1385 previously-measured bright sources.

1386

1387 This function first normalizes the calibration flux and then

1388 the full set of aperture corrections are measured relative

1389 to this normalized calibration flux.

1390

1391 Parameters

1392 ----------

1393 exposure : `lsst.afw.image.Exposure`

1394 Exposure to set the ApCorrMap on.

1395 bright_sources : `lsst.afw.table.SourceCatalog`

1396 Catalog of detected bright sources; modified to include columns

1397 necessary for point source determination for the aperture

1398 correction calculation.

1399 """

1400 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(

1401 exposure=exposure,

1402 catalog=bright_sources,

1403 ).ap_corr_map

1404

1405 ap_corr_map = self.measure_aperture_correction.run(exposure, bright_sources).apCorrMap

1406

1407 # Need to merge the aperture correction map from the normalization.

1408 for key in norm_ap_corr_map:

1409 ap_corr_map[key] = norm_ap_corr_map[key]

1410

1411 exposure.info.setApCorrMap(ap_corr_map)

1412

1413 def _find_stars(self, exposure, background, id_generator, background_to_photometric_ratio=None,

1414 adaptive_det_res_struct=None, num_astrometry_matches=None):

1415 """Detect stars on an exposure that has a PSF model, and measure their

1416 PSF, circular aperture, compensated gaussian fluxes.

1417

1418 Parameters

1419 ----------

1420 exposure : `lsst.afw.image.Exposure`

1421 Exposure to detect and measure stars on.

1422 background : `lsst.afw.math.BackgroundList`

1423 Background that was fit to the exposure during detection;

1424 modified in-place during subsequent detection.

1425 id_generator : `lsst.meas.base.IdGenerator`

1426 Object that generates source IDs and provides random seeds.

1427 background_to_photometric_ratio : `lsst.afw.image.Image`, optional

1428 Image to convert photometric-flattened image to

1429 background-flattened image.

1430

1431 Returns

1432 -------

1433 stars : `SourceCatalog`

1434 Sources that are very likely to be stars, with a limited set of

1435 measurements performed on them.

1436 """

1437 table = afwTable.SourceTable.make(self.initial_stars_schema.schema,

1438 id_generator.make_table_id_factory())

1439

1440 maxAdaptiveDetIter = 8

1441 threshFactor = 0.2

1442 if adaptive_det_res_struct is not None:

1443 for adaptiveDetIter in range(maxAdaptiveDetIter):

1444 adaptiveDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()

1445 adaptiveDetectionConfig.reEstimateBackground = False

1446 adaptiveDetectionConfig.includeThresholdMultiplier = 2.0

1447 psfThreshold = (

1448 adaptive_det_res_struct.thresholdValue*adaptive_det_res_struct.includeThresholdMultiplier

1449 )

1450 adaptiveDetectionConfig.thresholdValue = max(

1451 self.config.star_detection.thresholdValue,

1452 threshFactor*psfThreshold/adaptiveDetectionConfig.includeThresholdMultiplier

1453 )

1454 self.log.info("Using adaptive threshold detection (nIter = %d) with "

1455 "thresholdValue = %.2f and multiplier = %.1f",

1456 adaptiveDetIter, adaptiveDetectionConfig.thresholdValue,

1457 adaptiveDetectionConfig.includeThresholdMultiplier)

1458 adaptiveDetectionTask = lsst.meas.algorithms.SourceDetectionTask(

1459 config=adaptiveDetectionConfig

1460 )

1461 detections = adaptiveDetectionTask.run(

1462 table=table,

1463 exposure=exposure,

1464 background=background,

1465 backgroundToPhotometricRatio=background_to_photometric_ratio,

1466 )

1467 nFootprint = len(detections.sources)

1468 nPeak = 0

1469 nIsolated = 0

1470 for src in detections.sources:

1471 nPeakSrc = len(src.getFootprint().getPeaks())

1472 if nPeakSrc == 1:

1473 nIsolated += 1

1474 nPeak += nPeakSrc

1475 minIsolated = min(400, max(3, 0.005*nPeak, 0.6*num_astrometry_matches))

1476 if nFootprint > 0:

1477 self.log.info("Adaptive threshold detection _find_stars nIter: %d; "

1478 "nPeak/nFootprint = %.2f (max is 800); nIsolated = %d (min is %.1f).",

1479 adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated)

1480 if nPeak/nFootprint > 800 or nIsolated < minIsolated:

1481 threshFactor = max(0.01, 1.5*threshFactor)

1482 self.log.warning("nPeak/nFootprint = %.2f (max is 800); nIsolated = %d "

1483 "(min is %.1f).", nPeak/nFootprint, nIsolated, minIsolated)

1484 else:

1485 break

1486 else:

1487 threshFactor *= 0.75

1488 self.log.warning("No footprints detected on image. Decreasing threshold "

1489 "factor to %.2f. and rerunning.", threshFactor)

1490 else:

1491 # Re-estimate the background during this detection step, so that

1492 # measurement uses the most accurate background-subtraction.

1493 detections = self.star_detection.run(

1494 table=table,

1495 exposure=exposure,

1496 background=background,

1497 backgroundToPhotometricRatio=background_to_photometric_ratio,

1498 )

1499 sources = detections.sources

1500 self.star_sky_sources.run(exposure.mask, id_generator.catalog_id, sources)

1501

1502 n_sky_sources = np.sum(sources["sky_source"])

1503 if (self.config.do_downsample_footprints

1504 and (len(sources) - n_sky_sources) > self.config.downsample_max_footprints):

1505 if exposure.info.id is None:

1506 self.log.warning("Exposure does not have a proper id; using 0 seed for downsample.")

1507 seed = 0

1508 else:

1509 seed = exposure.info.id & 0xFFFFFFFF

1510

1511 gen = np.random.RandomState(seed)

1512

1513 # Isolate the sky sources before downsampling.

1514 indices = np.arange(len(sources))[~sources["sky_source"]]

1515 indices = gen.choice(

1516 indices,

1517 size=self.config.downsample_max_footprints,

1518 replace=False,

1519 )

1520 skyIndices, = np.where(sources["sky_source"])

1521 indices = np.concatenate((indices, skyIndices))

1522

1523 self.log.info("Downsampling from %d to %d non-sky-source footprints.", len(sources), len(indices))

1524

1525 sel = np.zeros(len(sources), dtype=bool)

1526 sel[indices] = True

1527 sources = sources[sel]

1528

1529 # TODO investigation: Could this deblender throw away blends of

1530 # non-PSF sources?

1531 self.star_deblend.run(exposure=exposure, sources=sources)

1532 # The deblender may not produce a contiguous catalog; ensure

1533 # contiguity for subsequent tasks.

1534 if not sources.isContiguous():

1535 sources = sources.copy(deep=True)

1536

1537 # Measure everything, and use those results to select only stars.

1538 self.star_measurement.run(sources, exposure)

1539 self.metadata["post_deblend_source_count"] = np.sum(~sources["sky_source"])

1540 self.metadata["failed_deblend_source_count"] = np.sum(

1541 ~sources["sky_source"] & sources["deblend_failed"]

1542 )

1543 self.metadata["saturated_source_count"] = np.sum(sources["base_PixelFlags_flag_saturated"])

1544 self.metadata["bad_source_count"] = np.sum(sources["base_PixelFlags_flag_bad"])

1545

1546 # Run the normalization calibration flux task to apply the

1547 # normalization correction to create normalized

1548 # calibration fluxes.

1549 self.star_normalized_calibration_flux.run(exposure=exposure, catalog=sources)

1550 self.star_apply_aperture_correction.run(sources, exposure.apCorrMap)

1551 self.star_catalog_calculation.run(sources)

1552 self.star_set_primary_flags.run(sources)

1553

1554 result = self.star_selector.run(sources)

1555 # The star selector may not produce a contiguous catalog.

1556 if not result.sourceCat.isContiguous():

1557 return result.sourceCat.copy(deep=True)

1558 else:

1559 return result.sourceCat

1560

1561 def _match_psf_stars(self, psf_stars, stars, psfSigma=None):

1562 """Match calibration stars to psf stars, to identify which were psf

1563 candidates, and which were used or reserved during psf measurement

1564 and the astrometric fit.

1565

1566 Parameters

1567 ----------

1568 psf_stars : `lsst.afw.table.SourceCatalog`

1569 PSF candidate stars that were sent to the psf determiner and

1570 used in the astrometric fit. Used to populate psf and astrometry

1571 related flag fields.

1572 stars : `lsst.afw.table.SourceCatalog`

1573 Stars that will be used for calibration; psf-related fields will

1574 be updated in-place.

1575

1576 Notes

1577 -----

1578 This code was adapted from CalibrateTask.copyIcSourceFields().

1579 """

1580 control = afwTable.MatchControl()

1581 matchRadius = 3.0 if psfSigma is None else max(3.0, psfSigma) # in pixels

1582 # Return all matched objects, to separate blends.

1583 control.findOnlyClosest = False

1584 matches = afwTable.matchXy(psf_stars, stars, matchRadius, control)

1585 deblend_key = stars.schema["deblend_nChild"].asKey()

1586 matches = [m for m in matches if m[1].get(deblend_key) == 0]

1587

1588 # Because we had to allow multiple matches to handle parents, we now

1589 # need to prune to the best (closest) matches.

1590 # Closest matches is a dict of psf_stars source ID to Match record

1591 # (psf_stars source, sourceCat source, distance in pixels).

1592 best = {}

1593 for match_psf, match_stars, d in matches:

1594 match = best.get(match_psf.getId())

1595 if match is None or d <= match[2]:

1596 best[match_psf.getId()] = (match_psf, match_stars, d)

1597 matches = list(best.values())

1598 # We'll use this to construct index arrays into each catalog.

1599 ids = np.array([(match_psf.getId(), match_stars.getId()) for match_psf, match_stars, d in matches]).T

1600

1601 if (n_matches := len(matches)) == 0:

1602 raise NoPsfStarsToStarsMatchError(n_psf_stars=len(psf_stars), n_stars=len(stars))

1603

1604 self.log.info("%d psf/astrometry stars out of %d matched %d calib stars",

1605 n_matches, len(psf_stars), len(stars))

1606 self.metadata["matched_psf_star_count"] = n_matches

1607

1608 # Check that no stars sources are listed twice; we already know

1609 # that each match has a unique psf_stars id, due to using as the key

1610 # in best above.

1611 n_unique = len(set(m[1].getId() for m in matches))

1612 if n_unique != n_matches:

1613 self.log.warning("%d psf_stars matched only %d stars", n_matches, n_unique)

1614

1615 # The indices of the IDs, so we can update the flag fields as arrays.

1616 idx_psf_stars = np.searchsorted(psf_stars["id"], ids[0])

1617 idx_stars = np.searchsorted(stars["id"], ids[1])

1618 for field in self.psf_fields:

1619 result = np.zeros(len(stars), dtype=bool)

1620 result[idx_stars] = psf_stars[field][idx_psf_stars]

1621 stars[field] = result

1622 stars['psf_id'][idx_stars] = psf_stars['id'][idx_psf_stars]

1623 stars['psf_max_value'][idx_stars] = psf_stars['psf_max_value'][idx_psf_stars]

1624

1625 def _fit_astrometry(self, exposure, stars, exposure_region=None):

1626 """Fit an astrometric model to the data and return the reference

1627 matches used in the fit, and the fitted WCS.

1628

1629 Parameters

1630 ----------

1631 exposure : `lsst.afw.image.Exposure`

1632 Exposure that is being fit, to get PSF and other metadata from.

1633 Modified to add the fitted skyWcs.

1634 stars : `SourceCatalog`

1635 Good stars selected for use in calibration, with RA/Dec coordinates

1636 computed from the pixel positions and fitted WCS.

1637 exposure_region : `lsst.sphgeom.Region`, optional

1638 The exposure region to use for the reference catalog filtering.

1639 If `None`, this region will be set as a padded bbox + current WCS

1640 of the exposure.

1641

1642 Returns

1643 -------

1644 matches : `list` [`lsst.afw.table.ReferenceMatch`]

1645 Reference/stars matches used in the fit.

1646 """

1647 initialWcs = exposure.wcs

1648 result = self.astrometry.run(stars, exposure, exposure_region=exposure_region)

1649

1650 # Record astrometry stats to metadata

1651 self.metadata["astrometry_matches_count"] = len(result.matches)

1652 if exposure.wcs is not None:

1653 if initialWcs is not None:

1654 center = exposure.getBBox().getCenter()

1655 self.metadata['initial_to_final_wcs'] = (

1656 initialWcs.pixelToSky(center).separation(

1657 exposure.wcs.pixelToSky(center)

1658 ).asArcseconds()

1659 )

1660 else:

1661 self.metadata['initial_to_final_wcs'] = float("nan")

1662 self.metadata['astrom_offset_mean'] = exposure.metadata['SFM_ASTROM_OFFSET_MEAN']

1663 self.metadata['astrom_offset_std'] = exposure.metadata['SFM_ASTROM_OFFSET_STD']

1664 self.metadata['astrom_offset_median'] = exposure.metadata['SFM_ASTROM_OFFSET_MEDIAN']

1665 else:

1666 self.metadata['initial_to_final_wcs'] = float("nan")

1667 self.metadata['astrom_offset_mean'] = float("nan")

1668 self.metadata['astrom_offset_std'] = float("nan")

1669 self.metadata['astrom_offset_median'] = float("nan")

1670

1671 return result.matches, result.matchMeta

1672

1673 def _fit_photometry(self, exposure, stars):

1674 """Fit a photometric model to the data and return the reference

1675 matches used in the fit, and the fitted PhotoCalib.

1676

1677 Parameters

1678 ----------

1679 exposure : `lsst.afw.image.Exposure`

1680 Exposure that is being fit, to get PSF and other metadata from.

1681 Has the fit `lsst.afw.image.PhotoCalib` attached, with pixel values

1682 unchanged.

1683 stars : `lsst.afw.table.SourceCatalog`

1684 Good stars selected for use in calibration.

1685 background : `lsst.afw.math.BackgroundList`

1686 Background that was fit to the exposure during detection of the

1687 above stars.

1688

1689 Returns

1690 -------

1691 calibrated_stars : `lsst.afw.table.SourceCatalog`

1692 Star catalog with flux/magnitude columns computed from the fitted

1693 photoCalib (instFlux columns are retained as well).

1694 matches : `list` [`lsst.afw.table.ReferenceMatch`]

1695 Reference/stars matches used in the fit.

1696 matchMeta : `lsst.daf.base.PropertyList`

1697 Metadata needed to unpersist matches, as returned by the matcher.

1698 photo_calib : `lsst.afw.image.PhotoCalib`

1699 Photometric calibration that was fit to the star catalog.

1700 """

1701 result = self.photometry.run(exposure, stars)

1702 calibrated_stars = result.photoCalib.calibrateCatalog(stars)

1703 exposure.setPhotoCalib(result.photoCalib)

1704 return calibrated_stars, result.matches, result.matchMeta, result.photoCalib

1705

1706 def _apply_photometry(self, exposure, background, background_to_photometric_ratio=None):

1707 """Apply the photometric model attached to the exposure to the

1708 exposure's pixels and an associated background model.

1709

1710 Parameters

1711 ----------

1712 exposure : `lsst.afw.image.Exposure`

1713 Exposure with the target `lsst.afw.image.PhotoCalib` attached.

1714 On return, pixel values will be calibrated and an identity

1715 photometric transform will be attached.

1716 background : `lsst.afw.math.BackgroundList`

1717 Background model to convert to nanojansky units in place.

1718 background_to_photometric_ratio : `lsst.afw.image.Image`, optional

1719 Image to convert photometric-flattened image to

1720 background-flattened image.

1721 """

1722 photo_calib = exposure.getPhotoCalib()

1723 exposure.maskedImage = photo_calib.calibrateImage(exposure.maskedImage)

1724 identity = afwImage.PhotoCalib(1.0,

1725 photo_calib.getCalibrationErr(),

1726 bbox=exposure.getBBox())

1727 exposure.setPhotoCalib(identity)

1728 exposure.metadata["BUNIT"] = "nJy"

1729

1730 assert photo_calib._isConstant, \

1731 "Background calibration assumes a constant PhotoCalib; PhotoCalTask should always return that."

1732

1733 for bg in background:

1734 # The statsImage is a view, but we can't assign to a function call

1735 # in python.

1736 binned_image = bg[0].getStatsImage()

1737 binned_image *= photo_calib.getCalibrationMean()

1738

1739 def _summarize(self, exposure, stars, background):

1740 """Compute summary statistics on the exposure and update in-place the

1741 calibrations attached to it.

1742

1743 Parameters

1744 ----------

1745 exposure : `lsst.afw.image.Exposure`

1746 Exposure that was calibrated, to get PSF and other metadata from.

1747 Should be in instrumental units with the photometric calibration

1748 attached.

1749 Modified to contain the computed summary statistics.

1750 stars : `SourceCatalog`

1751 Good stars selected used in calibration.

1752 background : `lsst.afw.math.BackgroundList`

1753 Background that was fit to the exposure during detection of the

1754 above stars. Should be in instrumental units.

1755 """

1756 summary = self.compute_summary_stats.run(exposure, stars, background)

1757 exposure.info.setSummaryStats(summary)

1758

1759 def _recordMaskedPixelFractions(self, exposure):

1760 """Record the fraction of all the pixels in an exposure

1761 that are masked with a given flag. Each fraction is

1762 recorded in the task metadata. One record per flag type.

1763

1764 Parameters

1765 ----------

1766 exposure : `lsst.afw.image.ExposureF`

1767 The target exposure to calculate masked pixel fractions for.

1768 """

1769

1770 mask = exposure.mask

1771 maskPlanes = list(mask.getMaskPlaneDict().keys())

1772 for maskPlane in maskPlanes:

1773 self.metadata[f"{maskPlane.lower()}_mask_fraction"] = (

1774 evaluateMaskFraction(mask, maskPlane)

1775 )

1776

1777 def _update_wcs_with_camera_model(self, exposure, cameraModel, exposure_record=None):

1778 """Replace the existing WCS with one generated using the pointing

1779 and the input camera model.

1780

1781 If the current detector does not exist in the cameraModel, the pointing

1782 will still get updated, but the original distortion model will be used.

1783

1784 Parameters

1785 ----------

1786 exposure : `lsst.afw.image.ExposureF`

1787 The exposure whose WCS will be updated.

1788 cameraModel : `lsst.afw.cameraGeom.Camera`

1789 Camera to be used when constructing updated WCS.

1790 exposure_record : `lsst.daf.butler.DimensionRecord`, optional

1791 Butler metadata for the ``exposure`` dimension. Used if

1792 constructing an updated WCS instead of the boresight and

1793 orientation angle from the visit info.

1794 """

1795 if cameraModel.get(exposure.detector.getId()):

1796 self.log.info("Updating WCS with the provided camera model.")

1797 detector = cameraModel[exposure.detector.getId()]

1798 else:

1799 self.log.warning(

1800 "useButlerCamera=True, but detector %s is not available in the provided camera. The "

1801 "astrometry fit will use the initial distortion model for this detector.",

1802 exposure.detector.getId())

1803 detector = exposure.detector

1804 if exposure_record is None:

1805 boresight = exposure.visitInfo.getBoresightRaDec()

1806 orientation = exposure.visitInfo.getBoresightRotAngle()

1807 else:

1808 boresight = SpherePoint(exposure_record.tracking_ra, exposure_record.tracking_dec, degrees)

1809 orientation = exposure_record.sky_angle * degrees

1810 self.log.info(

1811 "Pointing from exposure record is %.6f deg (%.3f arcsec) from initial pointing; "

1812 "orientation difference is %.6f deg (%.3f arcsec).",

1813 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asDegrees(),

1814 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asArcseconds(),

1815 (orientation - exposure.visitInfo.getBoresightRotAngle()).asDegrees(),

1816 (orientation - exposure.visitInfo.getBoresightRotAngle()).asArcseconds(),

1817 )

1818 newWcs = createInitialSkyWcsFromBoresight(boresight, orientation, detector)

1819 exposure.setWcs(newWcs)

1820

1821 def _update_wcs_to_exposure_record(self, exposure, exposure_record):

1822 """Replace the existing WCS with the one from the butler derived

1823 exposure record pointing.

1824

1825 Parameters

1826 ----------

1827 exposure : `lsst.afw.image.ExposureF`

1828 The exposure whose WCS will be updated.

1829 exposure_record : `lsst.daf.butler.DimensionRecord`, optional

1830 Butler metadata for the ``exposure`` dimension. Used if

1831 constructing an updated WCS instead of the boresight and

1832 orientation angle from the visit info.

1833 """

1834 detector = exposure.detector

1835 boresight = SpherePoint(exposure_record.tracking_ra, exposure_record.tracking_dec, degrees)

1836 orientation = exposure_record.sky_angle * degrees

1837 self.log.info(

1838 "Pointing from exposure record is %.6f deg (%.3f arcsec) from initial pointing; "

1839 "orientation difference is %.6f deg (%.3f arcsec).",

1840 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asDegrees(),

1841 boresight.separation(exposure.visitInfo.getBoresightRaDec()).asArcseconds(),

1842 (orientation - exposure.visitInfo.getBoresightRotAngle()).asDegrees(),

1843 (orientation - exposure.visitInfo.getBoresightRotAngle()).asArcseconds(),

1844 )

1845 newWcs = createInitialSkyWcsFromBoresight(boresight, orientation, detector)

1846 exposure.setWcs(newWcs)

1847

1848 def _compute_mask_fraction(self, mask, mask_planes, bad_mask_planes):

1849 """Evaluate the fraction of masked pixels in a (set of) mask plane(s).

1850

1851 Parameters

1852 ----------

1853 mask : `lsst.afw.image.Mask`

1854 The mask on which to evaluate the fraction.

1855 mask_planes : `list`, `str`

1856 The mask planes for which to evaluate the fraction.

1857 bad_mask_planes : `list`, `str`

1858 The mask planes to exclude from the computation.

1859

1860 Returns

1861 -------

1862 detected_fraction : `float`

1863 The calculated fraction of masked pixels

1864 """

1865 bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes)

1866 n_good_pix = np.sum(mask.array & bad_pixel_mask == 0)

1867 if n_good_pix == 0:

1868 detected_fraction = float("nan")

1869 return detected_fraction

1870 detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes)

1871 n_detected_pix = np.sum((mask.array & detected_pixel_mask != 0)

1872 & (mask.array & bad_pixel_mask == 0))

1873 detected_fraction = n_detected_pix/n_good_pix

1874 return detected_fraction

1875

1876 def _compute_per_amp_fraction(self, exposure, detected_fraction, mask_planes, bad_mask_planes):

1877 """Evaluate the maximum per-amplifier fraction of masked pixels.

1878

1879 Parameters

1880 ----------

1881 exposure : `lsst.afw.image.ExposureF`

1882 The exposure on which to compute the per-amp masked fraction.

1883 detected_fraction : `float`

1884 The current detected_fraction of the ``mask_planes`` for the

1885 full detector.

1886 mask_planes : `list`, `str`

1887 The mask planes for which to evaluate the fraction.

1888 bad_mask_planes : `list`, `str`

1889 The mask planes to exclude from the computation.

1890

1891 Returns

1892 -------

1893 n_above_max_per_amp : `int`

1894 The number of amplifiers with masked fractions above a maximum

1895 value (set by the current full-detector ``detected_fraction``).

1896 highest_detected_fraction_per_amp : `float`

1897 The highest value of the per-amplifier fraction of masked pixels.

1898 no_zero_det_amps : `bool`

1899 A boolean representing whether any of the amplifiers has zero

1900 masked pixels.

1901 """

1902 highest_detected_fraction_per_amp = -9.99

1903 n_above_max_per_amp = 0

1904 n_no_zero_det_amps = 0

1905 no_zero_det_amps = True

1906 amps = exposure.detector.getAmplifiers()

1907 if amps is not None:

1908 for ia, amp in enumerate(amps):

1909 amp_bbox = amp.getBBox()

1910 exp_bbox = exposure.getBBox()

1911 if not exp_bbox.contains(amp_bbox):

1912 self.log.info("Bounding box of amplifier (%s) does not fit in exposure's "

1913 "bounding box (%s). Skipping...", amp_bbox, exp_bbox)

1914 continue

1915 sub_image = exposure.subset(amp.getBBox())

1916 detected_fraction_amp = self._compute_mask_fraction(sub_image.mask,

1917 mask_planes,

1918 bad_mask_planes)

1919 self.log.debug("Current detected fraction for amplifier %s = %.3f",

1920 amp.getName(), detected_fraction_amp)

1921 if detected_fraction_amp < 0.002:

1922 n_no_zero_det_amps += 1

1923 if n_no_zero_det_amps > 2:

1924 no_zero_det_amps = False

1925 break

1926 highest_detected_fraction_per_amp = max(detected_fraction_amp,

1927 highest_detected_fraction_per_amp)

1928 if highest_detected_fraction_per_amp > min(0.998, max(0.8, 3.0*detected_fraction)):

1929 n_above_max_per_amp += 1

1930 if n_above_max_per_amp > 2:

1931 break

1932 else:

1933 self.log.info("No amplifier object for detector %d, so skipping per-amp "

1934 "detection fraction checks.", exposure.detector.getId())

1935 return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps

1936

1937 def _remeasure_star_background(self, result, background_to_photometric_ratio=None):

1938 """Remeasure the exposure's background with iterative adaptive

1939 threshold detection.

1940

1941 Parameters

1942 ----------

1943 result : `lsst.pipe.base.Struct`

1944 The current state of the result Strut from the run method of

1945 CalibrateImageTask. Will be modified in place.

1946 background_to_photometric_ratio : `lsst.afw.image.Image`, optional

1947 Image to convert photometric-flattened image to

1948 background-flattened image.

1949

1950 Returns

1951 -------

1952 result : `lsst.pipe.base.Struct`

1953 The modified result Struct with the new background subtracted.

1954 """

1955 # Restore the previously measured background and remeasure it

1956 # using an adaptive threshold detection iteration to ensure a

1957 # "Goldilocks Zone" for the fraction of detected pixels.

1958 backgroundOrig = result.background.clone()

1959 median_background_orig = np.median(backgroundOrig.getImage().array)

1960 self.log.info("Original median_background = %.2f", median_background_orig)

1961 result.exposure.image.array += result.background.getImage().array

1962 result.background = afwMath.BackgroundList()

1963

1964 origMask = result.exposure.mask.clone()

1965 bad_mask_planes = ["BAD", "EDGE", "NO_DATA"]

1966 detected_mask_planes = ["DETECTED", "DETECTED_NEGATIVE"]

1967 detected_fraction_orig = self._compute_mask_fraction(result.exposure.mask,

1968 detected_mask_planes,

1969 bad_mask_planes)

1970 minDetFracForFinalBg = 0.02

1971 maxDetFracForFinalBg = 0.93

1972 # Dilate the current detected mask planes and don't clear

1973 # them in the detection step.

1974 nPixToDilate = 10

1975 maxAdaptiveDetIter = 10

1976 for dilateIter in range(maxAdaptiveDetIter):

1977 dilatedMask = origMask.clone()

1978 for maskName in detected_mask_planes:

1979 # Compute the grown detection mask plane using SpanSet

1980 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)

1981 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)

1982 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)

1983 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())

1984 # Clear the detected mask plane

1985 detectedMask = dilatedMask.getMaskPlane(maskName)

1986 dilatedMask.clearMaskPlane(detectedMask)

1987 # Set the mask plane to the dilated one

1988 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)

1989

1990 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,

1991 detected_mask_planes,

1992 bad_mask_planes)

1993 if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1:

1994 break

1995 else:

1996 nPixToDilate -= 1

1997

1998 result.exposure.mask = dilatedMask

1999 self.log.debug("detected_fraction_orig = %.3f; detected_fraction_dilated = %.3f",

2000 detected_fraction_orig, detected_fraction_dilated)

2001 n_above_max_per_amp = -99

2002 highest_detected_fraction_per_amp = float("nan")

2003

2004 # Check the per-amplifier detection fractions.

2005 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \

2006 self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated,

2007 detected_mask_planes, bad_mask_planes)

2008 self.log.debug("Dilated mask: n_above_max_per_amp = %d; "

2009 "highest_detected_fraction_per_amp = %.3f",

2010 n_above_max_per_amp, highest_detected_fraction_per_amp)

2011

2012 bgIgnoreMasksToAdd = ["SAT", "SUSPECT", "SPIKE"]

2013 detected_fraction = 1.0

2014 nFootprintTemp = 1e12

2015 starBackgroundDetectionConfig = lsst.meas.algorithms.SourceDetectionConfig()

2016 for maskName in bgIgnoreMasksToAdd:

2017 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()

2018 and maskName not in starBackgroundDetectionConfig.background.ignoredPixelMask):

2019 starBackgroundDetectionConfig.background.ignoredPixelMask += [maskName]

2020 starBackgroundDetectionConfig.doTempLocalBackground = False

2021 starBackgroundDetectionConfig.nSigmaToGrow = 70.0

2022 starBackgroundDetectionConfig.reEstimateBackground = False

2023 starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0

2024 starBackgroundDetectionConfig.thresholdValue = max(2.0, 0.2*median_background_orig)

2025 starBackgroundDetectionConfig.thresholdType = "pixel_stdev"

2026

2027 n_above_max_per_amp = -99

2028 highest_detected_fraction_per_amp = float("nan")

2029 doCheckPerAmpDetFraction = True

2030

2031 minFootprints = self.config.star_background_min_footprints

2032 maxIter = 40

2033 for nIter in range(maxIter):

2034 currentThresh = starBackgroundDetectionConfig.thresholdValue

2035 nZeroEncountered = 0

2036 if nFootprintTemp == 0:

2037 zeroFactor = min(0.98, 0.9 + 0.01*nZeroEncountered)

2038 starBackgroundDetectionConfig.thresholdValue = zeroFactor*currentThresh

2039 self.log.warning("No footprints detected. Decreasing threshold to %.2f and rerunning.",

2040 starBackgroundDetectionConfig.thresholdValue)

2041 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(

2042 config=starBackgroundDetectionConfig)

2043 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)

2044 tempDetections = starBackgroundDetectionTask.run(

2045 table=table, exposure=result.exposure, clearMask=True)

2046 nFootprintTemp = len(tempDetections.sources)

2047 minFootprints = max(self.config.star_background_min_footprints,

2048 int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))

2049 minFootprints = min(200, minFootprints)

2050 nZeroEncountered += 1

2051 if nFootprintTemp >= minFootprints:

2052 detected_fraction = self._compute_mask_fraction(result.exposure.mask,

2053 detected_mask_planes,

2054 bad_mask_planes)

2055 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "

2056 "DETECTED_NEGATIVE in star_background_detection = %.3f "

2057 "(max is %.3f; min is %.3f) nFootprint = %d (current min is %d)",

2058 nIter, starBackgroundDetectionConfig.thresholdValue,

2059 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,

2060 nFootprintTemp, minFootprints)

2061 self.metadata['adaptive_threshold_value'] = starBackgroundDetectionConfig.thresholdValue

2062

2063 break

2064 else:

2065 # Still not enough footprints, so make sure this loop is

2066 # entered again.

2067 if nFootprintTemp > 0 and nFootprintTemp < minFootprints:

2068 nFootprintTemp = 0

2069 continue

2070 if detected_fraction > maxDetFracForFinalBg or nFootprintTemp <= minFootprints:

2071 starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh

2072 if nFootprintTemp < minFootprints and detected_fraction > 0.9*maxDetFracForFinalBg:

2073 if nFootprintTemp == 1:

2074 starBackgroundDetectionConfig.thresholdValue = 1.4*currentThresh

2075 else:

2076 starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh

2077

2078 if n_above_max_per_amp > 1:

2079 starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh

2080 if detected_fraction < minDetFracForFinalBg:

2081 starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh

2082 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(

2083 config=starBackgroundDetectionConfig)

2084 table = afwTable.SourceTable.make(self.initial_stars_schema.schema)

2085 tempDetections = starBackgroundDetectionTask.run(

2086 table=table, exposure=result.exposure, clearMask=True)

2087 result.exposure.mask |= dilatedMask

2088 nFootprintTemp = len(tempDetections.sources)

2089 minFootprints = max(self.config.star_background_min_footprints,

2090 int(self.config.star_background_peak_fraction*tempDetections.numPosPeaks))

2091 minFootprints = min(200, minFootprints)

2092 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,

2093 bad_mask_planes)

2094 self.log.info("nIter = %d, thresh = %.2f: Fraction of pixels marked as DETECTED or "

2095 "DETECTED_NEGATIVE in star_background_detection = %.3f "

2096 "(max is %.3f; min is %.3f); nFootprint = %d (current min is %d)",

2097 nIter, starBackgroundDetectionConfig.thresholdValue,

2098 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,

2099 nFootprintTemp, minFootprints)

2100 self.metadata['adaptive_threshold_value'] = starBackgroundDetectionConfig.thresholdValue

2101

2102 n_amp = len(result.exposure.detector.getAmplifiers())

2103 if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg:

2104 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \

2105 self._compute_per_amp_fraction(result.exposure, detected_fraction,

2106 detected_mask_planes, bad_mask_planes)

2107

2108 if not no_zero_det_amps:

2109 starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh

2110

2111 if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg

2112 and n_above_max_per_amp < int(0.75*n_amp)

2113 and no_zero_det_amps

2114 and nFootprintTemp >= minFootprints):

2115 if (n_above_max_per_amp < max(1, int(0.15*n_amp))

2116 or detected_fraction < 0.85*maxDetFracForFinalBg):

2117 break

2118 else:

2119 starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh

2120 self.log.debug("Number of amplifiers with detected fraction above the maximum "

2121 "(%.2f) excedes the maximum allowed (%d >= %d). Making a small "

2122 "tweak to the threshold (from %.2f to %.2f) and rerunning.",

2123 maxDetFracForFinalBg, n_above_max_per_amp, int(0.75*n_amp),

2124 currentThresh, starBackgroundDetectionConfig.thresholdValue)

2125 self.log.debug("n_above_max_per_amp = %d (abs max is %d)", n_above_max_per_amp, int(0.75*n_amp))

2126

2127 detected_fraction = self._compute_mask_fraction(result.exposure.mask, detected_mask_planes,

2128 bad_mask_planes)

2129 self.log.info("Fraction of pixels marked as DETECTED or DETECTED_NEGATIVE is now %.5f "

2130 "(highest per amp section = %.5f)",

2131 detected_fraction, highest_detected_fraction_per_amp)

2132

2133 if detected_fraction > maxDetFracForFinalBg:

2134 result.exposure.mask = dilatedMask

2135 self.log.warning("Final fraction of pixels marked as DETECTED or DETECTED_NEGATIVE "

2136 "was too large in star_background_detection = %.3f (max = %.3f). "

2137 "Reverting to dilated mask from PSF detection.",

2138 detected_fraction, maxDetFracForFinalBg)

2139 self.metadata['adaptive_threshold_value'] = float("nan")

2140 star_background = self.star_background.run(

2141 exposure=result.exposure,

2142 backgroundToPhotometricRatio=background_to_photometric_ratio,

2143 ).background

2144 result.background.append(star_background[0])

2145

2146 median_background = np.median(result.background.getImage().array)

2147 self.log.info("New initial median_background = %.2f", median_background)

2148

2149 # Perform one more round of background subtraction that is just an

2150 # overall pedestal (order = 0). This is intended to account for

2151 # any potential gross oversubtraction imposed by the higher-order

2152 # subtraction.

2153 # Dilate DETECTED mask a bit more if it's below 50% detected.

2154 nPixToDilate = 2

2155 if detected_fraction < 0.5:

2156 dilatedMask = result.exposure.mask.clone()

2157 for maskName in detected_mask_planes:

2158 # Compute the grown detection mask plane using SpanSet

2159 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)

2160 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)

2161 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)

2162 detectedMaskSpanSet = detectedMaskSpanSet.clippedTo(dilatedMask.getBBox())

2163 # Clear the detected mask plane

2164 detectedMask = dilatedMask.getMaskPlane(maskName)

2165 dilatedMask.clearMaskPlane(detectedMask)

2166 # Set the mask plane to the dilated one

2167 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)

2168

2169 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,

2170 detected_mask_planes,

2171 bad_mask_planes)

2172 result.exposure.mask = dilatedMask

2173 self.log.debug("detected_fraction_orig = %.3f; detected_fraction_dilated = %.3f",

2174 detected_fraction_orig, detected_fraction_dilated)

2175

2176 bbox = result.exposure.getBBox()

2177

2178 # Now measure a final background pedestal level (largely to account

2179 # for the over-subtraction often seen in the first pass, especially

2180 # in high fill-factor fields). Since the pedestal measurement can be

2181 # sensitive to the bin size, start an iteration with a small bin size,

2182 # then double it on each iteration until the relative or absolute

2183 # change criterion is met. If those are never achieved, the iteration

2184 # stops when the bin size gets bigger than the exposure's bounding box.

2185

2186 # Initialize the parameters for the pedestal iteration.

2187 pedestalBackgroundConfig = lsst.meas.algorithms.SubtractBackgroundConfig()

2188 for maskName in bgIgnoreMasksToAdd:

2189 if (maskName in result.exposure.mask.getMaskPlaneDict().keys()

2190 and maskName not in pedestalBackgroundConfig.ignoredPixelMask):

2191 pedestalBackgroundConfig.ignoredPixelMask += [maskName]

2192 pedestalBackgroundConfig.statisticsProperty = "MEDIAN"

2193 pedestalBackgroundConfig.approxOrderX = 0

2194 pedestalBackgroundConfig.binSize = min(32, int(math.ceil(max(bbox.width, bbox.height)/4.0)))

2195 self.log.info("Initial pedestal binSize = %d pixels", pedestalBackgroundConfig.binSize)

2196 cumulativePedestalLevel = 0.0

2197 # When the cumulative pedestal value changes by less than 5% from one

2198 # bin size to the next we assume we are converged.

2199 relativeStoppingCriterion = 0.05

2200 # When the cumulative pedestal value changes by less than 0.5 counts

2201 # (electrons or adu) from one bin size to the next, we assume we are

2202 # converged.

2203 absoluteStoppingCriterion = 0.5

2204 inPedestalIteration = True

2205 while inPedestalIteration:

2206 cumulativePedestalPrev = cumulativePedestalLevel

2207 pedestalBackgroundTask = lsst.meas.algorithms.SubtractBackgroundTask(

2208 config=pedestalBackgroundConfig)

2209 pedestalBackgroundList = pedestalBackgroundTask.run(

2210 exposure=result.exposure,

2211 background=result.background,

2212 backgroundToPhotometricRatio=background_to_photometric_ratio,

2213 ).background

2214 # Isolate the most recent pedestal background measured to log the

2215 # current computed value and add it to the cumulative value.

2216 pedestalBackground = afwMath.BackgroundList()

2217 pedestalBackground.append(pedestalBackgroundList[-1])

2218 pedestalBackgroundLevel = pedestalBackground.getImage().array[0, 0]

2219 cumulativePedestalLevel += pedestalBackgroundLevel

2220 absoluteDelta = abs(cumulativePedestalLevel - cumulativePedestalPrev)

2221 if cumulativePedestalPrev != 0.0:

2222 relativeDelta = abs(absoluteDelta/cumulativePedestalPrev)

2223 else:

2224 relativeDelta = 1.0

2225 self.log.info("Subtracted pedestalBackgroundLevel = %.4f (cumulative = %.4f; "

2226 "relativeDelta = %.4f; absoluteDelta = %.3f)",

2227 pedestalBackgroundLevel, cumulativePedestalLevel,

2228 relativeDelta, absoluteDelta)

2229

2230 if (relativeDelta < relativeStoppingCriterion

2231 or absoluteDelta < absoluteStoppingCriterion):

2232 # We are converged.

2233 inPedestalIteration = False

2234 else:

2235 # We have not yet converged; grow the bin size if possible.

2236 if pedestalBackgroundConfig.binSize*2 < 2*max(bbox.width, bbox.height):

2237 pedestalBackgroundConfig.binSize = int(pedestalBackgroundConfig.binSize*2)

2238 self.log.info("Increasing pedestal binSize to %d pixels and remeasuring.",

2239 pedestalBackgroundConfig.binSize)

2240 else:

2241 # We have reached the maximum bin size.

2242 self.log.warning("Reached maximum bin size. Exiting pedestal loop without meeting "

2243 "the convergence criteria (relativeDelta <= %.2f or "

2244 "absoluteDelta <= %.2f).",

2245 relativeStoppingCriterion, absoluteStoppingCriterion)

2246 inPedestalIteration = False

2247 self.log.info("Final subtracted cumulativePedestalBackgroundLevel = %.4f", cumulativePedestalLevel)

2248

2249 # Clear detected mask plane before final round of detection

2250 mask = result.exposure.mask

2251 for mp in detected_mask_planes:

2252 if mp not in mask.getMaskPlaneDict():

2253 mask.addMaskPlane(mp)

2254 mask &= ~mask.getPlaneBitMask(detected_mask_planes)

2255

2256 return result

Coverage for python / lsst / pipe / tasks / calibrateImage.py: 12%

806 statements