Acada lst integration tests

ctapipe/calib/camera/calibrator.py

@@ -123,9 +123,6 @@ def __init__(
         super().__init__(subarray=subarray, config=config, parent=parent, **kwargs)
         self.subarray = subarray
 
-        self._r1_empty_warn = False
-        self._dl0_empty_warn = False
-
         self.image_extractors = {}
 
         if image_extractor is None:
@@ -155,24 +152,12 @@ def __init__(
 
     def _check_r1_empty(self, waveforms):
         if waveforms is None:
-            if not self._r1_empty_warn:
-                warnings.warn(
-                    "Encountered an event with no R1 data. "
-                    "DL0 is unchanged in this circumstance."
-                )
-                self._r1_empty_warn = True
             return True
         else:
             return False
 
     def _check_dl0_empty(self, waveforms):
         if waveforms is None:
-            if not self._dl0_empty_warn:
-                warnings.warn(
-                    "Encountered an event with no DL0 data. "
-                    "DL1 is unchanged in this circumstance."
-                )
-                self._dl0_empty_warn = True
             return True
         else:
             return False

ctapipe/core/telescope_component.py

@@ -181,13 +181,6 @@ def attach_subarray(self, subarray):
                 ]
                 logger.debug(f"argument '{arg}' matched: {matched_tel_types}")
 
-                if len(matched_tel_types) == 0:
-                    logger.warning(
-                        "TelescopeParameter type argument '%s' did not match "
-                        "any known telescope types",
-                        arg,
-                    )
-
                 for tel_type in matched_tel_types:
                     self._value_for_type[tel_type] = value
                     for tel_id in subarray.get_tel_ids_for_type(tel_type):

ctapipe/core/tool.py

@@ -176,6 +176,11 @@ def main():
 
     quiet = Bool(default_value=False).tag(config=True)
     overwrite = Bool(default_value=False).tag(config=True)
+    debug = Bool(
+        default_value=False,
+        help="If true, exceptions are not caught and converted"
+        " to exit codes to enable using a debugger",
+    ).tag(config=True)
 
     _log_formatter_cls = ColoredFormatter
 
@@ -209,6 +214,11 @@ def __init__(self, **kwargs):
                 {"Tool": {"overwrite": True}},
                 "Overwrite existing output files without asking",
             ),
+            "debug": (
+                {"Tool": {"debug": True}},
+                "Disable the transformation of exceptions"
+                " to exit codes to enable debugger usage",
+            ),
         }
         self.flags = {**flags, **self.flags}
 
@@ -383,7 +393,7 @@ def finish(self):
         """
         self.log.info("Goodbye")
 
-    def run(self, argv=None, raises=False):
+    def run(self, argv=None, raises=None):
         """Run the tool.
 
         This automatically calls `Tool.initialize`, `Tool.start` and `Tool.finish`
@@ -410,6 +420,8 @@ def run(self, argv=None, raises=False):
                 Provenance().start_activity(self.name)
 
                 self.initialize(argv)
+                if raises is None:
+                    raises = self.debug
 
                 self.setup()
                 self.is_setup = True

ctapipe/io/hdf5merger.py

@@ -158,6 +158,15 @@ class HDF5Merger(Component):
         True, help="Whether to include processing statistics in merged output"
     ).tag(config=True)
 
+    single_ob = traits.Bool(
+        False,
+        help=(
+            "If true, input files are assumed to be multiple chunks from the same"
+            " observation block and the ob / sb blocks will only be copied from "
+            " the first input file"
+        ),
+    ).tag(config=True)
+
     def __init__(self, output_path=None, **kwargs):
         # enable using output_path as posarg
         if output_path not in {None, traits.Undefined}:
@@ -202,6 +211,7 @@ def __init__(self, output_path=None, **kwargs):
                 focal_length_choice=FocalLengthKind.EQUIVALENT,
             )
             self.required_nodes = _get_required_nodes(self.h5file)
