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 do_remeasure_star_background = pexConfig.Field(

300 dtype=bool,

301 default=True,

302 doc="Do iterative star background measurement (ignored if do_adaptive_threshold_detection is False).",

303 )

304 psf_adaptive_threshold_detection = pexConfig.ConfigurableField(

305 target=AdaptiveThresholdDetectionTask,

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

307 )

308 psf_source_measurement = pexConfig.ConfigurableField(

309 target=lsst.meas.base.SingleFrameMeasurementTask,

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

311 )

312 psf_measure_psf = pexConfig.ConfigurableField(

313 target=measurePsf.MeasurePsfTask,

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

315 )

316 psf_normalized_calibration_flux = pexConfig.ConfigurableField(

317 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,

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

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

320 )

321

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

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

324 measure_aperture_correction = pexConfig.ConfigurableField(

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

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

327 )

328

329 # subtasks used during star measurement

330 star_background = pexConfig.ConfigurableField(

331 target=lsst.meas.algorithms.SubtractBackgroundTask,

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

333 )

334 star_background_min_footprints = pexConfig.Field(

335 dtype=int,

336 default=3,

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

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

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

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

341 "threshold is met.",

342 )

343 star_background_peak_fraction = pexConfig.Field(

344 dtype=float,

345 default=0.01,

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

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

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

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

350 "threshold is met.",

351 )

352 star_detection = pexConfig.ConfigurableField(

353 target=lsst.meas.algorithms.SourceDetectionTask,

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

355 )

356 star_sky_sources = pexConfig.ConfigurableField(

357 target=lsst.meas.algorithms.SkyObjectsTask,

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

359 )

360 do_downsample_footprints = pexConfig.Field(

361 dtype=bool,

362 default=False,

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

364 )

365 downsample_max_footprints = pexConfig.Field(

366 dtype=int,

367 default=1000,

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

369 )

370 star_deblend = pexConfig.ConfigurableField(

371 target=lsst.meas.deblender.SourceDeblendTask,

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

373 )

374 star_measurement = pexConfig.ConfigurableField(

375 target=lsst.meas.base.SingleFrameMeasurementTask,

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

377 )

378 star_normalized_calibration_flux = pexConfig.ConfigurableField(

379 target=lsst.meas.algorithms.NormalizedCalibrationFluxTask,

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

381 "for the final output star catalog.",

382 )

383 star_apply_aperture_correction = pexConfig.ConfigurableField(

384 target=lsst.meas.base.ApplyApCorrTask,

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

386 )

387 star_catalog_calculation = pexConfig.ConfigurableField(

388 target=lsst.meas.base.CatalogCalculationTask,

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

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

391 )

392 star_set_primary_flags = pexConfig.ConfigurableField(

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

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

395 )

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

397 default="science",

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

399 )

400

401 # final calibrations and statistics

402 astrometry = pexConfig.ConfigurableField(

403 target=lsst.meas.astrom.AstrometryTask,

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

405 )

406 astrometry_ref_loader = pexConfig.ConfigField(

407 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,

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

409 )

410 photometry = pexConfig.ConfigurableField(

411 target=photoCal.PhotoCalTask,

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

413 )

414 photometry_ref_loader = pexConfig.ConfigField(

415 dtype=lsst.meas.algorithms.LoadReferenceObjectsConfig,

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

417 )

418 doMaskDiffractionSpikes = pexConfig.Field(

419 dtype=bool,

420 default=False,

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

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

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

424 )

425 diffractionSpikeMask = pexConfig.ConfigurableField(

426 target=diffractionSpikeMask.DiffractionSpikeMaskTask,

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

428 )

429

430 compute_summary_stats = pexConfig.ConfigurableField(

431 target=computeExposureSummaryStats.ComputeExposureSummaryStatsTask,

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

433 )

434

435 do_illumination_correction = pexConfig.Field(

436 dtype=bool,

437 default=False,

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

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

440 "for background subtraction.",

441 )

442

443 do_calibrate_pixels = pexConfig.Field(

444 dtype=bool,

445 default=True,

446 doc=(

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

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

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

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

451 "them to physical units."

452 )

453 )

454

455 do_include_astrometric_errors = pexConfig.Field(

456 dtype=bool,

457 default=True,

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

459 )

460

461 background_stats_ignored_pixel_masks = pexConfig.ListField(

462 dtype=str,

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

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

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

466 "meas.algorithms.SubtractBackgroundConfig algorithm."

467 )

468

469 run_sattle = pexConfig.Field(

470 dtype=bool,

471 default=False,

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

473 "in ip_diffim.detectAndMeasure alert verification."

474 )

475

476 sattle_historical = pexConfig.Field(

477 dtype=bool,

478 default=False,

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

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

481 "closest in time to the exposure.",

482 )

483

484 useButlerExposureRegion = pexConfig.Field(

485 dtype=bool,

486 default=False,

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

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

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

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

491 )

492

493 useButlerCamera = pexConfig.Field(

494 dtype=bool,

495 default=False,

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

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

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

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

500 "camera model from the obs_* package."

501 )

502

503 background_stats_flux_column = pexConfig.Field(

504 dtype=str,

505 default="base_CircularApertureFlux_12_0_flux",

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

507 )

508

509 def setDefaults(self):

510 super().setDefaults()

511

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

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

514 # like CRs for very good seeing images.

515 self.install_simple_psf.fwhm = 4

516

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

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

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

520 # thresholdValue * includeThresholdMultiplier. The low thresholdValue

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

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

523 self.psf_detection.thresholdValue = 10.0

524 self.psf_detection.includeThresholdMultiplier = 5.0

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

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

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

528 self.psf_detection.doTempLocalBackground = False

529 self.psf_detection.reEstimateBackground = False

530

531 # Minimal measurement plugins for PSF determination.

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

533 # shapeHSM/moments for star/galaxy separation.

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

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

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

537 "base_SdssCentroid",

538 "ext_shapeHSM_HsmSourceMoments",

539 "base_CircularApertureFlux",

540 "base_GaussianFlux",

541 "base_PsfFlux",

542 "base_CompensatedTophatFlux",

543 ]

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

545 # Only measure apertures we need for PSF measurement.

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

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

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

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

550 "base_CircularApertureFlux_12_0_instFlux"

551

552 # No extendedness information available: we need the aperture

553 # corrections to determine that.

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

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

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

557

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

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

560 self.star_detection.reEstimateBackground = False

561 self.star_detection.thresholdValue = 5.0

