Coverage for python / lsst / obs / base / defineVisits.py: 25%
420 statements
« prev ^ index » next coverage.py v7.13.1, created at 2026-01-10 01:42 -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 Struct, 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 ) -> Struct:
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 Returns
657 -------
658 result : `lsst.pipe.base.Struct`
659 Structure with the following attributes (all `int`):
661 - n_visits: total number of visits defined
662 - n_skipped: number of visits that were already present left alone
663 - n_new: number of new visits inserted
664 - n_fully_updated: number of existing visits fully updated
665 - n_partially_updated: number of visits with non-geometry updates.
667 Raises
668 ------
669 lsst.daf.butler.registry.ConflictingDefinitionError
670 Raised if a visit ID conflict is detected and the existing visit
671 differs from the new one.
672 """
673 # Normalize, expand, and deduplicate data IDs.
674 self.log.info("Preprocessing data IDs.")
675 dimensions = self.universe.conform(["exposure"])
676 data_id_set: set[DataCoordinate] = {
677 self.butler.registry.expandDataId(d, dimensions=dimensions) for d in dataIds
678 }
679 if not data_id_set:
680 raise RuntimeError("No exposures given.")
681 # Extract exposure DimensionRecords, check that there's only one
682 # instrument in play, and check for non-science exposures.
683 exposures = []
684 instruments = set()
685 instrument_cls_name: str | None = None
686 for dataId in data_id_set:
687 record = dataId.records["exposure"]
688 assert record is not None, "Guaranteed by expandDataIds call earlier."
689 # LSSTCam data can assign ra/dec to flats, and dome-closed
690 # engineering tests. Do not assign a visit if we know that
691 # can_see_sky is False. Treat None as True for this test.
692 can_see_sky = getattr(record, "can_see_sky", True)
693 if (
694 record.tracking_ra is None
695 or record.tracking_dec is None
696 or record.sky_angle is None
697 or can_see_sky is False
698 ):
699 if self.config.ignoreNonScienceExposures:
700 continue
701 else:
702 raise RuntimeError(
703 f"Input exposure {dataId} has observation_type "
704 f"{record.observation_type}, but is not on sky."
705 )
706 instruments.add(dataId["instrument"])
707 instrument_record = dataId.records["instrument"]
708 if instrument_record is not None:
709 instrument_cls_name = instrument_record.class_name
710 exposures.append(record)
711 if not exposures:
712 self.log.info("No on-sky exposures found after filtering.")
713 return Struct(n_visits=0, n_skipped=0, n_new=0, n_partially_updated=0, n_fully_updated=0)
714 if len(instruments) > 1:
715 raise RuntimeError(
716 "All data IDs passed to DefineVisitsTask.run must be "
717 f"from the same instrument; got {instruments}."
718 )
719 (instrument,) = instruments
721 # Might need the instrument class for later depending on universe
722 # and grouping scheme.
723 assert instrument_cls_name is not None, "Instrument must be defined in this dataId"
724 instrument_helper = Instrument.from_string(instrument_cls_name)
726 # Ensure the visit_system our grouping algorithm uses is in the
727 # registry, if it wasn't already.
728 visitSystems = self.groupExposures.getVisitSystems()
729 for visitSystem in visitSystems:
730 self.log.info("Registering visit_system %d: %s.", visitSystem.value, visitSystem)
731 self.butler.registry.syncDimensionData(
732 "visit_system",
733 {"instrument": instrument, "id": visitSystem.value, "name": str(visitSystem)},
734 )
736 # In true incremental we will be given the second snap on its
737 # own on the assumption that the previous snap was already handled.
738 # For correct grouping we need access to the other exposures in the
739 # visit.
740 if incremental:
741 exposures.extend(self.groupExposures.find_missing(exposures, self.butler.registry))
743 # Group exposures into visits, delegating to subtask.
744 self.log.info("Grouping %d exposure(s) into visits.", len(exposures))
745 definitions = list(self.groupExposures.group(exposures, instrument_helper))
746 # Iterate over visits, compute regions, and insert dimension data, one
747 # transaction per visit. If a visit already exists, we skip all other
748 # inserts.
749 self.log.info("Computing regions and other metadata for %d visit(s).", len(definitions))
750 n_skipped: int = 0
751 n_new: int = 0
752 n_fully_updated: int = 0
753 n_partially_updated: int = 0
754 for visitDefinition in self.progress.wrap(
755 definitions, total=len(definitions), desc="Computing regions and inserting visits"
756 ):
757 visitRecords = self._buildVisitRecords(visitDefinition, collections=collections)
758 with self.butler.registry.transaction():
759 inserted_or_updated = self.butler.registry.syncDimensionData(
760 "visit",
761 visitRecords.visit,
762 update=(update_records or incremental),
763 )
764 if inserted_or_updated or update_records:
765 if inserted_or_updated is True:
766 # This is a new visit, not an update to an existing
767 # one, so insert visit definition.
768 # We don't allow visit definitions to change even when
769 # asked to update, because we'd have to delete the old
770 # visit_definitions first and also worry about what
771 # this does to datasets that already use the visit.
772 self.butler.registry.insertDimensionData(
773 "visit_definition", *visitRecords.visit_definition
774 )
775 if visitRecords.visit_system_membership:
776 self.butler.registry.insertDimensionData(
777 "visit_system_membership", *visitRecords.visit_system_membership
778 )
779 elif incremental and len(visitRecords.visit_definition) > 1:
780 # The visit record was modified. This could happen
781 # if a multi-snap visit was redefined with an
782 # additional snap so play it safe and allow for the
783 # visit definition to be updated. We use update=False
784 # here since there should not be any rows updated,
785 # just additional rows added. update=True does not work
786 # correctly with multiple records. In incremental mode
787 # we assume that the caller wants the visit definition
788 # to be updated and has no worries about provenance
789 # with the previous definition.
790 for definition in visitRecords.visit_definition:
791 self.butler.registry.syncDimensionData("visit_definition", definition)
792 if inserted_or_updated is True:
793 # Insert visit-detector regions if the visit is new.
794 self.butler.registry.insertDimensionData(
795 "visit_detector_region",
796 *visitRecords.visit_detector_region,
797 replace=False,
798 )
799 self.log.verbose(
800 "Inserted %s visit_detector_region records for new visit %s.",
801 len(visitRecords.visit_detector_region),
802 visitRecords.visit.id,
803 )
804 n_new += 1
805 # Cast below is because MyPy can't determine that
806 # inserted_or_updated can only be False if update_records
807 # is True.
808 elif update_records or "region" in cast(dict, inserted_or_updated):
809 # Replace visit-detector regions if we were told to
810 # update records explicitly, or if the visit region
811 # changed in an incremental=True update.
812 self.butler.registry.insertDimensionData(
813 "visit_detector_region",
814 *visitRecords.visit_detector_region,
815 replace=True,
816 )
817 self.log.verbose(
818 "Re-inserted %s visit_detector_region records for updated visit %s.",
819 len(visitRecords.visit_detector_region),
820 visitRecords.visit.id,
821 )
822 n_fully_updated += 1
823 else:
824 self.log.verbose(
825 "Updated visit %s without modifying visit_detector_region records.",
826 visitRecords.visit.id,
827 )
828 n_partially_updated += 1
830 # Update obscore exposure records with region information
831 # from corresponding visits.
832 if self.config.updateObsCoreTable:
833 if obscore_manager := self.butler.registry.obsCoreTableManager:
834 obscore_updates: list[tuple[int, int, Region]] = []
835 exposure_ids = [rec.exposure for rec in visitRecords.visit_definition]
836 for record in visitRecords.visit_detector_region:
837 obscore_updates += [
838 (exposure, record.detector, record.region) for exposure in exposure_ids
839 ]
840 if obscore_updates:
841 obscore_manager.update_exposure_regions(
842 cast(str, instrument), obscore_updates
843 )
844 else:
845 self.log.verbose("Skipped already-existing visit %s.", visitRecords.visit.id)
846 n_skipped += 1
847 self.log.info(
848 "Finished writing database records for %d visit(s): %s left unchanged, %s new, "
849 "%s updated with new detector regions, %s updated without new detector regions.",
850 len(definitions),
851 n_skipped,
852 n_new,
853 n_fully_updated,
854 n_partially_updated,
855 )
856 return Struct(
857 n_visits=len(definitions),
858 n_skipped=n_skipped,
859 n_new=n_new,
860 n_fully_updated=n_fully_updated,
861 n_partially_updated=n_partially_updated,
862 )
865_T = TypeVar("_T")
868def _reduceOrNone(func: Callable[[_T, _T], _T | None], iterable: Iterable[_T | None]) -> _T | None:
869 """Apply a binary function to pairs of elements in an iterable until a
870 single value is returned, but return `None` if any element is `None` or
871 there are no elements.
872 """
873 r: _T | None = None
874 for v in iterable:
875 if v is None:
876 return None
877 if r is None:
878 r = v
879 else:
880 r = func(r, v)
881 return r
884def _value_if_equal(a: _T, b: _T) -> _T | None:
885 """Return either argument if they are equal, or `None` if they are not."""
886 return a if a == b else None
889def _calc_mean_angle(angles: list[float]) -> float:
890 """Calculate the mean angle, taking into account 0/360 wrapping.
892 Parameters
893 ----------
894 angles : `list` [`float`]
895 Angles to average together, in degrees.
897 Returns
898 -------
899 average : `float`
900 Average angle in degrees.
901 """
902 # Save on all the math if we only have one value.
903 if len(angles) == 1:
904 return angles[0]
906 # Convert polar coordinates of unit circle to complex values.
907 # Average the complex values.
908 # Convert back to a phase angle.
909 return math.degrees(cmath.phase(sum(cmath.rect(1.0, math.radians(d)) for d in angles) / len(angles)))
912class _GroupExposuresOneToOneConfig(GroupExposuresConfig):
913 visitSystemId: Field[int] = Field(
914 doc="Integer ID of the visit_system implemented by this grouping algorithm.",
915 dtype=int,
916 default=0,
917 deprecated="No longer used. Replaced by enum.",
918 )
919 visitSystemName: Field[str] = Field(
920 doc="String name of the visit_system implemented by this grouping algorithm.",
921 dtype=str,
922 default="one-to-one",
923 deprecated="No longer used. Replaced by enum.",
924 )
927@registerConfigurable("one-to-one", GroupExposuresTask.registry)
928class _GroupExposuresOneToOneTask(GroupExposuresTask, metaclass=ABCMeta):
929 """An exposure grouping algorithm that simply defines one visit for each
930 exposure, reusing the exposures identifiers for the visit.
931 """
933 ConfigClass = _GroupExposuresOneToOneConfig
935 def find_missing(
936 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
937 ) -> list[DimensionRecord]:
938 # By definition no exposures can be missing.
939 return []
941 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
942 # No grouping.
943 return {exposure.id: [exposure] for exposure in exposures}
945 def group(
946 self, exposures: list[DimensionRecord], instrument: Instrument
947 ) -> Iterable[VisitDefinitionData]:
948 # Docstring inherited from GroupExposuresTask.
949 visit_systems = {VisitSystem.from_name("one-to-one")}
950 for exposure in exposures:
951 yield VisitDefinitionData(
952 instrument=exposure.instrument,
953 id=exposure.id,
954 name=exposure.obs_id,
955 exposures=[exposure],
956 visit_systems=visit_systems,
957 )
959 def getVisitSystems(self) -> set[VisitSystem]:
960 # Docstring inherited from GroupExposuresTask.
961 return set(VisitSystem.from_names(["one-to-one"]))
964class _GroupExposuresByGroupMetadataConfig(GroupExposuresConfig):
965 visitSystemId: Field[int] = Field(
966 doc="Integer ID of the visit_system implemented by this grouping algorithm.",
967 dtype=int,
968 default=1,
969 deprecated="No longer used. Replaced by enum.",
970 )
971 visitSystemName: Field[str] = Field(
972 doc="String name of the visit_system implemented by this grouping algorithm.",
973 dtype=str,
974 default="by-group-metadata",
975 deprecated="No longer used. Replaced by enum.",
976 )
979@registerConfigurable("by-group-metadata", GroupExposuresTask.registry)
980class _GroupExposuresByGroupMetadataTask(GroupExposuresTask, metaclass=ABCMeta):
981 """An exposure grouping algorithm that uses the exposure group.
983 This algorithm uses the ``group`` dimension for modern universes and the
984 ``exposure.group_id`` for older universes.
986 This algorithm *assumes* group ID (generally populated from
987 `astro_metadata_translator.ObservationInfo.visit_id`) is not just unique,
988 but disjoint from all `ObservationInfo.exposure_id` values - if it isn't,
989 it will be impossible to ever use both this grouping algorithm and the
990 one-to-one algorithm for a particular camera in the same data repository.
991 """
993 ConfigClass = _GroupExposuresByGroupMetadataConfig
995 def find_missing(
996 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
997 ) -> list[DimensionRecord]:
998 groups = self.group_exposures(exposures)
999 # Determine which group implementation we are using.
1000 if "group" in registry.dimensions["exposure"].implied:
1001 group_key = "group"
1002 else:
1003 group_key = "group_name"
1004 missing_exposures: list[DimensionRecord] = []
1005 for exposures_in_group in groups.values():
1006 # We can not tell how many exposures are expected to be in the
1007 # visit so we have to query every time.
1008 first = exposures_in_group[0]
1009 records = set(
1010 registry.queryDimensionRecords(
1011 "exposure",
1012 where=f"exposure.{group_key} = groupnam",
1013 bind={"groupnam": getattr(first, group_key)},
1014 instrument=first.instrument,
1015 )
1016 )
1017 records.difference_update(set(exposures_in_group))
1018 missing_exposures.extend(list(records))
1019 return missing_exposures
1021 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
1022 groups = defaultdict(list)
1023 group_key = "group"
1024 if exposures and hasattr(exposures[0], "group_name"):
1025 group_key = "group_name"
1026 for exposure in exposures:
1027 groups[getattr(exposure, group_key)].append(exposure)
1028 return groups
1030 def group(
1031 self, exposures: list[DimensionRecord], instrument: Instrument
1032 ) -> Iterable[VisitDefinitionData]:
1033 # Docstring inherited from GroupExposuresTask.
1034 visit_systems = {VisitSystem.from_name("by-group-metadata")}
1035 groups = self.group_exposures(exposures)
1036 has_group_dimension: bool | None = None
1037 for visitName, exposuresInGroup in groups.items():
1038 instrument_name = exposuresInGroup[0].instrument
1039 assert instrument_name == instrument.getName(), "Inconsistency in instrument name"
1040 visit_ids: set[int] = set()
1041 if has_group_dimension is None:
1042 has_group_dimension = hasattr(exposuresInGroup[0], "group")
1043 if has_group_dimension:
1044 visit_ids = {instrument.group_name_to_group_id(e.group) for e in exposuresInGroup}
1045 else:
1046 visit_ids = {e.group_id for e in exposuresInGroup}
1047 assert len(visit_ids) == 1, "Grouping by exposure group does not yield consistent group IDs"
1048 yield VisitDefinitionData(
1049 instrument=instrument_name,
1050 id=visit_ids.pop(),
1051 name=visitName,
1052 exposures=exposuresInGroup,
1053 visit_systems=visit_systems,
1054 )
1056 def getVisitSystems(self) -> set[VisitSystem]:
1057 # Docstring inherited from GroupExposuresTask.
1058 return set(VisitSystem.from_names(["by-group-metadata"]))
1061class _GroupExposuresByCounterAndExposuresConfig(GroupExposuresConfig):
1062 visitSystemId: Field[int] = Field(
1063 doc="Integer ID of the visit_system implemented by this grouping algorithm.",
1064 dtype=int,
1065 default=2,
1066 deprecated="No longer used. Replaced by enum.",
1067 )
1068 visitSystemName: Field[str] = Field(
1069 doc="String name of the visit_system implemented by this grouping algorithm.",
1070 dtype=str,
1071 default="by-counter-and-exposures",
1072 deprecated="No longer used. Replaced by enum.",
1073 )
1076@registerConfigurable("one-to-one-and-by-counter", GroupExposuresTask.registry)
1077class _GroupExposuresByCounterAndExposuresTask(GroupExposuresTask, metaclass=ABCMeta):
1078 """An exposure grouping algorithm that uses the sequence start and
1079 sequence end metadata to create multi-exposure visits, but also
1080 creates one-to-one visits.
1082 This algorithm uses the exposure.seq_start and
1083 exposure.seq_end fields to collect related snaps.
1084 It also groups single exposures.
1085 """
1087 ConfigClass = _GroupExposuresByCounterAndExposuresConfig
1089 def find_missing(
1090 self, exposures: list[DimensionRecord], registry: lsst.daf.butler.Registry
1091 ) -> list[DimensionRecord]:
1092 """Analyze the exposures and return relevant exposures known to
1093 registry.
1094 """
1095 groups = self.group_exposures(exposures)
1096 missing_exposures: list[DimensionRecord] = []
1097 for exposures_in_group in groups.values():
1098 sorted_exposures = sorted(exposures_in_group, key=lambda e: e.seq_num)
1099 first = sorted_exposures[0]
1101 # Only need to look for the seq_nums that we don't already have.
1102 seq_nums = set(range(first.seq_start, first.seq_end + 1))
1103 seq_nums.difference_update({exp.seq_num for exp in sorted_exposures})
1105 if seq_nums:
1106 # Missing something. Check registry.
1107 records = list(
1108 registry.queryDimensionRecords(
1109 "exposure",
1110 where="exposure.seq_start = seq_start AND exposure.seq_end = seq_end AND "
1111 "exposure.seq_num IN (seq_nums)",
1112 bind={"seq_start": first.seq_start, "seq_end": first.seq_end, "seq_nums": seq_nums},
1113 instrument=first.instrument,
1114 )
1115 )
1116 missing_exposures.extend(records)
1118 return missing_exposures
1120 def group_exposures(self, exposures: list[DimensionRecord]) -> dict[Any, list[DimensionRecord]]:
1121 groups = defaultdict(list)
1122 for exposure in exposures:
1123 groups[exposure.day_obs, exposure.seq_start, exposure.seq_end].append(exposure)
1124 return groups
1126 def group(
1127 self, exposures: list[DimensionRecord], instrument: Instrument
1128 ) -> Iterable[VisitDefinitionData]:
1129 # Docstring inherited from GroupExposuresTask.
1130 system_one_to_one = VisitSystem.from_name("one-to-one")
1131 system_seq_start_end = VisitSystem.from_name("by-seq-start-end")
1133 groups = self.group_exposures(exposures)
1134 for visit_key, exposures_in_group in groups.items():
1135 instrument_name = exposures_in_group[0].instrument
1137 # It is possible that the first exposure in a visit has not
1138 # been ingested. This can be determined and if that is the case
1139 # we can not reliably define the multi-exposure visit.
1140 skip_multi = False
1141 sorted_exposures = sorted(exposures_in_group, key=lambda e: e.seq_num)
1142 first = sorted_exposures.pop(0)
1143 if first.seq_num != first.seq_start:
1144 # Special case seq_num == 0 since that implies that the
1145 # instrument has no counters and therefore no multi-exposure
1146 # visits.
1147 if first.seq_num != 0:
1148 self.log.warning(
1149 "First exposure for visit %s is not present. Skipping the multi-snap definition.",
1150 visit_key,
1151 )
1152 skip_multi = True
1154 multi_exposure = False
1155 if first.seq_start != first.seq_end:
1156 # This is a multi-exposure visit regardless of the number
1157 # of exposures present.
1158 multi_exposure = True
1160 # Define the one-to-one visits.
1161 for exposure in exposures_in_group:
1162 # Default is to use the exposure ID and name unless
1163 # this is the first exposure in a multi-exposure visit.
1164 visit_name = exposure.obs_id
1165 visit_id = exposure.id
1166 visit_systems = {system_one_to_one}
1168 if not multi_exposure:
1169 # This is also a by-counter visit.
1170 # It will use the same visit_name and visit_id.
1171 visit_systems.add(system_seq_start_end)
1173 elif not skip_multi and exposure == first:
1174 # This is the first legitimate exposure in a multi-exposure
1175 # visit. It therefore needs a modified visit name and ID
1176 # so it does not clash with the multi-exposure visit
1177 # definition.
1178 visit_name = f"{visit_name}_first"
1179 visit_id = int(f"9{visit_id}")
1181 yield VisitDefinitionData(
1182 instrument=instrument_name,
1183 id=visit_id,
1184 name=visit_name,
1185 exposures=[exposure],
1186 visit_systems=visit_systems,
1187 )
1189 # Multi-exposure visit.
1190 if not skip_multi and multi_exposure:
1191 # Define the visit using the first exposure
1192 visit_name = first.obs_id
1193 visit_id = first.id
1195 yield VisitDefinitionData(
1196 instrument=instrument_name,
1197 id=visit_id,
1198 name=visit_name,
1199 exposures=exposures_in_group,
1200 visit_systems={system_seq_start_end},
1201 )
1203 def getVisitSystems(self) -> set[VisitSystem]:
1204 # Docstring inherited from GroupExposuresTask.
1205 # Using a Config for this is difficult because what this grouping
1206 # algorithm is doing is using two visit systems.
1207 # One is using metadata (but not by-group) and the other is the
1208 # one-to-one. For now hard-code in class.
1209 return set(VisitSystem.from_names(["one-to-one", "by-seq-start-end"]))
1212class _ComputeVisitRegionsFromSingleRawWcsConfig(ComputeVisitRegionsConfig):
1213 mergeExposures: Field[bool] = Field(
1214 doc=(
1215 "If True, merge per-detector regions over all exposures in a "
1216 "visit (via convex hull) instead of using the first exposure and "
1217 "assuming its regions are valid for all others."
1218 ),
1219 dtype=bool,
1220 default=False,
1221 )
1222 detectorId: Field[int | None] = Field(
1223 doc=(
1224 "Load the WCS for the detector with this ID. If None, use an "
1225 "arbitrary detector (the first found in a query of the data "
1226 "repository for each exposure (or all exposures, if "
1227 "mergeExposures is True)."
1228 ),
1229 dtype=int,
1230 optional=True,
1231 default=None,
1232 )
1233 requireVersionedCamera: Field[bool] = Field(
1234 doc=(
1235 "If True, raise LookupError if version camera geometry cannot be "
1236 "loaded for an exposure. If False, use the nominal camera from "
1237 "the Instrument class instead."
1238 ),
1239 dtype=bool,
1240 optional=False,
1241 default=False,
1242 )
1245@registerConfigurable("single-raw-wcs", ComputeVisitRegionsTask.registry)
1246class _ComputeVisitRegionsFromSingleRawWcsTask(ComputeVisitRegionsTask):
1247 """A visit region calculator that uses a single raw WCS (recomputed from
1248 the ``exposure`` dimension record) and a camera to project the bounding
1249 boxes of all detectors onto the sky, relating different detectors by their
1250 positions in focal plane coordinates.
1252 Notes
1253 -----
1254 Most instruments should have their raw WCSs determined from a combination
1255 of boresight angle, rotator angle, and camera geometry, and hence this
1256 algorithm should produce stable results regardless of which detector the
1257 raw corresponds to. If this is not the case (e.g. because a per-file FITS
1258 WCS is used instead), either the ID of the detector should be fixed (see
1259 the ``detectorId`` config parameter) or a different algorithm used.
1260 """
1262 ConfigClass = _ComputeVisitRegionsFromSingleRawWcsConfig
1263 config: _ComputeVisitRegionsFromSingleRawWcsConfig
1265 def computeExposureBounds(
1266 self, exposure: DimensionRecord, *, collections: Any = None
1267 ) -> dict[int, list[UnitVector3d]]:
1268 """Compute the lists of unit vectors on the sphere that correspond to
1269 the sky positions of detector corners.
1271 Parameters
1272 ----------
1273 exposure : `DimensionRecord`
1274 Dimension record for the exposure.
1275 collections : Any, optional
1276 Collections to be searched for raws and camera geometry, overriding
1277 ``self.butler.collections.defaults``.
1278 Can be any of the types supported by the ``collections`` argument
1279 to butler construction.
1281 Returns
1282 -------
1283 bounds : `dict`
1284 Dictionary mapping detector ID to a list of unit vectors on the
1285 sphere representing that detector's corners projected onto the sky.
1286 """
1287 if collections is None:
1288 collections = list(self.butler.collections.defaults)
1289 camera, versioned = loadCamera(self.butler, exposure.dataId, collections=collections)
1290 if not versioned and self.config.requireVersionedCamera:
1291 raise LookupError(f"No versioned camera found for exposure {exposure.dataId}.")
1293 orientation = lsst.geom.Angle(exposure.sky_angle, lsst.geom.degrees)
1294 radec = lsst.geom.SpherePoint(
1295 lsst.geom.Angle(exposure.tracking_ra, lsst.geom.degrees),
1296 lsst.geom.Angle(exposure.tracking_dec, lsst.geom.degrees),
1297 )
1299 if self.config.detectorId is None:
1300 detectorId = next(camera.getIdIter())
1301 else:
1302 detectorId = self.config.detectorId
1303 wcsDetector = camera[detectorId]
1305 # Ask the raw formatter to create the relevant WCS
1306 # This allows flips to be taken into account
1307 instrument = self.getInstrument(exposure.instrument)
1308 rawFormatter = instrument.getRawFormatter({"detector": detectorId})
1310 try:
1311 wcs = rawFormatter.makeRawSkyWcsFromBoresight(radec, orientation, wcsDetector) # type: ignore
1312 except AttributeError:
1313 raise TypeError(
1314 f"Raw formatter is {get_full_type_name(rawFormatter)} but visit"
1315 " definition requires it to support 'makeRawSkyWcsFromBoresight'"
1316 ) from None
1318 fpToSky = wcsDetector.getTransform(FOCAL_PLANE, PIXELS).then(wcs.getTransform())
1319 bounds = {}
1320 for detector in camera:
1321 pixelsToSky = detector.getTransform(PIXELS, FOCAL_PLANE).then(fpToSky)
1322 pixCorners = Box2D(detector.getBBox().dilatedBy(self.config.padding)).getCorners()
1323 bounds[detector.getId()] = [
1324 skyCorner.getVector() for skyCorner in pixelsToSky.applyForward(pixCorners)
1325 ]
1326 return bounds
1328 def compute(
1329 self, visit: VisitDefinitionData, *, collections: Any = None
1330 ) -> tuple[Region, dict[int, Region]]:
1331 # Docstring inherited from ComputeVisitRegionsTask.
1332 if self.config.mergeExposures:
1333 detectorBounds: dict[int, list[UnitVector3d]] = defaultdict(list)
1334 for exposure in visit.exposures:
1335 exposureDetectorBounds = self.computeExposureBounds(exposure, collections=collections)
1336 for detectorId, bounds in exposureDetectorBounds.items():
1337 detectorBounds[detectorId].extend(bounds)
1338 else:
1339 detectorBounds = self.computeExposureBounds(visit.exposures[0], collections=collections)
1340 visitBounds = []
1341 detectorRegions = {}
1342 for detectorId, bounds in detectorBounds.items():
1343 detectorRegions[detectorId] = ConvexPolygon.convexHull(bounds)
1344 visitBounds.extend(bounds)
1345 return ConvexPolygon.convexHull(visitBounds), detectorRegions