+        self._n_merged = 0
 
     def __call__(self, other: Union[str, Path, tables.File]):
         """
@@ -213,7 +223,7 @@ def __call__(self, other: Union[str, Path, tables.File]):
 
         with exit_stack:
             # first file to be merged
-            if self.meta is None:
+            if self._n_merged == 0:
                 self.meta = self._read_meta(other)
                 self.data_model_version = self.meta.product.data_model_version
                 metadata.write_to_hdf5(self.meta.to_dict(), self.h5file)
@@ -229,6 +239,7 @@ def __call__(self, other: Union[str, Path, tables.File]):
                     self.log.info(
                         "Updated required nodes to %s", sorted(self.required_nodes)
                     )
+                self._n_merged += 1
             finally:
                 self._update_meta()
 
@@ -272,13 +283,15 @@ def _append(self, other):
         # Configuration
         self._append_subarray(other)
 
-        config_keys = [
-            "/configuration/observation/scheduling_block",
-            "/configuration/observation/observation_block",
-        ]
-        for key in config_keys:
-            if key in other.root:
-                self._append_table(other, other.root[key])
+        # in case of "single_ob", we only copy sb/ob blocks for the first file
+        if not self.single_ob or self._n_merged == 0:
+            config_keys = [
+                "/configuration/observation/scheduling_block",
+                "/configuration/observation/observation_block",
+            ]
+            for key in config_keys:
+                if key in other.root:
+                    self._append_table(other, other.root[key])
 
         # Simulation
         simulation_table_keys = [

ctapipe/io/pointing.py

@@ -0,0 +1,81 @@
+from typing import Any
+
+import astropy.units as u
+import numpy as np
+import tables
+from scipy.interpolate import interp1d
+
+from ctapipe.core import Component, traits
+
+from .astropy_helpers import read_table
+
+
+class PointingInterpolator(Component):
+    bounds_error = traits.Bool(default_value=True).tag(config=True)
+    extrapolate = traits.Bool(default_value=False).tag(config=True)
+
+    def __init__(self, h5file, **kwargs):
+        super().__init__(**kwargs)
+
+        if not isinstance(h5file, tables.File):
+            raise TypeError("h5file must be a tables.File")
+        self.h5file = h5file
+
+        self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False)
+        if self.bounds_error:
+            self.interp_options["bounds_error"] = True
+        elif self.extrapolate:
+            self.interp_options["bounds_error"] = False
+            self.interp_options["fill_value"] = "extrapolate"
+        else:
+            self.interp_options["bounds_error"] = False
+            self.interp_options["fill_value"] = np.nan
+
+        self._alt_interpolators = {}
+        self._az_interpolators = {}
+
+    def _read_pointing_table(self, tel_id):
+        pointing_table = read_table(
+            self.h5file,
+            f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}",
+        )
+
+        # sort first, so it's not done twice for each interpolator
+        pointing_table.sort("time")
+        mjd = pointing_table["time"].tai.mjd
+
+        az = pointing_table["azimuth"].quantity.to_value(u.rad)
+        # prepare azimuth for interpolation "unwrapping", i.e. turning 359, 1 into 359, 361
+        az = np.unwrap(az)
+        alt = pointing_table["altitude"].quantity.to_value(u.rad)
+
+        self._az_interpolators[tel_id] = interp1d(mjd, az, **self.interp_options)
+        self._alt_interpolators[tel_id] = interp1d(mjd, alt, **self.interp_options)
+
+    def __call__(self, tel_id, time):
+        """
+        Interpolate alt/az for given time and tel_id.
+
+        Parameters
+        ----------
+        tel_id : int
+            telescope id
+        time : astropy.time.Time
+            time for which to interpolate the pointing
+
+        Returns
+        -------
+        altitude : astropy.units.Quantity[deg]
+            interpolated altitude angle
+        azimuth : astropy.units.Quantity[deg]
+            interpolated azimuth angle
+
+
+        """
+        if tel_id not in self._az_interpolators:
+            self._read_pointing_table(tel_id)
+
+        mjd = time.tai.mjd
+        az = u.Quantity(self._az_interpolators[tel_id](mjd), u.rad, copy=False)
+        alt = u.Quantity(self._alt_interpolators[tel_id](mjd), u.rad, copy=False)
+        return alt, az

ctapipe/io/tableloader.py

@@ -13,6 +13,7 @@
 from astropy.utils.decorators import lazyproperty
 
 from ctapipe.instrument.optics import FocalLengthKind
+from ctapipe.io.pointing import PointingInterpolator
 
 from ..core import Component, Provenance, traits
 from ..instrument import SubarrayDescription
@@ -31,6 +32,7 @@
 SIMULATION_CONFIG_TABLE = "/configuration/simulation/run"
 SHOWER_DISTRIBUTION_TABLE = "/simulation/service/shower_distribution"
 OBSERVATION_TABLE = "/configuration/observation/observation_block"
+POINTING_GROUP = "/dl0/monitoring/telescope/pointing"
 
 DL2_SUBARRAY_GROUP = "/dl2/event/subarray"
 DL2_TELESCOPE_GROUP = "/dl2/event/telescope"
@@ -177,9 +179,11 @@ class TableLoader(Component):
     ).tag(config=True)
 
     load_dl1_images = traits.Bool(False, help="load extracted images").tag(config=True)
+
     load_dl1_parameters = traits.Bool(
         True, help="load reconstructed image parameters"
     ).tag(config=True)
+
     load_dl1_muons = traits.Bool(False, help="load muon ring parameters").tag(
         config=True
     )
@@ -206,6 +210,11 @@ class TableLoader(Component):
         False, help="join observation information to each event"
     ).tag(config=True)
 
+    interpolate_pointing = traits.Bool(
+        False,
+        help="If True, load monitoring pointing information and interpolate for each telescope event",
+    ).tag(config=True)
+
     focal_length_choice = traits.UseEnum(
         FocalLengthKind,
         default_value=FocalLengthKind.EFFECTIVE,
@@ -250,6 +259,13 @@ def __init__(self, input_url=None, h5file=None, **kwargs):
         self.instrument_table = None
         if self.load_instrument:
             self.instrument_table = self.subarray.to_table("joined")
+            # to prevent merging meta data onto event tables
+            self.instrument_table.meta.clear()
+
+        self._pointing_interpolator = PointingInterpolator(
+            h5file=self.h5file,
+            parent=self,
+        )
 
         groups = {
             "load_dl1_parameters": PARAMETERS_GROUP,
@@ -529,6 +545,11 @@ def _read_telescope_events_for_id(self, tel_id, start=None, stop=None):
             )
             table = _join_telescope_events(table, impacts)
 
+        if self.interpolate_pointing and POINTING_GROUP in self.h5file.root:
+            alt, az = self._pointing_interpolator(tel_id, table["time"])
+            table["telescope_pointing_altitude"] = alt
+            table["telescope_pointing_azimuth"] = az
+
         return table
 
     def _read_telescope_events_for_ids(self, tel_ids, start=None, stop=None):

ctapipe/io/tests/test_table_loader.py

@@ -1,3 +1,5 @@
+import shutil
+
 import astropy.units as u
 import numpy as np
 import pytest
@@ -436,3 +438,37 @@ def test_order_merged():
         for tel, table in tables.items():
             mask = np.isin(tel_trigger["tel_id"], loader.subarray.get_tel_ids(tel))
             check_equal_array_event_order(table, tel_trigger[mask])
+
+
+def test_interpolate_pointing(dl1_file, tmp_path):
+    from ctapipe.io import TableLoader, write_table
+
+    path = tmp_path / dl1_file.name
+    shutil.copy(dl1_file, path)
+
+    with TableLoader(dl1_file) as loader:
+        events = loader.read_telescope_events()
+        subarray = loader.subarray
+
+    # create some dummy monitoring data
+    time = events["time"]
+    start, stop = time[[0, -1]]
+    duration = (stop - start).to_value(u.s)
+
+    # start a bit before, go a bit longer
+    dt = np.arange(-1, duration + 2, 1) * u.s
+    time_mon = start + dt
+
+    alt = (69 + 2 * dt / dt[-1]) * u.deg
+    az = (180 + 5 * dt / dt[-1]) * u.deg
+
+    table = Table({"time": time_mon, "azimuth": az, "altitude": alt})
+
+    for tel_id in subarray.tel:
+        write_table(table, path, f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}")
+
+    with TableLoader(path, interpolate_pointing=True) as loader:
+        events = loader.read_telescope_events([4])
+        assert len(events) > 0
+        assert "telescope_pointing_azimuth" in events.colnames
+        assert "telescope_pointing_altitude" in events.colnames