562 self.star_detection.includeThresholdMultiplier = 2.0

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

564 "base_SdssCentroid",

565 "ext_shapeHSM_HsmSourceMoments",

566 "ext_shapeHSM_HsmPsfMoments",

567 "base_GaussianFlux",

568 "base_PsfFlux",

569 "base_CircularApertureFlux",

570 "base_ClassificationSizeExtendedness",

571 "base_CompensatedTophatFlux",

572 ]

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

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

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

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

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

578

579 # We measure and apply the normalization aperture correction with

580 # the psf_normalized_calibration_flux task, and we only apply the

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

582 self.star_normalized_calibration_flux.do_measure_ap_corr = False

583

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

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

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

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

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

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

590 # wanted for calibration.

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

592 # Set the flux and error fields

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

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

595

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

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

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

599 self.astrometry_ref_loader.anyFilterMapsToThis = "phot_g_mean"

600

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

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

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

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

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

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

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

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

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

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

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

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

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

614

615 def validate(self):

616 super().validate()

617

618 # Ensure that the normalization calibration flux tasks

619 # are configured correctly.

620 if not self.psf_normalized_calibration_flux.do_measure_ap_corr:

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

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

623 raise pexConfig.FieldValidationError(

624 CalibrateImageConfig.psf_normalized_calibration_flux, self, msg,

625 )

626 if self.star_normalized_calibration_flux.do_measure_ap_corr:

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

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

629 "fluxes.")

630 raise pexConfig.FieldValidationError(

631 CalibrateImageConfig.star_normalized_calibration_flux, self, msg,

632 )

633

634 # Ensure base_LocalPhotoCalib and base_LocalWcs plugins are not run,

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

636 # and WCS.

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

638 raise pexConfig.FieldValidationError(

639 CalibrateImageConfig.psf_source_measurement,

640 self,

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

642 )

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

644 raise pexConfig.FieldValidationError(

645 CalibrateImageConfig.star_measurement,

646 self,

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

648 )

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

650 raise pexConfig.FieldValidationError(

651 CalibrateImageConfig.psf_source_measurement,

652 self,

653 "base_LocalPhotoCalib cannot run CalibrateImageTask, "

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

655 )

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

657 raise pexConfig.FieldValidationError(

658 CalibrateImageConfig.star_measurement,

659 self,

660 "base_LocalPhotoCalib cannot run CalibrateImageTask, "

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

662 )

663

664 # Check for illumination correction and background consistency.

665 if self.do_illumination_correction:

666 if not self.psf_subtract_background.doApplyFlatBackgroundRatio:

667 raise pexConfig.FieldValidationError(

668 CalibrateImageConfig.psf_subtract_background,

669 self,

670 "CalibrateImageTask.psf_subtract_background must be configured with "

671 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

672 )

673 if self.psf_detection.reEstimateBackground:

674 if not self.psf_detection.doApplyFlatBackgroundRatio:

675 raise pexConfig.FieldValidationError(

676 CalibrateImageConfig.psf_detection,

677 self,

678 "CalibrateImageTask.psf_detection background must be configured with "

679 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

680 )

681 if self.star_detection.reEstimateBackground:

682 if not self.star_detection.doApplyFlatBackgroundRatio:

683 raise pexConfig.FieldValidationError(

684 CalibrateImageConfig.star_detection,

685 self,

686 "CalibrateImageTask.star_detection background must be configured with "

687 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

688 )

689 if not self.star_background.doApplyFlatBackgroundRatio:

690 raise pexConfig.FieldValidationError(

691 CalibrateImageConfig.star_background,

692 self,

693 "CalibrateImageTask.star_background background must be configured with "

694 "doApplyFlatBackgroundRatio=True if do_illumination_correction=True."

695 )

696

697 if self.run_sattle:

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

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

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

701

702 if not self.do_adaptive_threshold_detection:

703 if not self.psf_detection.reEstimateBackground:

704 raise pexConfig.FieldValidationError(

705 CalibrateImageConfig.psf_detection,

706 self,

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

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

709 "background must be configured with "

710 "reEstimateBackground = True to maintain original behavior."

711 )

712 if not self.star_detection.reEstimateBackground:

713 raise pexConfig.FieldValidationError(

714 CalibrateImageConfig.psf_detection,

715 self,

716 "If not using the adaptive threshold detection implementation "

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

718 "CalibrateImageTask.star_detection background must be configured "

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

720 )

721

722

723class CalibrateImageTask(pipeBase.PipelineTask):

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

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

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

727

728 Parameters

729 ----------

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

731 Schema of the initial_stars output catalog.

