Coverage for python / lsst / drp / tasks / update_visit_summary.py: 18%
367 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-18 00:45 -0800
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-18 00:45 -0800
1# This file is part of drp_tasks.
2#
3# Developed for the LSST Data Management System.
4# This product includes software developed by the LSST Project
5# (https://www.lsst.org).
6# See the COPYRIGHT file at the top-level directory of this distribution
7# for details of code ownership.
8#
9# This program is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 3 of the License, or
12# (at your option) any later version.
13#
14# This program is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with this program. If not, see <https://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = (
25 "UpdateVisitSummaryConnections",
26 "UpdateVisitSummaryConfig",
27 "UpdateVisitSummaryTask",
28 "PossiblyMultipleInput",
29 "PerTractInput",
30 "PerHealpixInput",
31 "GlobalInput",
32)
34import dataclasses
35import re
36from abc import ABC, abstractmethod
37from collections.abc import Iterable, Mapping
38from typing import Any
40import astropy.table
41import numpy as np
43import lsst.pipe.base.connectionTypes as cT
44from lsst.afw.cameraGeom import Camera
45from lsst.afw.geom import SkyWcs
46from lsst.afw.image import ExposureSummaryStats
47from lsst.afw.math import BackgroundList
48from lsst.afw.table import ExposureCatalog, ExposureRecord, SchemaMapper
49from lsst.daf.butler import Butler, DatasetRef, DeferredDatasetHandle
50from lsst.geom import Angle, Box2D, Box2I, SpherePoint, arcseconds, degrees
51from lsst.meas.astrom.fit_sip_approximation import FitSipApproximationTask
52from lsst.meas.astrom.refit_pointing import RefitPointingTask
53from lsst.pex.config import ChoiceField, ConfigurableField, Field, ListField
54from lsst.pipe.base import (
55 AlgorithmError,
56 AnnotatedPartialOutputsError,
57 InputQuantizedConnection,
58 InvalidQuantumError,
59 OutputQuantizedConnection,
60 PipelineTask,
61 PipelineTaskConfig,
62 PipelineTaskConnections,
63 QuantumContext,
64 Struct,
65)
66from lsst.pipe.tasks.computeExposureSummaryStats import ComputeExposureSummaryStatsTask
67from lsst.skymap import BaseSkyMap, TractInfo
68from lsst.skymap.detail import makeSkyPolygonFromBBox
69from lsst.sphgeom import HealpixPixelization
72def compute_center_for_detector_record(
73 record: ExposureRecord, bbox: Box2I | None = None, wcs: SkyWcs | None = None
74) -> SpherePoint | None:
75 """Compute the sky coordinate center for a detector to be used when
76 testing distance to tract center.
78 Parameters
79 ----------
80 record : `lsst.afw.table.ExposureRecord`
81 Exposure record to obtain WCS and bbox from if not provided.
82 bbox : `lsst.geom.Box2I`, optional
83 Bounding box for the detector in its own pixel coordinates.
84 wcs : `lsst.afw.geom.SkyWcs`, optional
85 WCS that maps the detector's pixel coordinate system to celestial
86 coordinates.
88 Returns
89 -------
90 center : `lsst.geom.SpherePoint` or `None`
91 Center of the detector in sky coordinates, or `None` if no WCS was
92 given or present in the given record.
93 """
94 if bbox is None:
95 bbox = record.getBBox()
96 if wcs is None:
97 wcs = record.getWcs()
98 if wcs is None:
99 return None
100 region = makeSkyPolygonFromBBox(bbox, wcs)
101 return SpherePoint(region.getCentroid())
104def _validate_wcs_provider(name: str) -> bool:
105 """A helper function to validate the WCS provider at config time.
107 Parameters
108 ----------
109 name : str
110 The source or dimension of the source of the WCS.
112 Returns
113 -------
114 is_valid : bool
115 Whether the provided source is valid.
116 """
117 # Check if ``name`` is a healpix dimension.
118 match_healpix = re.match(re.compile(r"healpix\d+"), name)
119 # Check if ``name`` is one of the other source options.
120 match_other = name in ["input_summary", "global", "tract"]
121 return match_healpix or match_other
124class PossiblyMultipleInput(ABC):
125 """A helper ABC for handling input `~lsst.afw.table.ExposureCatalog`
126 datasets that may be multiple (one per tract/visit combination) or
127 unique/global (one per visit).
128 """
130 @abstractmethod
131 def best_for_detector(
132 self,
133 detector_id: int,
134 center: SpherePoint | None = None,
135 bbox: Box2I | None = None,
136 ) -> tuple[int, ExposureRecord | None]:
137 """Return the exposure record for this detector that is the best match
138 for this detector.
140 Parameters
141 ----------
142 detector_id : `int`
143 Detector ID; used to find the right row in the catalog or catalogs.
144 center : `lsst.geom.SpherePoint` or `None`
145 Center of the detector in sky coordinates. If not provided, one
146 will be computed via `compute_center_for_detector_record`.
147 bbox : `lsst.geom.Box2I`, optional
148 Bounding box for the detector in its own pixel coordinates.
150 Returns
151 -------
152 tract_id : `int`
153 ID of the tract that supplied this record, or `-1` if ``record`` is
154 `None` or if the input was not per-tract.
155 record : `lsst.afw.table.ExposureRecord` or `None`
156 Best record for this detector, or `None` if there either were no
157 records for this detector or no WCS available to compute a center.
158 """
159 raise NotImplementedError()
162@dataclasses.dataclass
163class PerTractInput(PossiblyMultipleInput):
164 """Wrapper class for input `~lsst.afw.table.ExposureCatalog` datasets
165 that are per-tract.
167 This selects the best tract via the minimum average distance (on the sky)
168 from the detector's corners to the tract center.
169 """
171 catalogs_by_tract: list[tuple[TractInfo, ExposureCatalog]]
172 """List of tuples of catalogs and the tracts they correspond to
173 (`list` [`tuple` [`lsst.skymap.TractInfo`,
174 `lsst.afw.table.ExposureCatalog`]]).
175 """
177 @classmethod
178 def load(
179 cls,
180 butler: QuantumContext | Butler,
181 sky_map: BaseSkyMap,
182 refs: Iterable[DatasetRef],
183 ) -> PerTractInput:
184 """Load and wrap input catalogs.
186 Parameters
187 ----------
188 butler : `lsst.pipe.base.QuantumContext`
189 Butler proxy used in `~lsst.pipe.base.PipelineTask.runQuantum`.
190 sky_map : `lsst.skymap.BaseSkyMap`
191 Definition of tracts and patches.
192 refs : `~collections.abc.Iterable` [`lsst.daf.butler.DatasetRef`]
193 References to the catalog datasets to load.
195 Returns
196 -------
197 wrapper : `PerTractInput`
198 Wrapper object for the loaded catalogs.
199 """
200 catalogs_by_tract = []
201 for ref in refs:
202 tract_id = ref.dataId["tract"]
203 tract_info = sky_map[tract_id]
204 catalogs_by_tract.append(
205 (
206 tract_info,
207 butler.get(ref),
208 )
209 )
210 return cls(catalogs_by_tract)
212 def best_for_detector(
213 self,
214 detector_id: int,
215 center: SpherePoint | None = None,
216 bbox: Box2I | None = None,
217 ) -> tuple[int, ExposureRecord | None]:
218 # Docstring inherited.
219 best_result: tuple[int, ExposureRecord | None] = (-1, None)
220 best_distance: Angle = float("inf") * degrees
221 for tract_info, catalog in self.catalogs_by_tract:
222 record = catalog.find(detector_id)
223 if record is None:
224 continue
225 if center is None:
226 center_for_record = compute_center_for_detector_record(record, bbox=bbox)
227 if center_for_record is None:
228 continue
229 else:
230 center_for_record = center
231 center_distance = tract_info.ctr_coord.separation(center_for_record)
232 if best_distance > center_distance:
233 best_result = (tract_info.tract_id, record)
234 best_distance = center_distance
235 return best_result
238@dataclasses.dataclass
239class PerHealpixInput(PossiblyMultipleInput):
240 """Wrapper class for input `~lsst.afw.table.ExposureCatalog` datasets
241 that are per-healpix pixel.
243 This selects the best pixel based on which one contains the detector
244 centroid. Since the pixel regions are ConvexPolygons that do not perfectly
245 describe the pixel borders, in case of ambiguity, the pixel is chosen via
246 the minimum distance (on the sky) from the detector's centroid to the pixel
247 centroid.
248 """
250 catalogs_by_pixel: list[int, tuple[int, SpherePoint, ExposureCatalog]]
251 """List of tuples of catalogs and the pixels they correspond to
252 (`list` [`tuple` [`int`, `lsst.geom.SpherePoint`,
253 `lsst.afw.table.ExposureCatalog`]]).
254 """
256 @classmethod
257 def load(
258 cls,
259 butler: QuantumContext | Butler,
260 refs: Iterable[DatasetRef],
261 ) -> PerHealpixInput:
262 """Load and wrap input catalogs.
264 Parameters
265 ----------
266 butler : `lsst.pipe.base.QuantumContext`
267 Butler proxy used in `~lsst.pipe.base.PipelineTask.runQuantum`.
268 refs : `~collections.abc.Iterable` [`lsst.daf.butler.DatasetRef`]
269 References to the catalog datasets to load.
271 Returns
272 -------
273 wrapper : `PerHealpixInput`
274 Wrapper object for the loaded catalogs.
275 """
276 catalogs_by_pixel = []
277 pixelization = HealpixPixelization(3)
278 for ref in refs:
279 pixel_id = ref.dataId["healpix3"]
280 pixel_region = pixelization.pixel(pixel_id)
281 catalogs_by_pixel.append(
282 (
283 pixel_id,
284 pixel_region,
285 butler.get(ref),
286 )
287 )
288 return cls(catalogs_by_pixel)
290 def best_for_detector(
291 self,
292 detector_id: int,
293 center: SpherePoint | None = None,
294 bbox: Box2I | None = None,
295 ) -> tuple[int, ExposureRecord | None]:
296 # Docstring inherited.
297 best_result: tuple[int, ExposureRecord | None] = (-1, None)
298 inPixel = []
299 distances = []
300 results = []
301 for pixel_id, pixel_region, catalog in self.catalogs_by_pixel:
302 record = catalog.find(detector_id)
303 if record is None:
304 continue
305 if center is None:
306 center_for_record = compute_center_for_detector_record(record, bbox=bbox)
307 if center_for_record is None:
308 continue
309 else:
310 center_for_record = center
312 inPixel.append(pixel_region.contains(center_for_record.getVector()))
313 distances.append(SpherePoint(pixel_region.getCentroid()).separation(center_for_record))
314 results.append((pixel_id, record))
316 if len(results) == 0:
317 # There were no entries at all; return the null result.
318 return best_result
320 inPixel = np.array(inPixel)
321 results = np.array(results)
322 distances = np.array(distances)
323 if inPixel.sum() == 1:
324 # There was one match in the pixel; return it.
325 best_result = results[inPixel][0]
326 elif inPixel.sum() == 0:
327 # There were no entries in the pixel, but choose the closest.
328 best_result = results[distances.argmin()]
329 else:
330 # There were multiple entries in the pixel. Choose the closest.
331 best_result = results[inPixel][distances[inPixel].argmin()]
332 return best_result
335@dataclasses.dataclass
336class GlobalInput(PossiblyMultipleInput):
337 """Wrapper class for input `~lsst.afw.table.ExposureCatalog` datasets
338 that are not per-tract.
339 """
341 catalog: ExposureCatalog
342 """Loaded per-visit catalog dataset (`lsst.afw.table.ExposureCatalog`).
343 """
345 def best_for_detector(
346 self,
347 detector_id: int,
348 center: SpherePoint | None = None,
349 bbox: Box2I | None = None,
350 ) -> tuple[int, ExposureRecord | None]:
351 # Docstring inherited.
352 return -1, self.catalog.find(detector_id)
355class UpdateVisitSummaryConnections(
356 PipelineTaskConnections,
357 dimensions=("instrument", "visit"),
358 defaultTemplates={
359 "skyWcsName": "gbdesAstrometricFit",
360 "photoCalibName": "fgcm",
361 },
362):
363 camera = cT.PrerequisiteInput(
364 doc="Camera geometry.",
365 name="camera",
366 dimensions=("instrument",),
367 storageClass="Camera",
368 isCalibration=True,
369 )
370 sky_map = cT.Input(
371 doc="Description of tract/patch geometry.",
372 name=BaseSkyMap.SKYMAP_DATASET_TYPE_NAME,
373 dimensions=("skymap",),
374 storageClass="SkyMap",
375 )
376 input_summary_schema = cT.InitInput(
377 doc="Schema for input_summary_catalog.",
378 name="visitSummary_schema",
379 storageClass="ExposureCatalog",
380 )
381 input_summary_catalog = cT.Input(
382 doc="Visit summary table to load and modify.",
383 name="visitSummary",
384 dimensions=("instrument", "visit"),
385 storageClass="ExposureCatalog",
386 )
387 input_exposures = cT.Input(
388 doc=(
389 "Per-detector images to obtain image, mask, and variance from "
390 "(embedded summary stats and other components are ignored)."
391 ),
392 name="calexp",
393 dimensions=("instrument", "detector", "visit"),
394 storageClass="ExposureF",
395 multiple=True,
396 deferLoad=True,
397 deferGraphConstraint=True,
398 )
399 psf_overrides = cT.Input(
400 doc="Visit-level catalog of updated PSFs to use.",
401 name="finalized_psf_ap_corr_catalog",
402 dimensions=("instrument", "visit"),
403 storageClass="ExposureCatalog",
404 deferGraphConstraint=True,
405 )
406 psf_star_catalog = cT.Input(
407 doc="Per-visit table of PSF reserved- and used-star measurements.",
408 name="finalized_src_table",
409 dimensions=("instrument", "visit"),
410 storageClass="ArrowAstropy",
411 deferGraphConstraint=True,
412 )
413 ap_corr_overrides = cT.Input(
414 doc="Visit-level catalog of updated aperture correction maps to use.",
415 name="finalized_psf_ap_corr_catalog",
416 dimensions=("instrument", "visit"),
417 storageClass="ExposureCatalog",
418 deferGraphConstraint=True,
419 )
420 photo_calib_overrides_tract = cT.Input(
421 doc="Per-Tract visit-level catalog of updated photometric calibration objects to use.",
422 name="{photoCalibName}PhotoCalibCatalog",
423 dimensions=("instrument", "visit", "tract"),
424 storageClass="ExposureCatalog",
425 multiple=True,
426 deferGraphConstraint=True,
427 )
428 photo_calib_overrides_global = cT.Input(
429 doc="Global visit-level catalog of updated photometric calibration objects to use.",
430 name="{photoCalibName}PhotoCalibCatalog",
431 dimensions=("instrument", "visit"),
432 storageClass="ExposureCatalog",
433 deferGraphConstraint=True,
434 )
435 wcs_overrides_tract = cT.Input(
436 doc="Per-tract visit-level catalog of updated astrometric calibration objects to use.",
437 name="{skyWcsName}SkyWcsCatalog",
438 dimensions=("instrument", "visit", "tract"),
439 storageClass="ExposureCatalog",
440 multiple=True,
441 deferGraphConstraint=True,
442 )
443 wcs_overrides_skypix = cT.Input(
444 doc="Per-skypix pixel visit-level catalog of updated astrometric calibration objects to use.",
445 name="gbdesHealpix3AstrometricFitSkyWcsCatalog",
446 dimensions=("instrument", "visit", "healpix3"),
447 storageClass="ExposureCatalog",
448 multiple=True,
449 deferGraphConstraint=True,
450 )
451 wcs_overrides_global = cT.Input(
452 doc="Global visit-level catalog of updated astrometric calibration objects to use.",
453 name="{skyWcsName}SkyWcsCatalog",
454 dimensions=("instrument", "visit"),
455 storageClass="ExposureCatalog",
456 deferGraphConstraint=True,
457 )
458 background_originals = cT.Input(
459 doc="Per-detector original background that has already been subtracted from 'input_exposures'.",
460 name="calexpBackground",
461 dimensions=("instrument", "visit", "detector"),
462 storageClass="Background",
463 multiple=True,
464 deferLoad=True,
465 deferGraphConstraint=True,
466 )
467 background_overrides = cT.Input(
468 doc="Per-detector background that can be subtracted directly from 'input_exposures'.",
469 name="skyCorr",
470 dimensions=("instrument", "visit", "detector"),
471 storageClass="Background",
472 multiple=True,
473 deferLoad=True,
474 deferGraphConstraint=True,
475 )
476 output_summary_schema = cT.InitOutput(
477 doc="Schema of the output visit summary catalog.",
478 name="finalVisitSummary_schema",
479 storageClass="ExposureCatalog",
480 )
481 output_summary_catalog = cT.Output(
482 doc="Visit-level catalog summarizing all image characterizations and calibrations.",
483 name="finalVisitSummary",
484 dimensions=("instrument", "visit"),
485 storageClass="ExposureCatalog",
486 )
487 visit_geometry = cT.Output(
488 doc="Updated visit[, detector] regions that can be used to update butler dimensions records.",
489 name="visit_geometry",
490 dimensions=("instrument", "visit"),
491 storageClass="VisitGeometry",
492 )
494 def __init__(self, *, config: UpdateVisitSummaryConfig | None = None):
495 super().__init__(config=config)
496 match self.config.wcs_provider:
497 case "input_summary":
498 self.inputs.remove("wcs_overrides_tract")
499 self.inputs.remove("wcs_overrides_global")
500 self.inputs.remove("wcs_overrides_skypix")
501 case "tract":
502 self.inputs.remove("wcs_overrides_global")
503 self.inputs.remove("wcs_overrides_skypix")
504 case "global":
505 self.inputs.remove("wcs_overrides_tract")
506 self.inputs.remove("wcs_overrides_skypix")
507 case x if "healpix" in x:
508 self.inputs.remove("wcs_overrides_tract")
509 self.inputs.remove("wcs_overrides_global")
510 self.wcs_overrides_skypix = dataclasses.replace(
511 self.wcs_overrides_skypix, dimensions=("instrument", "visit", self.config.wcs_provider)
512 )
513 case bad:
514 raise ValueError(f"Invalid value wcs_provider={bad!r}; config was not validated.")
515 match self.config.photo_calib_provider:
516 case "input_summary":
517 self.inputs.remove("photo_calib_overrides_tract")
518 self.inputs.remove("photo_calib_overrides_global")
519 case "tract":
520 self.inputs.remove("photo_calib_overrides_global")
521 case "global":
522 self.inputs.remove("photo_calib_overrides_tract")
523 case bad:
524 raise ValueError(f"Invalid value photo_calib_provider={bad!r}; config was not validated.")
525 match self.config.background_provider:
526 case "input_summary":
527 self.inputs.remove("background_originals")
528 self.inputs.remove("background_overrides")
529 case "replacement":
530 pass
531 case bad:
532 raise ValueError(f"Invalid value background_provider={bad!r}; config was not validated.")
533 if self.config.do_refit_pointing:
534 self.camera = dataclasses.replace(self.camera, dimensions=config.camera_dimensions)
535 else:
536 del self.camera
537 if not self.config.do_write_visit_geometry:
538 del self.visit_geometry
541class UpdateVisitSummaryConfig(PipelineTaskConfig, pipelineConnections=UpdateVisitSummaryConnections):
542 """Configuration for UpdateVisitSummaryTask.
544 Notes
545 -----
546 The configuration defaults for this task reflect a simple or "least common
547 denominator" pipeline, not the more complete, more sophisticated pipeline
548 we run on the instruments we support best. The expectation is that the
549 various full pipeline definitions will generally import the simpler
550 definition, so making the defaults correspond to any full pipeline would
551 just lead to the simple pipeline setting them back to the simple-pipeline
552 values and the full pipeline still having to then override them to the
553 full-pipeline values.
554 """
556 compute_summary_stats = ConfigurableField(
557 doc="Subtask that computes summary statistics from Exposure components.",
558 target=ComputeExposureSummaryStatsTask,
559 )
560 wcs_provider = Field(
561 doc=(
562 "Which connection and behavior to use when applying WCS overrides."
563 "The string should be one of the following:\n"
564 "- input_summary: Propagate the WCS from the input visit summary "
565 "catalog and do not recompute WCS-based summary statistics.\n"
566 "- tract: Use the 'wcs_overrides_tract' connection to load an "
567 "`ExposureCatalog` with {visit, tract} dimensions and per-"
568 "detector rows, and recommpute WCS-based summary statistics.\n"
569 "- global: Use the 'wcs_overrides_global' connection to load an "
570 "`ExposureCatalog` with {visit} dimensions and per-detector rows,"
571 " and recommpute WCS-based summary statistics.\n"
572 "- healpixN: Use the 'wcs_overrides_skypix' connection to load an "
573 "`ExposureCatalog` with {visit, healpixN} dimensions and per-"
574 "detector rows, and recommpute WCS-based summary statistics."
575 ),
576 # If needed, we could add options here to propagate the WCS from
577 # the input exposures and/or transfer WCS-based summary statistics
578 # from them as well. Right now there's no use case for that, since
579 # the input visit summary is always produced after the last time we
580 # write a new Exposure.
581 dtype=str,
582 default="input_summary",
583 optional=False,
584 check=_validate_wcs_provider,
585 )
586 photo_calib_provider = ChoiceField(
587 doc="Which connection and behavior to use when applying photometric calibration overrides.",
588 dtype=str,
589 allowed={
590 "input_summary": (
591 "Propagate the PhotoCalib from the input visit summary catalog "
592 "and do not recompute photometric calibration summary "
593 "statistics."
594 ),
595 "tract": {
596 "Use the 'photo_calib_overrides_tract' connection to load an "
597 "`ExposureCatalog` with {visit, tract} dimensions and per-"
598 "detector rows, and recommpute photometric calibration summary "
599 "statistics."
600 },
601 "global": {
602 "Use the 'photo_calib_overrides_global' connection to load an "
603 "`ExposureCatalog` with {visit} dimensions and per-"
604 "detector rows, and recommpute photometric calibration summary "
605 "statistics."
606 },
607 # If needed, we could add options here to propagate the PhotoCalib
608 # from the input exposures and/or transfer photometric calibration
609 # summary statistics them as well. Right now there's no use case
610 # for that, since the input visit summary is always produced after
611 # the last time we write a new Exposure.
612 },
613 default="input_summary",
614 optional=False,
615 )
616 background_provider = ChoiceField(
617 doc="Which connection(s) and behavior to use when applying background overrides.",
618 dtype=str,
619 allowed={
620 "input_summary": (
621 "The input visit summary catalog already includes summary "
622 "statistics for the final backgrounds that can be used as-is."
623 ),
624 "replacement": {
625 "The 'background_originals' connection refers to a background "
626 "model that has been superseded by the model referred to by "
627 "the 'background_overrides' connection."
628 },
629 # Could also imagine an option in which there is no original
630 # background and the new one stands alone; can add later if needed.
631 },
632 default="input_summary",
633 optional=False,
634 )
635 # Could imagine an option here to say that the original background has not
636 # been subtracted from the input exposures, allowing postISRCCD to be used
637 # as input exposures. Can add later if needed.
638 wcs_consistency_threshold = Field(
639 doc=(
640 "Maximum distance (in arcseconds) between the on-sky corners of a ``wcs_provider`` WCS "
641 "and the ``input_summary`` WCS. When this threshold is exceeded, the WCS is set to None. "
642 "When this is not `None` and there is no input summary WCS, the WCS is also nulled."
643 ),
644 dtype=float,
645 optional=True,
646 default=None,
647 )
648 refit_pointing = ConfigurableField(
649 "A subtask for refitting the boresight pointing and orientation, and fitting WCS SIP approximations.",
650 target=RefitPointingTask,
651 )
652 camera_dimensions = ListField(
653 "The dimensions of the 'camera' prerequisite input connection.",
654 dtype=str,
655 default=["instrument"],
656 )
657 do_refit_pointing = Field(
658 "Whether to re-fit the pointing model and add TAN-SIP WCS approximations.",
659 dtype=bool,
660 default=True,
661 )
662 do_write_visit_geometry = Field(
663 "Whether to write refit-pointing regions that can be used to update butler dimension records.",
664 dtype=bool,
665 default=False,
666 )
667 do_fit_sip_approximations = Field(
668 "Whether to fit TAN-SIP WCS approximations for use in FITS headers. "
669 "This also causes the 'wcs_sip_delta_sky' and 'wcs_sip_delta_pixel' fields to be added "
670 "to track the precision of these approximations.",
671 dtype=bool,
672 default=True,
673 )
674 fit_sip_approximations = ConfigurableField(
675 "Configuration for fitting TAN-SIP WCS approximations for use in FITS headers.",
676 target=FitSipApproximationTask,
677 )
679 def validate(self):
680 super().validate()
681 if self.wcs_consistency_threshold is not None and self.wcs_provider == "input_summary":
682 raise ValueError("Cannot check for WCS consistency if there is only one WCS to consider.")
683 if self.do_write_visit_geometry and not self.do_refit_pointing:
684 raise ValueError("Cannot write visit_geometry without refitting the pointing.")
687class UpdateVisitSummaryTask(PipelineTask):
688 """A pipeline task that creates a new visit-summary table after all
689 `lsst.afw.image.Exposure` components have been finalized.
691 Notes
692 -----
693 This task is designed to be run just prior to making warps for coaddition,
694 as it aggregates all inputs other than the images and backgrounds into a
695 single ``ExposureCatalog`` dataset and recomputes summary statistics that
696 are useful in selecting which images should go into a coadd. Its output
697 can also be used to reconstruct a final processed visit image when combined
698 with a post-ISR image, the background model, and the final mask.
699 """
701 # The `run` method of this task can conditionally apply overrides for PSFs
702 # and aperture corrections, but its `PipelineTask` interface always applies
703 # them. We can always add the config options to make them optional later,
704 # if that turns out to be useful.
706 _DefaultName = "updateVisitSummary"
707 ConfigClass = UpdateVisitSummaryConfig
709 compute_summary_stats: ComputeExposureSummaryStatsTask
711 def __init__(self, *, initInputs: dict[str, Any] | None = None, **kwargs: Any):
712 super().__init__(initInputs=initInputs, **kwargs)
713 self.makeSubtask("compute_summary_stats")
714 if initInputs is None or "input_summary_schema" not in initInputs:
715 raise RuntimeError("Task requires 'input_summary_schema' in initInputs.")
716 input_summary_schema = initInputs["input_summary_schema"].schema
717 self.schema_mapper = SchemaMapper(input_summary_schema)
718 self.schema_mapper.addMinimalSchema(input_summary_schema)
719 self.schema = self.schema_mapper.getOutputSchema()
720 if self.config.wcs_provider == "tract":
721 self.schema.addField("wcsTractId", type="L", doc="ID of the tract that provided the WCS.")
722 elif "healpix" in self.config.wcs_provider:
723 self.schema.addField("wcsSkypixId", type="L", doc="ID of the skypix pixel that provided the WCS.")
724 if self.config.wcs_provider != "input_summary":
725 self.schema.addField(
726 "wcs_corner_max_offset",
727 type="Angle",
728 doc=(
729 "Maximum on-sky difference between the input-summary WCS and the wcs_provider WCS "
730 "at detector corners."
731 ),
732 )
733 if self.config.photo_calib_provider == "tract":
734 self.schema.addField(
735 "photoCalibTractId",
736 type="L",
737 doc="ID of the tract that provided the PhotoCalib.",
738 )
739 if self.config.do_refit_pointing:
740 self.makeSubtask("refit_pointing", schema=self.schema)
741 self.output_summary_schema = ExposureCatalog(self.schema)
742 self._wcs_consistency_threshold = (
743 self.config.wcs_consistency_threshold * arcseconds
744 if self.config.wcs_consistency_threshold is not None
745 else None
746 )
747 if self.config.do_fit_sip_approximations:
748 self.makeSubtask("fit_sip_approximations")
749 self.schema.addField(
750 "wcs_sip_delta_sky",
751 type="Angle",
752 doc="Maximum pixel-to-sky difference between the WCS and its FITS TAN-SIP approximation.",
753 )
754 self.schema.addField(
755 "wcs_sip_delta_pixel",
756 type=np.float64,
757 units="pixel",
758 doc="Maximum sky-to-pixl difference between the WCS and its FITS TAN-SIP approximation.",
759 )
761 def runQuantum(
762 self,
763 butlerQC: QuantumContext,
764 inputRefs: InputQuantizedConnection,
765 outputRefs: OutputQuantizedConnection,
766 ) -> None:
767 # Docstring inherited.
768 sky_map = butlerQC.get(inputRefs.sky_map)
769 del inputRefs.sky_map
770 inputs = {}
771 # Collapse the wcs_override_ and photo_calib_override_ connection pairs
772 # into individual inputs (either ExposureCatalog or PerTractInput
773 # objects).
774 match self.config.wcs_provider:
775 case "tract":
776 inputs["wcs_overrides"] = PerTractInput.load(butlerQC, sky_map, inputRefs.wcs_overrides_tract)
777 del inputRefs.wcs_overrides_tract
778 case "global":
779 inputs["wcs_overrides"] = GlobalInput(butlerQC.get(inputRefs.wcs_overrides_global))
780 del inputRefs.wcs_overrides_global
781 case "input_summary":
782 inputs["wcs_overrides"] = None
783 case x if "healpix" in x:
784 # TODO: match healpix of arbitrary int?
785 inputs["wcs_overrides"] = PerHealpixInput.load(butlerQC, inputRefs.wcs_overrides_skypix)
786 del inputRefs.wcs_overrides_skypix
787 match self.config.photo_calib_provider:
788 case "tract":
789 inputs["photo_calib_overrides"] = PerTractInput.load(
790 butlerQC, sky_map, inputRefs.photo_calib_overrides_tract
791 )
792 del inputRefs.photo_calib_overrides_tract
793 case "global":
794 inputs["photo_calib_overrides"] = GlobalInput(
795 butlerQC.get(inputRefs.photo_calib_overrides_global)
796 )
797 del inputRefs.photo_calib_overrides_global
798 case "input_summary":
799 inputs["photo_calib_overrides"] = None
800 # Load or make DeferredDatasetHandles for everything else.
801 inputs.update(butlerQC.get(inputRefs))
802 deferred_dataset_types = ["input_exposures"]
803 # Handle whether to look for background originals and overrides at all.
804 match self.config.background_provider:
805 case "replacement":
806 deferred_dataset_types.append("background_originals")
807 deferred_dataset_types.append("background_overrides")
808 # Transform the lists of DeferredDatasetHandles for the multiple=True,
809 # deferLoad=True connections into mappings keyed by detector ID.
810 for name in deferred_dataset_types:
811 handles_list = inputs[name]
812 inputs[name] = {handle.dataId["detector"]: handle for handle in handles_list}
813 for record in inputs["input_summary_catalog"]:
814 detector_id = record.getId()
815 if detector_id not in inputs[name]:
816 raise InvalidQuantumError(
817 f"No {name!r} with detector {detector_id} for visit "
818 f"{butlerQC.quantum.dataId['visit']} even though this detector is present "
819 "in the input visit summary catalog. "
820 "This is most likely to occur when the QuantumGraph that includes this task "
821 "was incorrectly generated with an explicit or implicit (from datasets) tract "
822 "constraint."
823 )
824 # Actually run the task and write the results.
825 result = Struct()
826 try:
827 self.run(**inputs, result=result)
828 except AlgorithmError as e:
829 error = AnnotatedPartialOutputsError.annotate(
830 e, self, result.output_summary_catalog, log=self.log
831 )
832 butlerQC.put(result, outputRefs)
833 raise error from e
835 butlerQC.put(result, outputRefs)
837 def run(
838 self,
839 *,
840 result: Struct,
841 input_summary_catalog: ExposureCatalog,
842 input_exposures: Mapping[int, DeferredDatasetHandle],
843 psf_overrides: ExposureCatalog,
844 psf_star_catalog: astropy.table.Table,
845 ap_corr_overrides: ExposureCatalog,
846 camera: Camera | None = None,
847 photo_calib_overrides: PossiblyMultipleInput | None = None,
848 wcs_overrides: PossiblyMultipleInput | None = None,
849 background_originals: Mapping[int, DeferredDatasetHandle] | None = None,
850 background_overrides: Mapping[int, DeferredDatasetHandle] | None = None,
851 ):
852 """Build an updated version of a visit summary catalog.
854 Parameters
855 ----------
856 result : `lsst.pipe.base.Struct`
857 Output struct to modify in-place. On successful return, this will
858 have the following attributes:
860 - ``output_summary_catalog`` (`lsst.afw.table.ExposureCatalog`)
861 The output visit summary catalog.
862 - ``visit_geometry`` (`lsst.obs.base.visit_geometry.VisitGeometry`)
863 Updated visit[, detector] regions that can be used to update
864 butler dimension records.
866 input_summary_catalog : `lsst.afw.table.ExposureCatalog`
867 Input catalog. Each row in this catalog will be used to produce
868 a row in the output catalog. Any override parameter that is `None`
869 will leave the corresponding values unchanged from those in this
870 input catalog.
871 input_exposures : `collections.abc.Mapping` [`int`,
872 `lsst.daf.butler.DeferredDatasetHandle`]
873 Deferred-load objects that fetch `lsst.afw.image.Exposure`
874 instances. Only the image, mask, and variance are used; all other
875 components are assumed to be superceded by at least
876 ``input_summary_catalog`` and probably some ``_overrides``
877 arguments as well. This usually corresponds to the ``calexp``
878 dataset.
879 psf_overrides : `lsst.afw.table.ExposureCatalog`
880 Catalog with attached `lsst.afw.detection.Psf` objects that
881 supersede the input catalog's PSFs.
882 psf_star_catalog : `astropy.table.Table`
883 Table containing PSF stars for use in computing PSF summary
884 statistics.
885 ap_corr_overrides : `lsst.afw.table.ExposureCatalog`
886 Catalog with attached `lsst.afw.image.ApCorrMap` objects that
887 supersede the input catalog's aperture corrections.
888 camera : `lsst.afw.cameraGeom.Camera`, optional
889 Camera geometry.
890 photo_calib_overrides : `PossiblyMultipleInput`, optional
891 Catalog wrappers with attached `lsst.afw.image.PhotoCalib`
892 objects that supersede the input catalog's photometric
893 calibrations.
894 wcs_overrides : `PossiblyMultipleInput`, optional
895 Catalog wrappers with attached `lsst.afw.geom.SkyWcs` objects
896 that supersede the input catalog's astrometric calibrations.
897 background_originals : `collections.abc.Mapping` [`int`,
898 `lsst.daf.butler.DeferredDatasetHandle`], optional
899 Deferred-load objects that fetch `lsst.afw.math.BackgroundList`
900 instances. These should correspond to the background already
901 subtracted from ``input_exposures``. If not provided and
902 ``background_overrides`` is, it is assumed that the background in
903 ``input_exposures`` has not been subtracted. If provided, all keys
904 in ``background_overrides`` must also be present in
905 ``background_originals``.
906 background_overrides : `collections.abc.Mapping` [`int`,
907 `lsst.daf.butler.DeferredDatasetHandle`], optional
908 Deferred-load objects that fetch `lsst.afw.math.BackgroundList`
909 instances. These should correspond to the background that should
910 now be subtracted from``input_exposures`` to yield the final
911 background-subtracted image.
913 Returns
914 -------
915 result : `lsst.pipe.base.Struct`
916 The same output struct passed in.
918 Notes
919 -----
920 If any override parameter is provided but does not have a value for a
921 particular detector, that component will be set to `None` in the
922 returned catalog for that detector and all summary statistics derived
923 from that component will be reset (usually to ``NaN``) as well. Not
924 passing an override parameter at all will instead pass through the
925 original component and values from the input catalog unchanged.
926 """
927 result.output_summary_catalog = ExposureCatalog(self.schema)
928 result.output_summary_catalog.setMetadata(input_summary_catalog.getMetadata())
929 for input_record in input_summary_catalog:
930 detector_id = input_record.getId()
931 output_record = result.output_summary_catalog.addNew()
933 # Make a new ExposureSummaryStats from the input record. These
934 # might be full of NaNs if the summary stats were not computed
935 # originally because initial astrometry or photometry failed, but
936 # if that's possible there should be a flag for it.
937 summary_stats = ExposureSummaryStats.from_record(input_record)
939 # Also copy the input record values to output record; this copies
940 # many of the same values just copied into `summary_stats` (which
941 # will be overridden later by summary_stats.update_record), but it
942 # also copies fields that aren't part of summary_stats, including
943 # the actual components like Psf, Wcs, etc.
944 output_record.assign(input_record, self.schema_mapper)
946 exposure = input_exposures[detector_id].get()
947 bbox = exposure.getBBox()
949 if wcs_overrides:
950 wcs_region, wcs_record = wcs_overrides.best_for_detector(detector_id, bbox=bbox)
951 if wcs_record is not None:
952 wcs = wcs_record.getWcs()
953 else:
954 wcs = None
955 if self.config.wcs_provider == "tract":
956 output_record["wcsTractId"] = wcs_region
957 elif "healpix" in self.config.wcs_provider:
958 output_record["wcsSkypixId"] = wcs_region
959 wcs = self._finish_wcs_assignment(input_record, output_record, wcs)
960 self.compute_summary_stats.update_wcs_stats(
961 summary_stats, wcs, bbox, output_record.getVisitInfo()
962 )
963 else:
964 wcs = input_record.getWcs()
966 sources = None
968 if (psf_record := psf_overrides.find(detector_id)) is not None:
969 psf = psf_record.getPsf()
970 else:
971 psf = None
972 output_record.setPsf(psf)
973 sources = psf_star_catalog[psf_star_catalog["detector"] == detector_id]
974 if len(sources) == 0:
975 sources = None
977 if (ap_corr_record := ap_corr_overrides.find(detector_id)) is not None:
978 ap_corr = ap_corr_record.getApCorrMap()
979 else:
980 ap_corr = None
981 output_record.setApCorrMap(ap_corr)
982 self.compute_summary_stats.update_psf_stats(
983 summary_stats,
984 psf,
985 bbox,
986 sources,
987 image_mask=exposure.mask,
988 image_ap_corr_map=ap_corr,
989 sources_is_astropy=True,
990 )
992 if photo_calib_overrides:
993 center = compute_center_for_detector_record(output_record, bbox, wcs)
994 (
995 photo_calib_tract,
996 photo_calib_record,
997 ) = photo_calib_overrides.best_for_detector(detector_id, center=center)
998 if photo_calib_record is not None:
999 photo_calib = photo_calib_record.getPhotoCalib()
1000 else:
1001 photo_calib = None
1002 if self.config.photo_calib_provider == "tract":
1003 output_record["photoCalibTractId"] = photo_calib_tract
1004 output_record.setPhotoCalib(photo_calib)
1005 self.compute_summary_stats.update_photo_calib_stats(summary_stats, photo_calib)
1007 if background_overrides is not None:
1008 if (handle := background_overrides.get(detector_id)) is not None:
1009 new_bkg = handle.get()
1010 if background_originals is not None:
1011 orig_bkg = background_originals[detector_id].get()
1012 else:
1013 orig_bkg = BackgroundList()
1015 full_bkg = orig_bkg.clone()
1016 for layer in new_bkg:
1017 full_bkg.append(layer)
1018 exposure.image -= new_bkg.getImage()
1019 self.compute_summary_stats.update_background_stats(summary_stats, full_bkg)
1020 self.compute_summary_stats.update_masked_image_stats(
1021 summary_stats, exposure.getMaskedImage()
1022 )
1024 # Update the estimated depth and effective exposure time.
1025 self.compute_summary_stats.update_magnitude_limit_stats(summary_stats, exposure)
1026 self.compute_summary_stats.update_effective_time_stats(summary_stats, exposure)
1028 summary_stats.update_record(output_record)
1029 del exposure
1031 if self.config.do_refit_pointing:
1032 refit_pointing_result = self.refit_pointing.run(
1033 catalog=result.output_summary_catalog,
1034 camera=camera,
1035 )
1036 result.visit_geometry = refit_pointing_result.regions
1038 if self.config.do_fit_sip_approximations:
1039 for output_record in result.output_summary_catalog:
1040 wcs = output_record.getWcs()
1041 if wcs is None:
1042 continue
1043 sip = self.fit_sip_approximations.run(wcs=wcs, bbox=output_record.getBBox())
1044 output_record.setWcs(sip.wcs)
1045 output_record["wcs_sip_delta_sky"] = sip.delta_sky
1046 output_record["wcs_sip_delta_pixel"] = sip.delta_pixel
1048 return result
1050 def _finish_wcs_assignment(
1051 self, input_record: ExposureRecord, output_record: ExposureRecord, updated: SkyWcs
1052 ) -> SkyWcs | None:
1053 preliminary = input_record.getWcs()
1054 if preliminary is None:
1055 # We can't compare to preliminary, so shortcut early.
1056 if self._wcs_consistency_threshold is not None:
1057 # We can't check consistency, so don't trust the updated WCS.
1058 output_record.setWcs(None)
1059 return None
1060 else:
1061 output_record.setWcs(updated)
1062 return updated
1063 if updated is None:
1064 # At present we assume all downstream code is "updated WCS or
1065 # bust", because we know some code (e.g. coaddition) is, and we
1066 # don't have a good way to mark a WCS as fallback preliminary
1067 # downstream.
1068 output_record.setWcs(None)
1069 return None
1070 pixel_corners = Box2D(input_record.getBBox()).getCorners()
1071 preliminary_sky_corners = preliminary.pixelToSky(pixel_corners)
1072 updated_sky_corners = updated.pixelToSky(pixel_corners)
1073 max_offset = max(a.separation(b) for a, b in zip(preliminary_sky_corners, updated_sky_corners))
1074 output_record["wcs_corner_max_offset"] = max_offset
1075 if self._wcs_consistency_threshold is not None and max_offset > self._wcs_consistency_threshold:
1076 output_record.setWcs(None)
1077 return None
1078 output_record.setWcs(updated)
1079 return updated