Coverage for python / lsst / obs / base / defineVisits.py: 25%
419 statements
« prev ^ index » next coverage.py v7.12.0, created at 2025-12-19 01:15 -0800
1# This file is part of obs_base.
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 <http://www.gnu.org/licenses/>.
22from __future__ import annotations
24__all__ = [
25 "DefineVisitsConfig",
26 "DefineVisitsTask",
27 "GroupExposuresConfig",
28 "GroupExposuresTask",
29 "VisitDefinitionData",
30 "VisitSystem",
31]
33import cmath
34import dataclasses
35import enum
36import math
37import operator
38from abc import ABCMeta, abstractmethod
39from collections import defaultdict
40from collections.abc import Callable, Iterable, Sequence
41from typing import Any, ClassVar, TypeVar, cast
43import lsst.geom
44from lsst.afw.cameraGeom import FOCAL_PLANE, PIXELS
45from lsst.daf.butler import Butler, DataCoordinate, DataId, DimensionRecord, Progress, Timespan
46from lsst.geom import Box2D
47from lsst.pex.config import Config, Field, makeRegistry, registerConfigurable
48from lsst.pipe.base import Task
49from lsst.sphgeom import ConvexPolygon, Region, UnitVector3d
50from lsst.utils.introspection import get_full_type_name
52from ._instrument import Instrument, loadCamera
55class VisitSystem(enum.Enum):
56 """Enumeration used to label different visit systems."""
58 ONE_TO_ONE = 0
59 """Each exposure is assigned to its own visit."""
61 BY_GROUP_METADATA = 1
62 """Visit membership is defined by the value of the group dimension or, for
63 older dimension universes, exposure.group_id."""
65 BY_SEQ_START_END = 2
66 """Visit membership is defined by the values of the ``exposure.day_obs``,
67 ``exposure.seq_start``, and ``exposure.seq_end`` values.
68 """
70 @classmethod
71 def all(cls) -> frozenset[VisitSystem]:
72 """Return a `frozenset` containing all members."""
73 return frozenset(cls.__members__.values())
75 @classmethod
76 def from_name(cls, external_name: str) -> VisitSystem:
77 """Construct the enumeration from given name."""
78 name = external_name.upper()
79 name = name.replace("-", "_")
80 try:
81 return cls.__members__[name]
82 except KeyError:
83 raise KeyError(f"Visit system named '{external_name}' not known.") from None
85 @classmethod
86 def from_names(cls, names: Iterable[str] | None) -> frozenset[VisitSystem]:
87 """Return a `frozenset` of all the visit systems matching the supplied
88 names.
90 Parameters
91 ----------
92 names : iterable of `str`, or `None`
93 Names of visit systems. Case insensitive. If `None` or empty, all
94 the visit systems are returned.
96 Returns
97 -------
98 systems : `frozenset` of `VisitSystem`
99 The matching visit systems.
100 """
101 if not names:
102 return cls.all()
104 return frozenset({cls.from_name(name) for name in names})
106 def __str__(self) -> str:
107 name = self.name.lower()
108 name = name.replace("_", "-")
109 return name
112@dataclasses.dataclass
113class VisitDefinitionData:
114 """Struct representing a group of exposures that will be used to define a
115 visit.
116 """
118 instrument: str
119 """Name of the instrument this visit will be associated with.
120 """
122 id: int
123 """Integer ID of the visit.
125 This must be unique across all visit systems for the instrument.
126 """
128 name: str
129 """String name for the visit.
131 This must be unique across all visit systems for the instrument.
132 """
134 visit_systems: set[VisitSystem]
135 """All the visit systems associated with this visit."""
137 exposures: list[DimensionRecord] = dataclasses.field(default_factory=list)
138 """Dimension records for the exposures that are part of this visit.
139 """
142@dataclasses.dataclass
143class _VisitRecords:
144 """Struct containing the dimension records associated with a visit."""
146 visit: DimensionRecord
147 """Record for the 'visit' dimension itself.
148 """
150 visit_definition: list[DimensionRecord]
151 """Records for 'visit_definition', which relates 'visit' to 'exposure'.
152 """
154 visit_detector_region: list[DimensionRecord]
155 """Records for 'visit_detector_region', which associates the combination
156 of a 'visit' and a 'detector' with a region on the sky.
157 """
159 visit_system_membership: list[DimensionRecord]
160 """Records relating visits to an associated visit system."""
163class GroupExposuresConfig(Config):
164 """Configure exposure grouping."""
167class GroupExposuresTask(Task, metaclass=ABCMeta):
168 """Abstract base class for the subtask of `DefineVisitsTask` that is
169 responsible for grouping exposures into visits.
171 Subclasses should be registered with `GroupExposuresTask.registry` to
172 enable use by `DefineVisitsTask`, and should generally correspond to a
173 particular 'visit_system' dimension value. They are also responsible for
174 defining visit IDs and names that are unique across all visit systems in
175 use by an instrument.
177 Parameters
178 ----------
179 config : `GroupExposuresConfig`
180 Configuration information.
181 **kwargs
182 Additional keyword arguments forwarded to the `Task` constructor.
183 """
185 def __init__(self, config: GroupExposuresConfig, **kwargs: Any):
186 Task.__init__(self, config=config, **kwargs)
188 ConfigClass = GroupExposuresConfig
190 _DefaultName = "groupExposures"
192 registry = makeRegistry(
193 doc="Registry of algorithms for grouping exposures into visits.",
194 configBaseType=GroupExposuresConfig,
195 )
197 @abstractmethod
198 def find_missing(
199 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
200 ) -> list[DimensionRecord]:
201 """Determine, if possible, which exposures might be missing.
203 Parameters
204 ----------
205 exposures : `list` of `lsst.daf.butler.DimensionRecord`
206 The exposure records to analyze.
207 registry : `lsst.daf.butler.Registry`
208 A butler registry that contains these exposure records.
210 Returns
211 -------
212 missing : `list` of `lsst.daf.butler.DimensionRecord`
213 Any exposure records present in registry that were related to
214 the given exposures but were missing from that list and deemed
215 to be relevant.
217 Notes
218 -----
219 Only some grouping schemes are able to find missing exposures. It
220 is acceptable to return an empty list.
221 """
222 raise NotImplementedError()
224 @abstractmethod
225 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
226 """Group the exposures in a way most natural for this visit definition.
228 Parameters
229 ----------
230 exposures : `list` of `lsst.daf.butler.DimensionRecord`
231 The exposure records to group.
233 Returns
234 -------
235 groups : `dict` [Any, `list` of `DimensionRecord`]
236 Groupings of exposure records. The key type is relevant to the
237 specific visit definition and could be a string or a tuple.
238 """
239 raise NotImplementedError()
241 @abstractmethod
242 def group(
243 self, exposures: list[DimensionRecord], instrument: Instrument
244 ) -> Iterable[VisitDefinitionData]:
245 """Group the given exposures into visits.
247 Parameters
248 ----------
249 exposures : `list` [ `DimensionRecord` ]
250 DimensionRecords (for the 'exposure' dimension) describing the
251 exposures to group.
252 instrument : `~lsst.obs.base.Instrument`
253 Instrument specification that can be used to optionally support
254 some visit ID definitions.
256 Returns
257 -------
258 visits : `Iterable` [ `VisitDefinitionData` ]
259 Structs identifying the visits and the exposures associated with
260 them. This may be an iterator or a container.
261 """
262 raise NotImplementedError()
264 def getVisitSystems(self) -> set[VisitSystem]:
265 """Return identifiers for the 'visit_system' dimension this
266 algorithm implements.
268 Returns
269 -------
270 visit_systems : `Set` [`VisitSystem`]
271 The visit systems used by this algorithm.
272 """
273 raise NotImplementedError()
276class ComputeVisitRegionsConfig(Config):
277 """Configure visit region calculations."""
279 padding: Field[int] = Field(
280 dtype=int,
281 default=250,
282 doc=(
283 "Pad raw image bounding boxes with specified number of pixels "
284 "when calculating their (conservatively large) region on the "
285 "sky. Note that the config value for pixelMargin of the "
286 "reference object loaders in meas_algorithms should be <= "
287 "the value set here."
288 ),
289 )
292class ComputeVisitRegionsTask(Task, metaclass=ABCMeta):
293 """Abstract base class for the subtask of `DefineVisitsTask` that is
294 responsible for extracting spatial regions for visits and visit+detector
295 combinations.
297 Subclasses should be registered with `ComputeVisitRegionsTask.registry` to
298 enable use by `DefineVisitsTask`.
300 Parameters
301 ----------
302 config : `ComputeVisitRegionsConfig`
303 Configuration information.
304 butler : `lsst.daf.butler.Butler`
305 The butler to use.
306 **kwargs
307 Additional keyword arguments forwarded to the `Task` constructor.
308 """
310 def __init__(self, config: ComputeVisitRegionsConfig, *, butler: Butler, **kwargs: Any):
311 Task.__init__(self, config=config, **kwargs)
312 self.butler = butler
313 self.instrumentMap: dict[str, Instrument] = {}
315 ConfigClass = ComputeVisitRegionsConfig
317 _DefaultName = "computeVisitRegions"
319 registry = makeRegistry(
320 doc="Registry of algorithms for computing on-sky regions for visits and visit+detector combinations.",
321 configBaseType=ComputeVisitRegionsConfig,
322 )
324 def getInstrument(self, instrumentName: str) -> Instrument:
325 """Retrieve an `~lsst.obs.base.Instrument` associated with this
326 instrument name.
328 Parameters
329 ----------
330 instrumentName : `str`
331 The name of the instrument.
333 Returns
334 -------
335 instrument : `~lsst.obs.base.Instrument`
336 The associated instrument object.
338 Notes
339 -----
340 The result is cached.
341 """
342 instrument = self.instrumentMap.get(instrumentName)
343 if instrument is None:
344 instrument = Instrument.fromName(instrumentName, self.butler.registry)
345 self.instrumentMap[instrumentName] = instrument
346 return instrument
348 @abstractmethod
349 def compute(
350 self,
351 visit: VisitDefinitionData,
352 *,
353 collections: Sequence[str] | str | None = None,
354 ) -> tuple[Region, dict[int, Region]]:
355 """Compute regions for the given visit and all detectors in that visit.
357 Parameters
358 ----------
359 visit : `VisitDefinitionData`
360 Struct describing the visit and the exposures associated with it.
361 collections : `Sequence` [ `str` ] or `str` or `None`
362 Collections to be searched for camera geometry, overriding
363 ``self.butler.collections.defaults``. Can be any of the types
364 supported by the ``collections`` argument to butler construction.
366 Returns
367 -------
368 visitRegion : `lsst.sphgeom.Region`
369 Region for the full visit.
370 visitDetectorRegions : `dict` [ `int`, `lsst.sphgeom.Region` ]
371 Dictionary mapping detector ID to the region for that detector.
372 Should include all detectors in the visit.
373 """
374 raise NotImplementedError()
377class DefineVisitsConfig(Config):
378 """Configure visit definition."""
380 groupExposures = GroupExposuresTask.registry.makeField(
381 doc="Algorithm for grouping exposures into visits.",
382 default="one-to-one-and-by-counter",
383 )
384 computeVisitRegions = ComputeVisitRegionsTask.registry.makeField(
385 doc="Algorithm from computing visit and visit+detector regions.",
386 default="single-raw-wcs",
387 )
388 ignoreNonScienceExposures: Field[bool] = Field(
389 doc=(
390 "If True, silently ignore input exposures that do not have "
391 "observation_type=SCIENCE. If False, raise an exception if one "
392 "encountered."
393 ),
394 dtype=bool,
395 optional=False,
396 default=True,
397 )
398 updateObsCoreTable: Field[bool] = Field(
399 doc=(
400 "If True, update exposure regions in obscore table after visits "
401 "are defined. If False, do not update obscore table."
402 ),
403 dtype=bool,
404 default=True,
405 )
408class DefineVisitsTask(Task):
409 """Driver Task for defining visits (and their spatial regions) in Gen3
410 Butler repositories.
412 Parameters
413 ----------
414 config : `DefineVisitsConfig`
415 Configuration for the task.
416 butler : `~lsst.daf.butler.Butler`
417 Writeable butler instance. Will be used to read ``camera`` datasets
418 and insert/sync dimension data.
419 **kwargs
420 Additional keyword arguments are forwarded to the `lsst.pipe.base.Task`
421 constructor.
423 Notes
424 -----
425 Each instance of `DefineVisitsTask` reads from / writes to the same Butler.
426 Each invocation of `DefineVisitsTask.run` processes an independent group of
427 exposures into one or more new visits, all belonging to the same visit
428 system and instrument.
430 The actual work of grouping exposures and computing regions is delegated to
431 pluggable subtasks (`GroupExposuresTask` and `ComputeVisitRegionsTask`),
432 respectively. The defaults are to create one visit for every exposure, and
433 to use exactly one (arbitrary) detector-level raw dataset's WCS along with
434 camera geometry to compute regions for all detectors, but the raw WCS is
435 recomputed from the ``exposure`` dimension record's rotation angle and
436 boresight rather than by loading the ``raw.wcs`` dataset directly. Other
437 implementations can be created and configured for instruments for which
438 these choices are unsuitable (e.g. because visits and exposures are not
439 one-to-one, or because ``raw.wcs`` datasets for different detectors may not
440 be consistent with camera geometry).
442 Defining the same visit the same way multiple times (e.g. via multiple
443 invocations of this task on the same exposures, with the same
444 configuration) is safe, but it may be inefficient, as most of the work must
445 be done before new visits can be compared to existing visits.
446 """
448 def __init__(self, config: DefineVisitsConfig, *, butler: Butler, **kwargs: Any):
449 config.validate() # Not a CmdlineTask nor PipelineTask, so have to validate the config here.
450 super().__init__(config, **kwargs)
451 self.butler = butler
452 self.universe = self.butler.dimensions
453 self.progress = Progress("obs.base.DefineVisitsTask")
454 self.makeSubtask("groupExposures")
455 self.makeSubtask("computeVisitRegions", butler=self.butler)
457 def _reduce_kwargs(self) -> dict:
458 # Add extra parameters to pickle
459 return dict(**super()._reduce_kwargs(), butler=self.butler)
461 ConfigClass: ClassVar[type[Config]] = DefineVisitsConfig
463 _DefaultName: ClassVar[str] = "defineVisits"
465 config: DefineVisitsConfig
466 groupExposures: GroupExposuresTask
467 computeVisitRegions: ComputeVisitRegionsTask
469 def _buildVisitRecords(
470 self, definition: VisitDefinitionData, *, collections: Sequence[str] | str | None = None
471 ) -> _VisitRecords:
472 """Build the DimensionRecords associated with a visit.
474 Parameters
475 ----------
476 definition : `VisitDefinitionData`
477 Struct with identifiers for the visit and records for its
478 constituent exposures.
479 collections : `Sequence` [ `str` ] or `str` or `None`
480 Collections to be searched for camera geometry, overriding
481 ``self.butler.collections.defaults``. Can be any of the types
482 supported by the ``collections`` argument to butler construction.
484 Results
485 -------
486 records : `_VisitRecords`
487 Struct containing DimensionRecords for the visit, including
488 associated dimension elements.
489 """
490 dimension = self.universe["visit"]
492 # Some registries support additional items.
493 supported = {meta.name for meta in dimension.metadata}
495 # Compute all regions.
496 visitRegion, visitDetectorRegions = self.computeVisitRegions.compute(
497 definition, collections=collections
498 )
499 # Aggregate other exposure quantities.
500 timespan = Timespan(
501 begin=_reduceOrNone(min, (e.timespan.begin for e in definition.exposures)),
502 end=_reduceOrNone(max, (e.timespan.end for e in definition.exposures)),
503 )
504 exposure_time = _reduceOrNone(operator.add, (e.exposure_time for e in definition.exposures))
505 physical_filter = _reduceOrNone(_value_if_equal, (e.physical_filter for e in definition.exposures))
506 target_name = _reduceOrNone(_value_if_equal, (e.target_name for e in definition.exposures))
507 science_program = _reduceOrNone(_value_if_equal, (e.science_program for e in definition.exposures))
509 # observing day for a visit is defined by the earliest observation
510 # of the visit
511 observing_day = _reduceOrNone(min, (e.day_obs for e in definition.exposures))
512 observation_reason = _reduceOrNone(
513 _value_if_equal, (e.observation_reason for e in definition.exposures)
514 )
515 if observation_reason is None:
516 # Be explicit about there being multiple reasons
517 observation_reason = "various"
519 # Use the mean zenith angle as an approximation
520 zenith_angle = _reduceOrNone(operator.add, (e.zenith_angle for e in definition.exposures))
521 if zenith_angle is not None:
522 zenith_angle /= len(definition.exposures)
524 # New records that may not be supported.
525 extras: dict[str, Any] = {}
526 if "seq_num" in supported:
527 extras["seq_num"] = _reduceOrNone(min, (e.seq_num for e in definition.exposures))
528 if "azimuth" in supported:
529 # Must take into account 0/360 problem.
530 extras["azimuth"] = _calc_mean_angle([e.azimuth for e in definition.exposures])
532 # visit_system handling changed. This is the logic for visit/exposure
533 # that has support for seq_start/seq_end.
534 if "seq_num" in supported:
535 # Map visit to exposure.
536 visit_definition = [
537 self.universe["visit_definition"].RecordClass(
538 instrument=definition.instrument,
539 visit=definition.id,
540 exposure=exposure.id,
541 )
542 for exposure in definition.exposures
543 ]
545 # Map visit to visit system.
546 visit_system_membership = []
547 for visit_system in self.groupExposures.getVisitSystems():
548 if visit_system in definition.visit_systems:
549 record = self.universe["visit_system_membership"].RecordClass(
550 instrument=definition.instrument,
551 visit=definition.id,
552 visit_system=visit_system.value,
553 )
554 visit_system_membership.append(record)
556 else:
557 # The old approach can only handle one visit system at a time.
558 # If we have been configured with multiple options, prefer the
559 # one-to-one.
560 visit_systems = self.groupExposures.getVisitSystems()
561 if len(visit_systems) > 1:
562 one_to_one = VisitSystem.from_name("one-to-one")
563 if one_to_one not in visit_systems:
564 raise ValueError(
565 f"Multiple visit systems specified ({visit_systems}) for use with old"
566 " dimension universe but unable to find one-to-one."
567 )
568 visit_system = one_to_one
569 else:
570 visit_system = visit_systems.pop()
572 extras["visit_system"] = visit_system.value
574 # The old visit_definition included visit system.
575 visit_definition = [
576 self.universe["visit_definition"].RecordClass(
577 instrument=definition.instrument,
578 visit=definition.id,
579 exposure=exposure.id,
580 visit_system=visit_system.value,
581 )
582 for exposure in definition.exposures
583 ]
585 # This concept does not exist in old schema.
586 visit_system_membership = []
588 # Construct the actual DimensionRecords.
589 return _VisitRecords(
590 visit=dimension.RecordClass(
591 instrument=definition.instrument,
592 id=definition.id,
593 name=definition.name,
594 physical_filter=physical_filter,
595 target_name=target_name,
596 science_program=science_program,
597 observation_reason=observation_reason,
598 day_obs=observing_day,
599 zenith_angle=zenith_angle,
600 exposure_time=exposure_time,
601 timespan=timespan,
602 region=visitRegion,
603 # TODO: no seeing value in exposure dimension records, so we
604 # can't set that here. But there are many other columns that
605 # both dimensions should probably have as well.
606 **extras,
607 ),
608 visit_definition=visit_definition,
609 visit_system_membership=visit_system_membership,
610 visit_detector_region=[
611 self.universe["visit_detector_region"].RecordClass(
612 instrument=definition.instrument,
613 visit=definition.id,
614 detector=detectorId,
615 region=detectorRegion,
616 )
617 for detectorId, detectorRegion in visitDetectorRegions.items()
618 ],
619 )
621 def run(
622 self,
623 dataIds: Iterable[DataId],
624 *,
625 collections: Sequence[str] | str | None = None,
626 update_records: bool = False,
627 incremental: bool = False,
628 ) -> None:
629 """Add visit definitions to the registry for the given exposures.
631 Parameters
632 ----------
633 dataIds : `Iterable` [ `dict` or `~lsst.daf.butler.DataCoordinate` ]
634 Exposure-level data IDs. These must all correspond to the same
635 instrument, and are expected to be on-sky science exposures.
636 collections : `Sequence` [ `str` ] or `str` or `None`
637 Collections to be searched for camera geometry, overriding
638 ``self.butler.collections.defaults``. Can be any of the types
639 supported by the ``collections`` argument to butler construction.
640 update_records : `bool`, optional
641 If `True` (`False` is default), update existing ``visit`` records
642 and ``visit_detector_region`` records. THIS IS AN ADVANCED OPTION
643 THAT SHOULD ONLY BE USED TO FIX REGIONS AND/OR METADATA THAT ARE
644 KNOWN TO BE BAD, AND IT CANNOT BE USED TO REMOVE EXPOSURES OR
645 DETECTORS FROM A VISIT.
646 incremental : `bool`, optional
647 If `True` indicate that exposures are being ingested incrementally
648 and visit definition will be run on partial visits. This will
649 allow the ``visit`` record to be updated if it already exists, but
650 (unlike ``update_records=True``) it will only update the
651 ``visit_detector_region`` records if the ``visit`` record's region
652 changes. If there is any risk that files are being ingested
653 incrementally it is critical that this parameter is set to `True`
654 and not to rely on ``update_records``.
656 Raises
657 ------
658 lsst.daf.butler.registry.ConflictingDefinitionError
659 Raised if a visit ID conflict is detected and the existing visit
660 differs from the new one.
661 """
662 # Normalize, expand, and deduplicate data IDs.
663 self.log.info("Preprocessing data IDs.")
664 dimensions = self.universe.conform(["exposure"])
665 data_id_set: set[DataCoordinate] = {
666 self.butler.registry.expandDataId(d, dimensions=dimensions) for d in dataIds
667 }
668 if not data_id_set:
669 raise RuntimeError("No exposures given.")
670 # Extract exposure DimensionRecords, check that there's only one
671 # instrument in play, and check for non-science exposures.
672 exposures = []
673 instruments = set()
674 instrument_cls_name: str | None = None
675 for dataId in data_id_set:
676 record = dataId.records["exposure"]
677 assert record is not None, "Guaranteed by expandDataIds call earlier."
678 # LSSTCam data can assign ra/dec to flats, and dome-closed
679 # engineering tests. Do not assign a visit if we know that
680 # can_see_sky is False. Treat None as True for this test.
681 can_see_sky = getattr(record, "can_see_sky", True)
682 if (
683 record.tracking_ra is None
684 or record.tracking_dec is None
685 or record.sky_angle is None
686 or can_see_sky is False
687 ):
688 if self.config.ignoreNonScienceExposures:
689 continue
690 else:
691 raise RuntimeError(
692 f"Input exposure {dataId} has observation_type "
693 f"{record.observation_type}, but is not on sky."
694 )
695 instruments.add(dataId["instrument"])
696 instrument_record = dataId.records["instrument"]
697 if instrument_record is not None:
698 instrument_cls_name = instrument_record.class_name
699 exposures.append(record)
700 if not exposures:
701 self.log.info("No on-sky exposures found after filtering.")
702 return
703 if len(instruments) > 1:
704 raise RuntimeError(
705 "All data IDs passed to DefineVisitsTask.run must be "
706 f"from the same instrument; got {instruments}."
707 )
708 (instrument,) = instruments
710 # Might need the instrument class for later depending on universe
711 # and grouping scheme.
712 assert instrument_cls_name is not None, "Instrument must be defined in this dataId"
713 instrument_helper = Instrument.from_string(instrument_cls_name)
715 # Ensure the visit_system our grouping algorithm uses is in the
716 # registry, if it wasn't already.
717 visitSystems = self.groupExposures.getVisitSystems()
718 for visitSystem in visitSystems:
719 self.log.info("Registering visit_system %d: %s.", visitSystem.value, visitSystem)
720 self.butler.registry.syncDimensionData(
721 "visit_system",
722 {"instrument": instrument, "id": visitSystem.value, "name": str(visitSystem)},
723 )
725 # In true incremental we will be given the second snap on its
726 # own on the assumption that the previous snap was already handled.
727 # For correct grouping we need access to the other exposures in the
728 # visit.
729 if incremental:
730 exposures.extend(self.groupExposures.find_missing(exposures, self.butler.registry))
732 # Group exposures into visits, delegating to subtask.
733 self.log.info("Grouping %d exposure(s) into visits.", len(exposures))
734 definitions = list(self.groupExposures.group(exposures, instrument_helper))
735 # Iterate over visits, compute regions, and insert dimension data, one
736 # transaction per visit. If a visit already exists, we skip all other
737 # inserts.
738 self.log.info("Computing regions and other metadata for %d visit(s).", len(definitions))
739 n_skipped: int = 0
740 n_new: int = 0
741 n_fully_updated: int = 0
742 n_partially_updated: int = 0
743 for visitDefinition in self.progress.wrap(
744 definitions, total=len(definitions), desc="Computing regions and inserting visits"
745 ):
746 visitRecords = self._buildVisitRecords(visitDefinition, collections=collections)
747 with self.butler.registry.transaction():
748 inserted_or_updated = self.butler.registry.syncDimensionData(
749 "visit",
750 visitRecords.visit,
751 update=(update_records or incremental),
752 )
753 if inserted_or_updated or update_records:
754 if inserted_or_updated is True:
755 # This is a new visit, not an update to an existing
756 # one, so insert visit definition.
757 # We don't allow visit definitions to change even when
758 # asked to update, because we'd have to delete the old
759 # visit_definitions first and also worry about what
760 # this does to datasets that already use the visit.
761 self.butler.registry.insertDimensionData(
762 "visit_definition", *visitRecords.visit_definition
763 )
764 if visitRecords.visit_system_membership:
765 self.butler.registry.insertDimensionData(
766 "visit_system_membership", *visitRecords.visit_system_membership
767 )
768 elif incremental and len(visitRecords.visit_definition) > 1:
769 # The visit record was modified. This could happen
770 # if a multi-snap visit was redefined with an
771 # additional snap so play it safe and allow for the
772 # visit definition to be updated. We use update=False
773 # here since there should not be any rows updated,
774 # just additional rows added. update=True does not work
775 # correctly with multiple records. In incremental mode
776 # we assume that the caller wants the visit definition
777 # to be updated and has no worries about provenance
778 # with the previous definition.
779 for definition in visitRecords.visit_definition:
780 self.butler.registry.syncDimensionData("visit_definition", definition)
781 if inserted_or_updated is True:
782 # Insert visit-detector regions if the visit is new.
783 self.butler.registry.insertDimensionData(
784 "visit_detector_region",
785 *visitRecords.visit_detector_region,
786 replace=False,
787 )
788 self.log.verbose(
789 "Inserted %s visit_detector_region records for new visit %s.",
790 len(visitRecords.visit_detector_region),
791 visitRecords.visit.id,
792 )
793 n_new += 1
794 # Cast below is because MyPy can't determine that
795 # inserted_or_updated can only be False if update_records
796 # is True.
797 elif update_records or "region" in cast(dict, inserted_or_updated):
798 # Replace visit-detector regions if we were told to
799 # update records explicitly, or if the visit region
800 # changed in an incremental=True update.
801 self.butler.registry.insertDimensionData(
802 "visit_detector_region",
803 *visitRecords.visit_detector_region,
804 replace=True,
805 )
806 self.log.verbose(
807 "Re-inserted %s visit_detector_region records for updated visit %s.",
808 len(visitRecords.visit_detector_region),
809 visitRecords.visit.id,
810 )
811 n_fully_updated += 1
812 else:
813 self.log.verbose(
814 "Updated visit %s without modifying visit_detector_region records.",
815 visitRecords.visit.id,
816 )
817 n_partially_updated += 1
819 # Update obscore exposure records with region information
820 # from corresponding visits.
821 if self.config.updateObsCoreTable:
822 if obscore_manager := self.butler.registry.obsCoreTableManager:
823 obscore_updates: list[tuple[int, int, Region]] = []
824 exposure_ids = [rec.exposure for rec in visitRecords.visit_definition]
825 for record in visitRecords.visit_detector_region:
826 obscore_updates += [
827 (exposure, record.detector, record.region) for exposure in exposure_ids
828 ]
829 if obscore_updates:
830 obscore_manager.update_exposure_regions(
831 cast(str, instrument), obscore_updates
832 )
833 else:
834 self.log.verbose("Skipped already-existing visit %s.", visitRecords.visit.id)
835 n_skipped += 1
836 self.log.info(
837 "Finished writing database records for %d visit(s): %s left unchanged, %s new, "
838 "%s updated with new detector regions, %s updated without new detector regions.",
839 len(definitions),
840 n_skipped,
841 n_new,
842 n_fully_updated,
843 n_partially_updated,
844 )
847_T = TypeVar("_T")
850def _reduceOrNone(func: Callable[[_T, _T], _T | None], iterable: Iterable[_T | None]) -> _T | None:
851 """Apply a binary function to pairs of elements in an iterable until a
852 single value is returned, but return `None` if any element is `None` or
853 there are no elements.
854 """
855 r: _T | None = None
856 for v in iterable:
857 if v is None:
858 return None
859 if r is None:
860 r = v
861 else:
862 r = func(r, v)
863 return r
866def _value_if_equal(a: _T, b: _T) -> _T | None:
867 """Return either argument if they are equal, or `None` if they are not."""
868 return a if a == b else None
871def _calc_mean_angle(angles: list[float]) -> float:
872 """Calculate the mean angle, taking into account 0/360 wrapping.
874 Parameters
875 ----------
876 angles : `list` [`float`]
877 Angles to average together, in degrees.
879 Returns
880 -------
881 average : `float`
882 Average angle in degrees.
883 """
884 # Save on all the math if we only have one value.
885 if len(angles) == 1:
886 return angles[0]
888 # Convert polar coordinates of unit circle to complex values.
889 # Average the complex values.
890 # Convert back to a phase angle.
891 return math.degrees(cmath.phase(sum(cmath.rect(1.0, math.radians(d)) for d in angles) / len(angles)))
894class _GroupExposuresOneToOneConfig(GroupExposuresConfig):
895 visitSystemId: Field[int] = Field(
896 doc="Integer ID of the visit_system implemented by this grouping algorithm.",
897 dtype=int,
898 default=0,
899 deprecated="No longer used. Replaced by enum.",
900 )
901 visitSystemName: Field[str] = Field(
902 doc="String name of the visit_system implemented by this grouping algorithm.",
903 dtype=str,
904 default="one-to-one",
905 deprecated="No longer used. Replaced by enum.",
906 )
909@registerConfigurable("one-to-one", GroupExposuresTask.registry)
910class _GroupExposuresOneToOneTask(GroupExposuresTask, metaclass=ABCMeta):
911 """An exposure grouping algorithm that simply defines one visit for each
912 exposure, reusing the exposures identifiers for the visit.
913 """
915 ConfigClass = _GroupExposuresOneToOneConfig
917 def find_missing(
918 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
919 ) -> list[DimensionRecord]:
920 # By definition no exposures can be missing.
921 return []
923 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
924 # No grouping.
925 return {exposure.id: [exposure] for exposure in exposures}
927 def group(
928 self, exposures: list[DimensionRecord], instrument: Instrument
929 ) -> Iterable[VisitDefinitionData]:
930 # Docstring inherited from GroupExposuresTask.
931 visit_systems = {VisitSystem.from_name("one-to-one")}
932 for exposure in exposures:
933 yield VisitDefinitionData(
934 instrument=exposure.instrument,
935 id=exposure.id,
936 name=exposure.obs_id,
937 exposures=[exposure],
938 visit_systems=visit_systems,
939 )
941 def getVisitSystems(self) -> set[VisitSystem]:
942 # Docstring inherited from GroupExposuresTask.
943 return set(VisitSystem.from_names(["one-to-one"]))
946class _GroupExposuresByGroupMetadataConfig(GroupExposuresConfig):
947 visitSystemId: Field[int] = Field(
948 doc="Integer ID of the visit_system implemented by this grouping algorithm.",
949 dtype=int,
950 default=1,
951 deprecated="No longer used. Replaced by enum.",
952 )
953 visitSystemName: Field[str] = Field(
954 doc="String name of the visit_system implemented by this grouping algorithm.",
955 dtype=str,
956 default="by-group-metadata",
957 deprecated="No longer used. Replaced by enum.",
958 )
961@registerConfigurable("by-group-metadata", GroupExposuresTask.registry)
962class _GroupExposuresByGroupMetadataTask(GroupExposuresTask, metaclass=ABCMeta):
963 """An exposure grouping algorithm that uses the exposure group.
965 This algorithm uses the ``group`` dimension for modern universes and the
966 ``exposure.group_id`` for older universes.
968 This algorithm *assumes* group ID (generally populated from
969 `astro_metadata_translator.ObservationInfo.visit_id`) is not just unique,
970 but disjoint from all `ObservationInfo.exposure_id` values - if it isn't,
971 it will be impossible to ever use both this grouping algorithm and the
972 one-to-one algorithm for a particular camera in the same data repository.
973 """
975 ConfigClass = _GroupExposuresByGroupMetadataConfig
977 def find_missing(
978 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
979 ) -> list[DimensionRecord]:
980 groups = self.group_exposures(exposures)
981 # Determine which group implementation we are using.
982 if "group" in registry.dimensions["exposure"].implied:
983 group_key = "group"
984 else:
985 group_key = "group_name"
986 missing_exposures: list[DimensionRecord] = []
987 for exposures_in_group in groups.values():
988 # We can not tell how many exposures are expected to be in the
989 # visit so we have to query every time.
990 first = exposures_in_group[0]
991 records = set(
992 registry.queryDimensionRecords(
993 "exposure",
994 where=f"exposure.{group_key} = groupnam",
995 bind={"groupnam": getattr(first, group_key)},
996 instrument=first.instrument,
997 )
998 )
999 records.difference_update(set(exposures_in_group))
1000 missing_exposures.extend(list(records))
1001 return missing_exposures
1003 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
1004 groups = defaultdict(list)
1005 group_key = "group"
1006 if exposures and hasattr(exposures[0], "group_name"):
1007 group_key = "group_name"
1008 for exposure in exposures:
1009 groups[getattr(exposure, group_key)].append(exposure)
1010 return groups
1012 def group(
1013 self, exposures: list[DimensionRecord], instrument: Instrument
1014 ) -> Iterable[VisitDefinitionData]:
1015 # Docstring inherited from GroupExposuresTask.
1016 visit_systems = {VisitSystem.from_name("by-group-metadata")}
1017 groups = self.group_exposures(exposures)
1018 has_group_dimension: bool | None = None
1019 for visitName, exposuresInGroup in groups.items():
1020 instrument_name = exposuresInGroup[0].instrument
1021 assert instrument_name == instrument.getName(), "Inconsistency in instrument name"
1022 visit_ids: set[int] = set()
1023 if has_group_dimension is None:
1024 has_group_dimension = hasattr(exposuresInGroup[0], "group")
1025 if has_group_dimension:
1026 visit_ids = {instrument.group_name_to_group_id(e.group) for e in exposuresInGroup}
1027 else:
1028 visit_ids = {e.group_id for e in exposuresInGroup}
1029 assert len(visit_ids) == 1, "Grouping by exposure group does not yield consistent group IDs"
1030 yield VisitDefinitionData(
1031 instrument=instrument_name,
1032 id=visit_ids.pop(),
1033 name=visitName,
1034 exposures=exposuresInGroup,
1035 visit_systems=visit_systems,
1036 )
1038 def getVisitSystems(self) -> set[VisitSystem]:
1039 # Docstring inherited from GroupExposuresTask.
1040 return set(VisitSystem.from_names(["by-group-metadata"]))
1043class _GroupExposuresByCounterAndExposuresConfig(GroupExposuresConfig):
1044 visitSystemId: Field[int] = Field(
1045 doc="Integer ID of the visit_system implemented by this grouping algorithm.",
1046 dtype=int,
1047 default=2,
1048 deprecated="No longer used. Replaced by enum.",
1049 )
1050 visitSystemName: Field[str] = Field(
1051 doc="String name of the visit_system implemented by this grouping algorithm.",
1052 dtype=str,
1053 default="by-counter-and-exposures",
1054 deprecated="No longer used. Replaced by enum.",
1055 )
1058@registerConfigurable("one-to-one-and-by-counter", GroupExposuresTask.registry)
1059class _GroupExposuresByCounterAndExposuresTask(GroupExposuresTask, metaclass=ABCMeta):
1060 """An exposure grouping algorithm that uses the sequence start and
1061 sequence end metadata to create multi-exposure visits, but also
1062 creates one-to-one visits.
1064 This algorithm uses the exposure.seq_start and
1065 exposure.seq_end fields to collect related snaps.
1066 It also groups single exposures.
1067 """
1069 ConfigClass = _GroupExposuresByCounterAndExposuresConfig
1071 def find_missing(
1072 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
1073 ) -> list[DimensionRecord]:
1074 """Analyze the exposures and return relevant exposures known to
1075 registry.
1076 """
1077 groups = self.group_exposures(exposures)
1078 missing_exposures: list[DimensionRecord] = []
1079 for exposures_in_group in groups.values():
1080 sorted_exposures = sorted(exposures_in_group, key=lambda e: e.seq_num)
1081 first = sorted_exposures[0]
1083 # Only need to look for the seq_nums that we don't already have.
1084 seq_nums = set(range(first.seq_start, first.seq_end + 1))
1085 seq_nums.difference_update({exp.seq_num for exp in sorted_exposures})
1087 if seq_nums:
1088 # Missing something. Check registry.
1089 records = list(
1090 registry.queryDimensionRecords(
1091 "exposure",
1092 where="exposure.seq_start = seq_start AND exposure.seq_end = seq_end AND "
1093 "exposure.seq_num IN (seq_nums)",
1094 bind={"seq_start": first.seq_start, "seq_end": first.seq_end, "seq_nums": seq_nums},
1095 instrument=first.instrument,
1096 )
1097 )
1098 missing_exposures.extend(records)
1100 return missing_exposures
1102 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
1103 groups = defaultdict(list)
1104 for exposure in exposures:
1105 groups[exposure.day_obs, exposure.seq_start, exposure.seq_end].append(exposure)
1106 return groups
1108 def group(
1109 self, exposures: list[DimensionRecord], instrument: Instrument
1110 ) -> Iterable[VisitDefinitionData]:
1111 # Docstring inherited from GroupExposuresTask.
1112 system_one_to_one = VisitSystem.from_name("one-to-one")
1113 system_seq_start_end = VisitSystem.from_name("by-seq-start-end")
1115 groups = self.group_exposures(exposures)
1116 for visit_key, exposures_in_group in groups.items():
1117 instrument_name = exposures_in_group[0].instrument
1119 # It is possible that the first exposure in a visit has not
1120 # been ingested. This can be determined and if that is the case
1121 # we can not reliably define the multi-exposure visit.
1122 skip_multi = False
1123 sorted_exposures = sorted(exposures_in_group, key=lambda e: e.seq_num)
1124 first = sorted_exposures.pop(0)
1125 if first.seq_num != first.seq_start:
1126 # Special case seq_num == 0 since that implies that the
1127 # instrument has no counters and therefore no multi-exposure
1128 # visits.
1129 if first.seq_num != 0:
1130 self.log.warning(
1131 "First exposure for visit %s is not present. Skipping the multi-snap definition.",
1132 visit_key,
1133 )
1134 skip_multi = True
1136 multi_exposure = False
1137 if first.seq_start != first.seq_end:
1138 # This is a multi-exposure visit regardless of the number
1139 # of exposures present.
1140 multi_exposure = True
1142 # Define the one-to-one visits.
1143 for exposure in exposures_in_group:
1144 # Default is to use the exposure ID and name unless
1145 # this is the first exposure in a multi-exposure visit.
1146 visit_name = exposure.obs_id
1147 visit_id = exposure.id
1148 visit_systems = {system_one_to_one}
1150 if not multi_exposure:
1151 # This is also a by-counter visit.
1152 # It will use the same visit_name and visit_id.
1153 visit_systems.add(system_seq_start_end)
1155 elif not skip_multi and exposure == first:
1156 # This is the first legitimate exposure in a multi-exposure
1157 # visit. It therefore needs a modified visit name and ID
1158 # so it does not clash with the multi-exposure visit
1159 # definition.
1160 visit_name = f"{visit_name}_first"
1161 visit_id = int(f"9{visit_id}")
1163 yield VisitDefinitionData(
1164 instrument=instrument_name,
1165 id=visit_id,
1166 name=visit_name,
1167 exposures=[exposure],
1168 visit_systems=visit_systems,
1169 )
1171 # Multi-exposure visit.
1172 if not skip_multi and multi_exposure:
1173 # Define the visit using the first exposure
1174 visit_name = first.obs_id
1175 visit_id = first.id
1177 yield VisitDefinitionData(
1178 instrument=instrument_name,
1179 id=visit_id,
1180 name=visit_name,
1181 exposures=exposures_in_group,
1182 visit_systems={system_seq_start_end},
1183 )
1185 def getVisitSystems(self) -> set[VisitSystem]:
1186 # Docstring inherited from GroupExposuresTask.
1187 # Using a Config for this is difficult because what this grouping
1188 # algorithm is doing is using two visit systems.
1189 # One is using metadata (but not by-group) and the other is the
1190 # one-to-one. For now hard-code in class.
1191 return set(VisitSystem.from_names(["one-to-one", "by-seq-start-end"]))
1194class _ComputeVisitRegionsFromSingleRawWcsConfig(ComputeVisitRegionsConfig):
1195 mergeExposures: Field[bool] = Field(
1196 doc=(
1197 "If True, merge per-detector regions over all exposures in a "
1198 "visit (via convex hull) instead of using the first exposure and "
1199 "assuming its regions are valid for all others."
1200 ),
1201 dtype=bool,
1202 default=False,
1203 )
1204 detectorId: Field[int | None] = Field(
1205 doc=(
1206 "Load the WCS for the detector with this ID. If None, use an "
1207 "arbitrary detector (the first found in a query of the data "
1208 "repository for each exposure (or all exposures, if "
1209 "mergeExposures is True)."
1210 ),
1211 dtype=int,
1212 optional=True,
1213 default=None,
1214 )
1215 requireVersionedCamera: Field[bool] = Field(
1216 doc=(
1217 "If True, raise LookupError if version camera geometry cannot be "
1218 "loaded for an exposure. If False, use the nominal camera from "
1219 "the Instrument class instead."
1220 ),
1221 dtype=bool,
1222 optional=False,
1223 default=False,
1224 )
1227@registerConfigurable("single-raw-wcs", ComputeVisitRegionsTask.registry)
1228class _ComputeVisitRegionsFromSingleRawWcsTask(ComputeVisitRegionsTask):
1229 """A visit region calculator that uses a single raw WCS (recomputed from
1230 the ``exposure`` dimension record) and a camera to project the bounding
1231 boxes of all detectors onto the sky, relating different detectors by their
1232 positions in focal plane coordinates.
1234 Notes
1235 -----
1236 Most instruments should have their raw WCSs determined from a combination
1237 of boresight angle, rotator angle, and camera geometry, and hence this
1238 algorithm should produce stable results regardless of which detector the
1239 raw corresponds to. If this is not the case (e.g. because a per-file FITS
1240 WCS is used instead), either the ID of the detector should be fixed (see
1241 the ``detectorId`` config parameter) or a different algorithm used.
1242 """
1244 ConfigClass = _ComputeVisitRegionsFromSingleRawWcsConfig
1245 config: _ComputeVisitRegionsFromSingleRawWcsConfig
1247 def computeExposureBounds(
1248 self, exposure: DimensionRecord, *, collections: Any = None
1249 ) -> dict[int, list[UnitVector3d]]:
1250 """Compute the lists of unit vectors on the sphere that correspond to
1251 the sky positions of detector corners.
1253 Parameters
1254 ----------
1255 exposure : `DimensionRecord`
1256 Dimension record for the exposure.
1257 collections : Any, optional
1258 Collections to be searched for raws and camera geometry, overriding
1259 ``self.butler.collections.defaults``.
1260 Can be any of the types supported by the ``collections`` argument
1261 to butler construction.
1263 Returns
1264 -------
1265 bounds : `dict`
1266 Dictionary mapping detector ID to a list of unit vectors on the
1267 sphere representing that detector's corners projected onto the sky.
1268 """
1269 if collections is None:
1270 collections = list(self.butler.collections.defaults)
1271 camera, versioned = loadCamera(self.butler, exposure.dataId, collections=collections)
1272 if not versioned and self.config.requireVersionedCamera:
1273 raise LookupError(f"No versioned camera found for exposure {exposure.dataId}.")
1275 orientation = lsst.geom.Angle(exposure.sky_angle, lsst.geom.degrees)
1276 radec = lsst.geom.SpherePoint(
1277 lsst.geom.Angle(exposure.tracking_ra, lsst.geom.degrees),
1278 lsst.geom.Angle(exposure.tracking_dec, lsst.geom.degrees),
1279 )
1281 if self.config.detectorId is None:
1282 detectorId = next(camera.getIdIter())
1283 else:
1284 detectorId = self.config.detectorId
1285 wcsDetector = camera[detectorId]
1287 # Ask the raw formatter to create the relevant WCS
1288 # This allows flips to be taken into account
1289 instrument = self.getInstrument(exposure.instrument)
1290 rawFormatter = instrument.getRawFormatter({"detector": detectorId})
1292 try:
1293 wcs = rawFormatter.makeRawSkyWcsFromBoresight(radec, orientation, wcsDetector) # type: ignore
1294 except AttributeError:
1295 raise TypeError(
1296 f"Raw formatter is {get_full_type_name(rawFormatter)} but visit"
1297 " definition requires it to support 'makeRawSkyWcsFromBoresight'"
1298 ) from None
1300 fpToSky = wcsDetector.getTransform(FOCAL_PLANE, PIXELS).then(wcs.getTransform())
1301 bounds = {}
1302 for detector in camera:
1303 pixelsToSky = detector.getTransform(PIXELS, FOCAL_PLANE).then(fpToSky)
1304 pixCorners = Box2D(detector.getBBox().dilatedBy(self.config.padding)).getCorners()
1305 bounds[detector.getId()] = [
1306 skyCorner.getVector() for skyCorner in pixelsToSky.applyForward(pixCorners)
1307 ]
1308 return bounds
1310 def compute(
1311 self, visit: VisitDefinitionData, *, collections: Any = None
1312 ) -> tuple[Region, dict[int, Region]]:
1313 # Docstring inherited from ComputeVisitRegionsTask.
1314 if self.config.mergeExposures:
1315 detectorBounds: dict[int, list[UnitVector3d]] = defaultdict(list)
1316 for exposure in visit.exposures:
1317 exposureDetectorBounds = self.computeExposureBounds(exposure, collections=collections)
1318 for detectorId, bounds in exposureDetectorBounds.items():
1319 detectorBounds[detectorId].extend(bounds)
1320 else:
1321 detectorBounds = self.computeExposureBounds(visit.exposures[0], collections=collections)
1322 visitBounds = []
1323 detectorRegions = {}
1324 for detectorId, bounds in detectorBounds.items():
1325 detectorRegions[detectorId] = ConvexPolygon.convexHull(bounds)
1326 visitBounds.extend(bounds)
1327 return ConvexPolygon.convexHull(visitBounds), detectorRegions