732 """

733 _DefaultName = "calibrateImage"

734 ConfigClass = CalibrateImageConfig

735

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

737 super().__init__(**kwargs)

738 self.makeSubtask("snap_combine")

739

740 # PSF determination subtasks

741 self.makeSubtask("install_simple_psf")

742 self.makeSubtask("psf_repair")

743 self.makeSubtask("psf_subtract_background")

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

745 self.psf_schema.addField(

746 'psf_max_value',

747 type=np.float32,

748 doc="PSF max value.",

749 )

750 afwTable.CoordKey.addErrorFields(self.psf_schema)

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

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

753 self.makeSubtask("psf_adaptive_threshold_detection")

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

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

756

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

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

759

760 # star measurement subtasks

761 if initial_stars_schema is None:

762 initial_stars_schema = afwTable.SourceTable.makeMinimalSchema()

763

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

765 # astrometric fitting, and aperture correction calculations.

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

767 "calib_astrometry_used",

768 # TODO DM-39203:

769 # these can be removed once apcorr is gone.

770 "apcorr_slot_CalibFlux_used", "apcorr_base_GaussianFlux_used",

771 "apcorr_base_PsfFlux_used",)

772 for field in self.psf_fields:

773 item = self.psf_schema.find(field)

774 initial_stars_schema.addField(item.getField())

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

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

777 initial_stars_schema.addField("psf_id",

778 type=id_type,

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

780 initial_stars_schema.addField("psf_max_value",

781 type=psf_max_value_type,

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

783

784 afwTable.CoordKey.addErrorFields(initial_stars_schema)

785 self.makeSubtask("star_background")

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

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

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

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

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

791

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

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

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

795 self.makeSubtask("star_selector")

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

797 if self.config.doMaskDiffractionSpikes:

798 self.makeSubtask("diffractionSpikeMask")

799 self.makeSubtask("compute_summary_stats")

800

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

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

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

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

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

806 # - note that calibrateCatalog will happily reuse existing output

807 # columns.

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

809 self.initial_stars_schema = dummy_photo_calib.calibrateCatalog(

810 afwTable.SourceCatalog(initial_stars_schema)

811 )

812

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

814 inputs = butlerQC.get(inputRefs)

815 exposures = inputs.pop("exposures")

816

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

818 if self.config.useButlerExposureRegion:

819 exposure_region = butlerQC.quantum.dataId.region

820 else:

821 exposure_region = None

822

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

824

825 astrometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(

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

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

828 name=self.config.connections.astrometry_ref_cat,

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

830 self.astrometry.setRefObjLoader(astrometry_loader)

831

832 photometry_loader = lsst.meas.algorithms.ReferenceObjectLoader(

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

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

835 name=self.config.connections.photometry_ref_cat,

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

837 self.photometry.match.setRefObjLoader(photometry_loader)

838

839 if self.config.doMaskDiffractionSpikes:

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

841 # mask.

842 self.diffractionSpikeMask.setRefObjLoader(photometry_loader)

843

844 if self.config.do_illumination_correction:

845 background_flat = inputs.pop("background_flat")

846 illumination_correction = inputs.pop("illumination_correction")

847 else:

848 background_flat = None

849 illumination_correction = None

850

851 camera_model = None

852 if self.config.useButlerCamera:

853 if "camera_model" in inputs:

854 camera_model = inputs.pop("camera_model")

855 else:

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

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

858 exposures[0].filter.bandLabel)

859

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

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

862

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

864 # exist, even as None.

865 result = pipeBase.Struct(

866 exposure=None,

867 stars_footprints=None,

868 psf_stars_footprints=None,

869 background_to_photometric_ratio=None,

870 )

871 try:

872 self.run(

873 exposures=exposures,

874 result=result,

875 id_generator=id_generator,

876 background_flat=background_flat,

877 illumination_correction=illumination_correction,

878 camera_model=camera_model,

879 exposure_record=exposure_record,

880 exposure_region=exposure_region,

881 )

882 except pipeBase.AlgorithmError as e:

883 error = pipeBase.AnnotatedPartialOutputsError.annotate(

884 e,

885 self,

886 result.exposure,

887 result.psf_stars_footprints,

888 result.stars_footprints,

889 log=self.log

890 )

891 butlerQC.put(result, outputRefs)

892 raise error from e

893

894 butlerQC.put(result, outputRefs)

895

896 @timeMethod

897 def run(

898 self,

899 *,

900 exposures,

901 id_generator=None,

902 result=None,

903 background_flat=None,

904 illumination_correction=None,

905 camera_model=None,

906 exposure_record=None,

907 exposure_region=None,

908 ):

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

910 and measurement and calibrate astrometry and photometry from that.

911

912 Parameters

913 ----------

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

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

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

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

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

919 before doing further processing.

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

921 Object that generates source IDs and provides random seeds.

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

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

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

925 this is also returned.

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

927 Background flat-field image.

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

929 Illumination correction image.

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

931 Camera to be used if constructing updated WCS.

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

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

934 constructing an updated WCS instead of the boresight and

935 orientation angle from the visit info.

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

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

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

939 of the exposure.

940

941 Returns

942 -------

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

944 Results as a struct with attributes:

945

946 ``exposure``

947 Calibrated exposure, with pixels in nJy units.

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

949 ``stars``

950 Stars that were used to calibrate the exposure, with

951 calibrated fluxes and magnitudes.

952 (`astropy.table.Table`)

953 ``stars_footprints``

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

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

956 ``psf_stars``

957 Stars that were used to determine the image PSF.

958 (`astropy.table.Table`)

959 ``psf_stars_footprints``

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

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

962 ``background``

963 Background that was fit to the exposure when detecting

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

965 ``applied_photo_calib``

966 Photometric calibration that was fit to the star catalog and

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

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

969 ``astrometry_matches``

970 Reference catalog stars matches used in the astrometric fit.

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

972 `lsst.afw.table.BaseCatalog`).

973 ``photometry_matches``

974 Reference catalog stars matches used in the photometric fit.

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

976 `lsst.afw.table.BaseCatalog`).

977 ``mask``

978 Copy of the mask plane of `exposure`.

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

980 """

981 if result is None:

982 result = pipeBase.Struct()

983 if id_generator is None:

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

985

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

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

988

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

990

991 # Check input image processing.

992 if self.config.do_illumination_correction:

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

994 raise pipeBase.InvalidQuantumError(

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

996 )

997

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

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

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

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

1002 if camera_model:

1003 self._update_wcs_with_camera_model(

1004 result.exposure, camera_model, exposure_record=exposure_record

1005 )

1006 elif exposure_record is not None:

1007 self._update_wcs_to_exposure_record(result.exposure, exposure_record)

1008

1009 result.background = None

1010 result.background_to_photometric_ratio = None

1011 summary_stat_catalog = None

1012 # Some exposure components are set to initial placeholder objects

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

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

1015 # masquerade as the real thing.

1016 have_fit_psf = False

1017 have_fit_astrometry = False

1018 have_fit_photometry = False

1019 try:

1020 result.background_to_photometric_ratio = self._apply_illumination_correction(

1021 result.exposure,

1022 background_flat,

1023 illumination_correction,

1024 )

1025

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

1027 result,

1028 id_generator,

1029 )

1030 have_fit_psf = True

1031

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

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

1034 # model is informative.

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

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

1037 raise AllCentroidsFlaggedError(

1038 n_sources=len(result.psf_stars_footprints),

1039 psf_shape_ixx=psf_shape.getIxx(),

1040 psf_shape_iyy=psf_shape.getIyy(),

1041 psf_shape_ixy=psf_shape.getIxy(),

1042 psf_size=psf_shape.getDeterminantRadius(),

1043 )

1044

1045 if result.background is None:

1046 result.background = afwMath.BackgroundList()

1047

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

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

1050 # Run astrometry using PSF candidate stars.

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

1052 afwTable.updateSourceCoords(

1053 result.exposure.wcs,

1054 sourceList=result.psf_stars_footprints,

1055 include_covariance=self.config.do_include_astrometric_errors,

1056 )

1057 astrometry_matches, astrometry_meta = self._fit_astrometry(

1058 result.exposure, result.psf_stars_footprints, exposure_region=exposure_region

1059 )

1060 num_astrometry_matches = len(astrometry_matches)

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

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

1063 astrometry_meta)

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

1065

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

1067 # in the adaptive background estimation.

1068 if self.config.doMaskDiffractionSpikes:

1069 self.diffractionSpikeMask.run(result.exposure)

1070

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

1072 if self.config.do_remeasure_star_background and self.config.do_adaptive_threshold_detection:

1073 self._remeasure_star_background(

1074 result,

1075 background_to_photometric_ratio=result.background_to_photometric_ratio,

1076 )

1077

1078 # Run the stars_detection subtask for the photometric calibration.

1079 result.stars_footprints = self._find_stars(

1080 result.exposure,

1081 result.background,

1082 id_generator,

1083 background_to_photometric_ratio=result.background_to_photometric_ratio,

1084 adaptive_det_res_struct=adaptive_det_res_struct,

1085 num_astrometry_matches=num_astrometry_matches,

1086 )

1087 psf = result.exposure.getPsf()

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

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

1090 psfSigma=psfSigma)

1091

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

1093 afwTable.updateSourceCoords(

1094 result.exposure.wcs,

1095 sourceList=result.stars_footprints,

1096 include_covariance=self.config.do_include_astrometric_errors,

1097 )

1098

1099 summary_stat_catalog = result.stars_footprints

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

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

1102

1103 # Validate the astrometric fit. Send in the stars_footprints

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

1105 # a failure.

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

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

1108 have_fit_astrometry = True

1109

1110 result.stars_footprints, photometry_matches, \

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

1112 have_fit_photometry = True

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

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

1115 # table view.

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

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

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

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

1120 summary_stat_catalog = result.stars_footprints

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

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

1123 photometry_meta)

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

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

1126 except pipeBase.AlgorithmError:

1127 if not have_fit_psf:

1128 result.exposure.setPsf(None)

1129 if not have_fit_astrometry:

1130 result.exposure.setWcs(None)

1131 if not have_fit_photometry:

1132 result.exposure.setPhotoCalib(None)

1133 # Summary stat calculations can handle missing components

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

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

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

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

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

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

1140 raise

1141 else:

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

1143

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

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

1146 # by the subtraction.

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

1148 for maskName in self.config.background_stats_ignored_pixel_masks:

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

1150 and maskName not in bgIgnoredPixelMasks):

1151 bgIgnoredPixelMasks += [maskName]

1152

1153 statsCtrl = afwMath.StatisticsControl()

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

1155

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

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

1158 self.metadata["bg_subtracted_skyPixel_instFlux_median"] = median_bg

1159

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

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

1162 self.metadata["bg_subtracted_skyPixel_instFlux_stdev"] = stdev_bg

1163

1164 self.metadata["bg_subtracted_skySource_flux_median"] = (

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

1166 )

1167 self.metadata["bg_subtracted_skySource_flux_stdev"] = (

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

1169 )

1170

1171 if self.config.do_calibrate_pixels:

1172 self._apply_photometry(

1173 result.exposure,

1174 result.background,

1175 background_to_photometric_ratio=result.background_to_photometric_ratio,

1176 )

1177 result.applied_photo_calib = photo_calib

1178 else:

1179 result.applied_photo_calib = None

1180

1181 self._recordMaskedPixelFractions(result.exposure)

1182

1183 if self.config.run_sattle:

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

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

1186 try:

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

1188 historical=self.config.sattle_historical)

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

1190 except requests.exceptions.HTTPError:

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

1192

1193 return result

1194

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

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

1197

1198 Parameters

1199 ----------

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

1201 Exposure to convert to a photometric-flattened image.

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

1203 Flat image that had previously been applied to exposure.

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

1205 Illumination correction image to convert to photometric-flattened

1206 image.

1207

1208 Returns

1209 -------

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

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

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

1213 configured to use the illumination correction.

1214 """

1215 if not self.config.do_illumination_correction:

1216 return None

1217

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

1219 # bfi = image / background_flat

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

1221 # pfi = image / reference_flux_flat

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

1223 # where the illumination correction contains the jacobian

1224 # of the wcs, converting to fluence units.

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

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

1227 # flattened image to the photometric-flattened image:

1228 # bfi / pfi = illumination_correction.

1229

1230 background_to_photometric_ratio = illumination_correction.image.clone()

1231

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

1233 # a photometric-flattened image.

1234 exposure.maskedImage /= background_to_photometric_ratio

1235

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

1237

1238 return background_to_photometric_ratio

1239

1240 def _compute_psf(self, result, id_generator):

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

1242 them, repairing likely cosmic rays before detection.

1243

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

1245 model does not include contributions from cosmic rays.

1246

1247 Parameters

1248 ----------

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

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

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

1252 attributes:

1253

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

1255 Exposure to detect and measure bright stars on.

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

1257 Background that was fit to the exposure during detection.

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

1259 Image to convert photometric-flattened image to

1260 background-flattened image.

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

1262 Object that generates source IDs and provides random seeds.

1263

1264 Returns

1265 -------

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

1267 Catalog of detected bright sources.

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

1269 PSF candidates returned by the psf determiner.

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

1271 Result struct from the adaptive threshold detection.

1272

1273 Notes

1274 -----

1275 This method modifies the exposure, background and

1276 background_to_photometric_ratio attributes of the result struct

1277 in-place.

1278 """

1279 def log_psf(msg, addToMetadata=False):

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

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

1282 metadata.

1283

1284 Parameters

1285 ----------

1286 msg : `str`

1287 Message to prepend the log info with.

1288 addToMetadata : `bool`, optional

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

1290 metadata (the default is False).

1291 """

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

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

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

1295 if result.background is not None:

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

1297 else:

1298 median_background = 0.0

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

1300 msg, sigma, dimensions, median_background)

1301 if addToMetadata:

1302 self.metadata["final_psf_sigma"] = sigma

1303

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

1305 self.config.install_simple_psf.fwhm)

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

1307

1308 result.background = self.psf_subtract_background.run(

1309 exposure=result.exposure,

1310 backgroundToPhotometricRatio=result.background_to_photometric_ratio,

1311 ).background

1312 log_psf("Initial PSF:")

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

1314

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

1316 if not self.config.do_adaptive_threshold_detection:

1317 adaptive_det_res_struct = None

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

1319 # measurement uses the most accurate background-subtraction.

1320 detections = self.psf_detection.run(

1321 table=table,

1322 exposure=result.exposure,

1323 background=result.background,

1324 backgroundToPhotometricRatio=result.background_to_photometric_ratio,

1325 )

1326 else:

1327 initialThreshold = self.config.psf_detection.thresholdValue

1328 initialThresholdMultiplier = self.config.psf_detection.includeThresholdMultiplier

1329 adaptive_det_res_struct = self.psf_adaptive_threshold_detection.run(

1330 table, result.exposure,

1331 initialThreshold=initialThreshold,

1332 initialThresholdMultiplier=initialThresholdMultiplier,

1333 )

1334 detections = adaptive_det_res_struct.detections

1335

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

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

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

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

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

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

1342

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

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

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

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

1347 # revert to it.

1348 if not self.config.do_adaptive_threshold_detection:

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

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

1351 # converge.

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

1353

1354 log_psf("Rerunning with simple PSF:")

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

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

1357 # to use the fitted PSF?

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

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

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

1361 # relevant once DM-39203 is merged?

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

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

1364 # measurement uses the most accurate background-subtraction.

1365 detections = self.psf_detection.run(

1366 table=table,

1367 exposure=result.exposure,

1368 background=result.background,

1369 backgroundToPhotometricRatio=result.background_to_photometric_ratio,

1370 )

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

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

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

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

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

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

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

1378

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

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

1381 # Final measurement with the CRs removed.

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

1383

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

1385 # calibration flux normalization and aperture corrections.

1386 return detections.sources, psf_result.cellSet, adaptive_det_res_struct

1387

1388 def _measure_aperture_correction(self, exposure, bright_sources):

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

1390 previously-measured bright sources.

1391

1392 This function first normalizes the calibration flux and then

1393 the full set of aperture corrections are measured relative

1394 to this normalized calibration flux.

1395

1396 Parameters

1397 ----------

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

1399 Exposure to set the ApCorrMap on.

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

1401 Catalog of detected bright sources; modified to include columns

1402 necessary for point source determination for the aperture

1403 correction calculation.

1404 """

1405 norm_ap_corr_map = self.psf_normalized_calibration_flux.run(

1406 exposure=exposure,

1407 catalog=bright_sources,

1408 ).ap_corr_map

1409

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

1411

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

1413 for key in norm_ap_corr_map:

1414 ap_corr_map[key] = norm_ap_corr_map[key]

1415

1416 exposure.info.setApCorrMap(ap_corr_map)

1417

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

1419 adaptive_det_res_struct=None, num_astrometry_matches=None):

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

1421 PSF, circular aperture, compensated gaussian fluxes.

1422

1423 Parameters

1424 ----------

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

1426 Exposure to detect and measure stars on.

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

1428 Background that was fit to the exposure during detection;

1429 modified in-place during subsequent detection.

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

1431 Object that generates source IDs and provides random seeds.

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

1433 Image to convert photometric-flattened image to

1434 background-flattened image.

1435

1436 Returns

1437 -------

1438 stars : `SourceCatalog`

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

1440 measurements performed on them.

1441 """

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

1443 id_generator.make_table_id_factory())

1444

1445 maxAdaptiveDetIter = 8

1446 threshFactor = 0.2

1447 if adaptive_det_res_struct is not None:

1448 for adaptiveDetIter in range(maxAdaptiveDetIter):

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

1450 if self.config.do_remeasure_star_background:

1451 adaptiveDetectionConfig.reEstimateBackground = False

1452 else:

1453 adaptiveDetectionConfig.reEstimateBackground = True

1454 adaptiveDetectionConfig.includeThresholdMultiplier = 2.0

1455 psfThreshold = (

1456 adaptive_det_res_struct.thresholdValue*adaptive_det_res_struct.includeThresholdMultiplier

1457 )

1458 adaptiveDetectionConfig.thresholdValue = max(

1459 self.config.star_detection.thresholdValue,

1460 threshFactor*psfThreshold/adaptiveDetectionConfig.includeThresholdMultiplier

1461 )

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

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

1464 adaptiveDetIter, adaptiveDetectionConfig.thresholdValue,

1465 adaptiveDetectionConfig.includeThresholdMultiplier)

1466 adaptiveDetectionTask = lsst.meas.algorithms.SourceDetectionTask(

1467 config=adaptiveDetectionConfig

1468 )

1469 detections = adaptiveDetectionTask.run(

1470 table=table,

1471 exposure=exposure,

1472 background=background,

1473 backgroundToPhotometricRatio=background_to_photometric_ratio,

1474 )

1475 nFootprint = len(detections.sources)

1476 nPeak = 0

1477 nIsolated = 0

1478 for src in detections.sources:

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

1480 if nPeakSrc == 1:

1481 nIsolated += 1

1482 nPeak += nPeakSrc

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

1484 if nFootprint > 0:

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

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

1487 adaptiveDetIter, nPeak/nFootprint, nIsolated, minIsolated)

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

1489 threshFactor = max(0.01, 1.5*threshFactor)

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

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

1492 else:

1493 break

1494 else:

1495 threshFactor *= 0.75

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

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

1498 else:

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

1500 # measurement uses the most accurate background-subtraction.

1501 detections = self.star_detection.run(

1502 table=table,

1503 exposure=exposure,

1504 background=background,

1505 backgroundToPhotometricRatio=background_to_photometric_ratio,

1506 )

1507 sources = detections.sources

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

1509

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

1511 if (self.config.do_downsample_footprints

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

1513 if exposure.info.id is None:

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

1515 seed = 0

1516 else:

1517 seed = exposure.info.id & 0xFFFFFFFF

1518

1519 gen = np.random.RandomState(seed)

1520

1521 # Isolate the sky sources before downsampling.

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

1523 indices = gen.choice(

1524 indices,

1525 size=self.config.downsample_max_footprints,

1526 replace=False,

1527 )

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

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

1530

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

1532

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

1534 sel[indices] = True

1535 sources = sources[sel]

1536

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

1538 # non-PSF sources?

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

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

1541 # contiguity for subsequent tasks.

1542 if not sources.isContiguous():

1543 sources = sources.copy(deep=True)

1544

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

1546 self.star_measurement.run(sources, exposure)

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

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

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

1550 )

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

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

1553

1554 # Run the normalization calibration flux task to apply the

1555 # normalization correction to create normalized

1556 # calibration fluxes.

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

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

1559 self.star_catalog_calculation.run(sources)

1560 self.star_set_primary_flags.run(sources)

1561

1562 result = self.star_selector.run(sources)

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

1564 if not result.sourceCat.isContiguous():

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

1566 else:

1567 return result.sourceCat

1568

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

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

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

1572 and the astrometric fit.

1573

1574 Parameters

1575 ----------

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

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

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

1579 related flag fields.

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

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

1582 be updated in-place.

1583

1584 Notes

1585 -----

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

1587 """

1588 control = afwTable.MatchControl()

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

1590 # Return all matched objects, to separate blends.

1591 control.findOnlyClosest = False

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

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

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

1595

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

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

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

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

1600 best = {}

1601 for match_psf, match_stars, d in matches:

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

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

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

1605 matches = list(best.values())

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

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

1608

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

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

1611

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

1613 n_matches, len(psf_stars), len(stars))

1614 self.metadata["matched_psf_star_count"] = n_matches

1615

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

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

1618 # in best above.

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

1620 if n_unique != n_matches:

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

1622

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

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

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

1626 for field in self.psf_fields:

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

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

1629 stars[field] = result

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

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

1632

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

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

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

1636

1637 Parameters

1638 ----------

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

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

1641 Modified to add the fitted skyWcs.

1642 stars : `SourceCatalog`

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

1644 computed from the pixel positions and fitted WCS.

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

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

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

1648 of the exposure.

1649

1650 Returns

1651 -------

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

1653 Reference/stars matches used in the fit.

1654 """

1655 initialWcs = exposure.wcs

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

1657

1658 # Record astrometry stats to metadata

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

1660 if exposure.wcs is not None:

1661 if initialWcs is not None:

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

1663 self.metadata['initial_to_final_wcs'] = (

1664 initialWcs.pixelToSky(center).separation(

1665 exposure.wcs.pixelToSky(center)

1666 ).asArcseconds()

1667 )

1668 else:

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

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

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

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

1673 else:

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

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

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

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

1678

1679 return result.matches, result.matchMeta

1680

1681 def _fit_photometry(self, exposure, stars):

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

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

1684

1685 Parameters

1686 ----------

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

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

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

1690 unchanged.

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

1692 Good stars selected for use in calibration.

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

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

1695 above stars.

1696

1697 Returns

1698 -------

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

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

1701 photoCalib (instFlux columns are retained as well).

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

1703 Reference/stars matches used in the fit.

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

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

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

1707 Photometric calibration that was fit to the star catalog.

1708 """

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

1710 calibrated_stars = result.photoCalib.calibrateCatalog(stars)

1711 exposure.setPhotoCalib(result.photoCalib)

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

1713

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

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

1716 exposure's pixels and an associated background model.

1717

1718 Parameters

1719 ----------

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

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

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

1723 photometric transform will be attached.

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

1725 Background model to convert to nanojansky units in place.

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

1727 Image to convert photometric-flattened image to

1728 background-flattened image.

1729 """

1730 photo_calib = exposure.getPhotoCalib()

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

1732 identity = afwImage.PhotoCalib(1.0,

1733 photo_calib.getCalibrationErr(),

1734 bbox=exposure.getBBox())

1735 exposure.setPhotoCalib(identity)

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

1737

1738 assert photo_calib._isConstant, \

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

1740

1741 for bg in background:

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

1743 # in python.

1744 binned_image = bg[0].getStatsImage()

1745 binned_image *= photo_calib.getCalibrationMean()

1746

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

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

1749 calibrations attached to it.

1750

1751 Parameters

1752 ----------

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

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

1755 Should be in instrumental units with the photometric calibration

1756 attached.

1757 Modified to contain the computed summary statistics.

1758 stars : `SourceCatalog`

1759 Good stars selected used in calibration.

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

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

1762 above stars. Should be in instrumental units.

1763 """

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

1765 exposure.info.setSummaryStats(summary)

1766

1767 def _recordMaskedPixelFractions(self, exposure):

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

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

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

1771

1772 Parameters

1773 ----------

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

1775 The target exposure to calculate masked pixel fractions for.

1776 """

1777

1778 mask = exposure.mask

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

1780 for maskPlane in maskPlanes:

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

1782 evaluateMaskFraction(mask, maskPlane)

1783 )

1784

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

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

1787 and the input camera model.

1788

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

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

1791

1792 Parameters

1793 ----------

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

1795 The exposure whose WCS will be updated.

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

1797 Camera to be used when constructing updated WCS.

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

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

1800 constructing an updated WCS instead of the boresight and

1801 orientation angle from the visit info.

1802 """

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

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

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

1806 else:

1807 self.log.warning(

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

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

1810 exposure.detector.getId())

1811 detector = exposure.detector

1812 if exposure_record is None:

1813 boresight = exposure.visitInfo.getBoresightRaDec()

1814 orientation = exposure.visitInfo.getBoresightRotAngle()

1815 else:

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

1817 orientation = exposure_record.sky_angle * degrees

1818 self.log.info(

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

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

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

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

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

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

1825 )

1826 newWcs = createInitialSkyWcsFromBoresight(boresight, orientation, detector)

1827 exposure.setWcs(newWcs)

1828

1829 def _update_wcs_to_exposure_record(self, exposure, exposure_record):

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

1831 exposure record pointing.

1832

1833 Parameters

1834 ----------

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

1836 The exposure whose WCS will be updated.

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

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

1839 constructing an updated WCS instead of the boresight and

1840 orientation angle from the visit info.

1841 """

1842 detector = exposure.detector

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

1844 orientation = exposure_record.sky_angle * degrees

1845 self.log.info(

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

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

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

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

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

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

1852 )

1853 newWcs = createInitialSkyWcsFromBoresight(boresight, orientation, detector)

1854 exposure.setWcs(newWcs)

1855

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

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

1858

1859 Parameters

1860 ----------

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

1862 The mask on which to evaluate the fraction.

1863 mask_planes : `list`, `str`

1864 The mask planes for which to evaluate the fraction.

1865 bad_mask_planes : `list`, `str`

1866 The mask planes to exclude from the computation.

1867

1868 Returns

1869 -------

1870 detected_fraction : `float`

1871 The calculated fraction of masked pixels

1872 """

1873 bad_pixel_mask = afwImage.Mask.getPlaneBitMask(bad_mask_planes)

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

1875 if n_good_pix == 0:

1876 detected_fraction = float("nan")

1877 return detected_fraction

1878 detected_pixel_mask = afwImage.Mask.getPlaneBitMask(mask_planes)

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

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

1881 detected_fraction = n_detected_pix/n_good_pix

1882 return detected_fraction

1883

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

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

1886

1887 Parameters

1888 ----------

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

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

1891 detected_fraction : `float`

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

1893 full detector.

1894 mask_planes : `list`, `str`

1895 The mask planes for which to evaluate the fraction.

1896 bad_mask_planes : `list`, `str`

1897 The mask planes to exclude from the computation.

1898

1899 Returns

1900 -------

1901 n_above_max_per_amp : `int`

1902 The number of amplifiers with masked fractions above a maximum

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

1904 highest_detected_fraction_per_amp : `float`

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

1906 no_zero_det_amps : `bool`

1907 A boolean representing whether any of the amplifiers has zero

1908 masked pixels.

1909 """

1910 highest_detected_fraction_per_amp = -9.99

1911 n_above_max_per_amp = 0

1912 n_no_zero_det_amps = 0

1913 no_zero_det_amps = True

1914 amps = exposure.detector.getAmplifiers()

1915 if amps is not None:

1916 for ia, amp in enumerate(amps):

1917 amp_bbox = amp.getBBox()

1918 exp_bbox = exposure.getBBox()

1919 if not exp_bbox.contains(amp_bbox):

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

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

1922 continue

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

1924 detected_fraction_amp = self._compute_mask_fraction(sub_image.mask,

1925 mask_planes,

1926 bad_mask_planes)

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

1928 amp.getName(), detected_fraction_amp)

1929 if detected_fraction_amp < 0.002:

1930 n_no_zero_det_amps += 1

1931 if n_no_zero_det_amps > 2:

1932 no_zero_det_amps = False

1933 break

1934 highest_detected_fraction_per_amp = max(detected_fraction_amp,

1935 highest_detected_fraction_per_amp)

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

1937 n_above_max_per_amp += 1

1938 if n_above_max_per_amp > 2:

1939 break

1940 else:

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

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

1943 return n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps

1944

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

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

1947 threshold detection.

1948

1949 Parameters

1950 ----------

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

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

1953 CalibrateImageTask. Will be modified in place.

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

1955 Image to convert photometric-flattened image to

1956 background-flattened image.

1957

1958 Returns

1959 -------

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

1961 The modified result Struct with the new background subtracted.

1962 """

1963 # Restore the previously measured background and remeasure it

1964 # using an adaptive threshold detection iteration to ensure a

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

1966 backgroundOrig = result.background.clone()

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

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

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

1970 result.background = afwMath.BackgroundList()

1971

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

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

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

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

1976 detected_mask_planes,

1977 bad_mask_planes)

1978 minDetFracForFinalBg = 0.02

1979 maxDetFracForFinalBg = 0.93

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

1981 # them in the detection step.

1982 nPixToDilate = 10

1983 maxAdaptiveDetIter = 10

1984 for dilateIter in range(maxAdaptiveDetIter):

1985 dilatedMask = origMask.clone()

1986 for maskName in detected_mask_planes:

1987 # Compute the grown detection mask plane using SpanSet

1988 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)

1989 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)

1990 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)

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

1992 # Clear the detected mask plane

1993 detectedMask = dilatedMask.getMaskPlane(maskName)

1994 dilatedMask.clearMaskPlane(detectedMask)

1995 # Set the mask plane to the dilated one

1996 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)

1997

1998 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,

1999 detected_mask_planes,

2000 bad_mask_planes)

2001 if detected_fraction_dilated < maxDetFracForFinalBg or nPixToDilate == 1:

2002 break

2003 else:

2004 nPixToDilate -= 1

2005

2006 result.exposure.mask = dilatedMask

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

2008 detected_fraction_orig, detected_fraction_dilated)

2009 n_above_max_per_amp = -99

2010 highest_detected_fraction_per_amp = float("nan")

2011

2012 # Check the per-amplifier detection fractions.

2013 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \

2014 self._compute_per_amp_fraction(result.exposure, detected_fraction_dilated,

2015 detected_mask_planes, bad_mask_planes)

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

2017 "highest_detected_fraction_per_amp = %.3f",

2018 n_above_max_per_amp, highest_detected_fraction_per_amp)

2019

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

2021 detected_fraction = 1.0

2022 nFootprintTemp = 1e12

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

2024 for maskName in bgIgnoreMasksToAdd:

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

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

2027 starBackgroundDetectionConfig.background.ignoredPixelMask += [maskName]

2028 starBackgroundDetectionConfig.doTempLocalBackground = False

2029 starBackgroundDetectionConfig.nSigmaToGrow = 70.0

2030 starBackgroundDetectionConfig.reEstimateBackground = False

2031 starBackgroundDetectionConfig.includeThresholdMultiplier = 1.0

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

2033 starBackgroundDetectionConfig.thresholdType = "pixel_stdev"

2034

2035 n_above_max_per_amp = -99

2036 highest_detected_fraction_per_amp = float("nan")

2037 doCheckPerAmpDetFraction = True

2038

2039 minFootprints = self.config.star_background_min_footprints

2040 maxIter = 40

2041 for nIter in range(maxIter):

2042 currentThresh = starBackgroundDetectionConfig.thresholdValue

2043 nZeroEncountered = 0

2044 if nFootprintTemp == 0:

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

2046 starBackgroundDetectionConfig.thresholdValue = zeroFactor*currentThresh

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

2048 starBackgroundDetectionConfig.thresholdValue)

2049 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(

2050 config=starBackgroundDetectionConfig)

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

2052 tempDetections = starBackgroundDetectionTask.run(

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

2054 nFootprintTemp = len(tempDetections.sources)

2055 minFootprints = max(self.config.star_background_min_footprints,

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

2057 minFootprints = min(200, minFootprints)

2058 nZeroEncountered += 1

2059 if nFootprintTemp >= minFootprints:

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

2061 detected_mask_planes,

2062 bad_mask_planes)

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

2064 "DETECTED_NEGATIVE in star_background_detection = %.3f "

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

2066 nIter, starBackgroundDetectionConfig.thresholdValue,

2067 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,

2068 nFootprintTemp, minFootprints)

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

2070

2071 break

2072 else:

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

2074 # entered again.

2075 if nFootprintTemp > 0 and nFootprintTemp < minFootprints:

2076 nFootprintTemp = 0

2077 continue

2078 if detected_fraction > maxDetFracForFinalBg or nFootprintTemp <= minFootprints:

2079 starBackgroundDetectionConfig.thresholdValue = 1.07*currentThresh

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

2081 if nFootprintTemp == 1:

2082 starBackgroundDetectionConfig.thresholdValue = 1.4*currentThresh

2083 else:

2084 starBackgroundDetectionConfig.thresholdValue = 1.2*currentThresh

2085

2086 if n_above_max_per_amp > 1:

2087 starBackgroundDetectionConfig.thresholdValue = 1.1*currentThresh

2088 if detected_fraction < minDetFracForFinalBg:

2089 starBackgroundDetectionConfig.thresholdValue = 0.8*currentThresh

2090 starBackgroundDetectionTask = lsst.meas.algorithms.SourceDetectionTask(

2091 config=starBackgroundDetectionConfig)

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

2093 tempDetections = starBackgroundDetectionTask.run(

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

2095 result.exposure.mask |= dilatedMask

2096 nFootprintTemp = len(tempDetections.sources)

2097 minFootprints = max(self.config.star_background_min_footprints,

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

2099 minFootprints = min(200, minFootprints)

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

2101 bad_mask_planes)

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

2103 "DETECTED_NEGATIVE in star_background_detection = %.3f "

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

2105 nIter, starBackgroundDetectionConfig.thresholdValue,

2106 detected_fraction, maxDetFracForFinalBg, minDetFracForFinalBg,

2107 nFootprintTemp, minFootprints)

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

2109

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

2111 if doCheckPerAmpDetFraction: # detected_fraction < maxDetFracForFinalBg:

2112 n_above_max_per_amp, highest_detected_fraction_per_amp, no_zero_det_amps = \

2113 self._compute_per_amp_fraction(result.exposure, detected_fraction,

2114 detected_mask_planes, bad_mask_planes)

2115

2116 if not no_zero_det_amps:

2117 starBackgroundDetectionConfig.thresholdValue = 0.95*currentThresh

2118

2119 if (detected_fraction < maxDetFracForFinalBg and detected_fraction > minDetFracForFinalBg

2120 and n_above_max_per_amp < int(0.75*n_amp)

2121 and no_zero_det_amps

2122 and nFootprintTemp >= minFootprints):

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

2124 or detected_fraction < 0.85*maxDetFracForFinalBg):

2125 break

2126 else:

2127 starBackgroundDetectionConfig.thresholdValue = 1.05*currentThresh

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

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

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

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

2132 currentThresh, starBackgroundDetectionConfig.thresholdValue)

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

2134

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

2136 bad_mask_planes)

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

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

2139 detected_fraction, highest_detected_fraction_per_amp)

2140

2141 if detected_fraction > maxDetFracForFinalBg:

2142 result.exposure.mask = dilatedMask

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

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

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

2146 detected_fraction, maxDetFracForFinalBg)

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

2148 star_background = self.star_background.run(

2149 exposure=result.exposure,

2150 backgroundToPhotometricRatio=background_to_photometric_ratio,

2151 ).background

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

2153

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

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

2156

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

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

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

2160 # subtraction.

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

2162 nPixToDilate = 2

2163 if detected_fraction < 0.5:

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

2165 for maskName in detected_mask_planes:

2166 # Compute the grown detection mask plane using SpanSet

2167 detectedMaskBit = dilatedMask.getPlaneBitMask(maskName)

2168 detectedMaskSpanSet = SpanSet.fromMask(dilatedMask, detectedMaskBit)

2169 detectedMaskSpanSet = detectedMaskSpanSet.dilated(nPixToDilate)

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

2171 # Clear the detected mask plane

2172 detectedMask = dilatedMask.getMaskPlane(maskName)

2173 dilatedMask.clearMaskPlane(detectedMask)

2174 # Set the mask plane to the dilated one

2175 detectedMaskSpanSet.setMask(dilatedMask, detectedMaskBit)

2176

2177 detected_fraction_dilated = self._compute_mask_fraction(dilatedMask,

2178 detected_mask_planes,

2179 bad_mask_planes)

2180 result.exposure.mask = dilatedMask

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

2182 detected_fraction_orig, detected_fraction_dilated)

2183

2184 bbox = result.exposure.getBBox()

2185

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

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

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

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

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

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

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

2193

2194 # Initialize the parameters for the pedestal iteration.

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

2196 for maskName in bgIgnoreMasksToAdd:

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

2198 and maskName not in pedestalBackgroundConfig.ignoredPixelMask):

2199 pedestalBackgroundConfig.ignoredPixelMask += [maskName]

2200 pedestalBackgroundConfig.statisticsProperty = "MEDIAN"

2201 pedestalBackgroundConfig.approxOrderX = 0

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

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

2204 cumulativePedestalLevel = 0.0

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

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

2207 relativeStoppingCriterion = 0.05

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

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

2210 # converged.

2211 absoluteStoppingCriterion = 0.5

2212 inPedestalIteration = True

2213 while inPedestalIteration:

2214 cumulativePedestalPrev = cumulativePedestalLevel

2215 pedestalBackgroundTask = lsst.meas.algorithms.SubtractBackgroundTask(

2216 config=pedestalBackgroundConfig)

2217 pedestalBackgroundList = pedestalBackgroundTask.run(

2218 exposure=result.exposure,

2219 background=result.background,

2220 backgroundToPhotometricRatio=background_to_photometric_ratio,

2221 ).background

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

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

2224 pedestalBackground = afwMath.BackgroundList()

2225 pedestalBackground.append(pedestalBackgroundList[-1])

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

2227 cumulativePedestalLevel += pedestalBackgroundLevel

2228 absoluteDelta = abs(cumulativePedestalLevel - cumulativePedestalPrev)

2229 if cumulativePedestalPrev != 0.0:

2230 relativeDelta = abs(absoluteDelta/cumulativePedestalPrev)

2231 else:

2232 relativeDelta = 1.0

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

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

2235 pedestalBackgroundLevel, cumulativePedestalLevel,

2236 relativeDelta, absoluteDelta)

2237

2238 if (relativeDelta < relativeStoppingCriterion

2239 or absoluteDelta < absoluteStoppingCriterion):

2240 # We are converged.

2241 inPedestalIteration = False

2242 else:

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

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

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

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

2247 pedestalBackgroundConfig.binSize)

2248 else:

2249 # We have reached the maximum bin size.

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

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

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

2253 relativeStoppingCriterion, absoluteStoppingCriterion)

2254 inPedestalIteration = False

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

2256

2257 # Clear detected mask plane before final round of detection

2258 mask = result.exposure.mask

2259 for mp in detected_mask_planes:

2260 if mp not in mask.getMaskPlaneDict():

2261 mask.addMaskPlane(mp)

2262 mask &= ~mask.getPlaneBitMask(detected_mask_planes)

2263

2264 return result

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

809 statements