Merge pull request #2409 from cta-observatory/interpolate_pointing_ta…

…bleloader

Interpolate pointing in TableLoader and HDF5EventSource

docs/changes/2409.feature.rst

@@ -0,0 +1,5 @@
+Add support for interpolating a monitoring pointing table
+in ``TableLoader``. The corresponding table is not yet written by ``ctapipe``,
+but can be written by external tools.
+This is to enable analysis of real observations, where the pointing changes over time in
+alt/az.

src/ctapipe/conftest.py

@@ -1,12 +1,15 @@
 """
 common pytest fixtures for tests in ctapipe
 """
-
+import shutil
 from copy import deepcopy
 
 import astropy.units as u
+import numpy as np
 import pytest
+import tables
 from astropy.coordinates import EarthLocation
+from astropy.table import Table
 from pytest_astropy_header.display import PYTEST_HEADER_MODULES
 
 from ctapipe.core import run_tool
@@ -646,3 +649,38 @@ def disp_reconstructor_path(model_tmp_path, gamma_train_clf):
 def reference_location():
     """a dummy EarthLocation to use for SubarrayDescriptions"""
     return EarthLocation(lon=-17 * u.deg, lat=28 * u.deg, height=2200 * u.m)
+
+
+@pytest.fixture(scope="session")
+def dl1_mon_pointing_file(dl1_file, dl1_tmp_path):
+    from ctapipe.instrument import SubarrayDescription
+    from ctapipe.io import read_table, write_table
+
+    path = dl1_tmp_path / "dl1_mon_ponting.dl1.h5"
+    shutil.copy(dl1_file, path)
+
+    events = read_table(path, "/dl1/event/subarray/trigger")
+    subarray = SubarrayDescription.from_hdf(path)
+
+    # 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}")
+
+    # remove static pointing table
+    with tables.open_file(path, "r+") as f:
+        f.remove_node("/configuration/telescope/pointing", recursive=True)
+
+    return path

src/ctapipe/io/hdf5eventsource.py

@@ -52,7 +52,8 @@
 from .datalevels import DataLevel
 from .eventsource import EventSource
 from .hdf5tableio import HDF5TableReader
-from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP
+from .pointing import PointingInterpolator
+from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP, POINTING_GROUP
 
 __all__ = ["HDF5EventSource"]
 
@@ -581,6 +582,13 @@ def _generator(self):
             ignore_columns={"trigger_pixels"},
         )
 
+        pointing_interpolator = None
+        if POINTING_GROUP in self.file_.root:
+            pointing_interpolator = PointingInterpolator(
+                h5file=self.file_,
+                parent=self,
+            )
+
         counter = 0
         for trigger, index in events:
             data = ArrayEventContainer(
@@ -630,7 +638,7 @@ def _generator(self):
                 continue
 
             self._fill_array_pointing(data)
-            self._fill_telescope_pointing(data)
+            self._fill_telescope_pointing(data, pointing_interpolator)
 
             for tel_id in data.trigger.tel.keys():
                 key = f"tel_{tel_id:03d}"
@@ -715,10 +723,19 @@ def _fill_array_pointing(self, event):
         else:
             raise ValueError(f"Unsupported pointing frame: {frame}")
 
-    def _fill_telescope_pointing(self, event):
+    def _fill_telescope_pointing(self, event, tel_pointing_interpolator=None):
         """
         Fill the telescope pointing information of a given event
         """
+        if tel_pointing_interpolator is not None:
+            for tel_id, trigger in event.trigger.tel.items():
+                alt, az = tel_pointing_interpolator(tel_id, trigger.time)
+                event.pointing.tel[tel_id] = TelescopePointingContainer(
+                    altitude=alt,
+                    azimuth=az,
+                )
+            return
+
         for tel_id in event.trigger.tels_with_trigger:
             tel_pointing = self._constant_telescope_pointing.get(tel_id)
             if tel_pointing is None:

src/ctapipe/io/pointing.py

@@ -0,0 +1,137 @@
+from typing import Any
+
+import astropy.units as u
+import numpy as np
+import tables
+from astropy.time import Time
+from scipy.interpolate import interp1d
+
+from ctapipe.core import Component, traits
+
+from .astropy_helpers import read_table
+
+
+class PointingInterpolator(Component):
+    """
+    Interpolate pointing from a monitoring table to a given timestamp.
+
+    Parameters
+    ----------
+    h5file : None | tables.File
+        A open hdf5 file with read access.
+        The monitoring table is expected to be stored in that file at
+        ``/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}``
+
+        If not given, monitoring tables can be added via `PointingInterpolator.add_table`.
+    """
+
+    bounds_error = traits.Bool(
+        default_value=True,
+        help="If true, raises an exception when trying to extrapolate out of the given table",
+    ).tag(config=True)
+
+    extrapolate = traits.Bool(
+        help="If bounds_error is False, this flag will specify whether values outside"
+        "the available values are filled with nan (False) or extrapolated (True).",
+        default_value=False,
+    ).tag(config=True)
+
+    def __init__(self, h5file=None, **kwargs):
+        super().__init__(**kwargs)
+
+        if h5file is not None and 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 add_table(self, tel_id, pointing_table):
+        """
+        Add a table to this interpolator
+
+        Parameters
+        ----------
+        tel_id : int
+            Telescope id
+        pointing_table : astropy.table.Table
+            Table of pointing values, expected columns
+            are ``time`` as ``Time`` column, ``azimuth`` and ``altitude``
+            as quantity columns.
+        """
+        missing = {"time", "azimuth", "altitude"} - set(pointing_table.colnames)
+        if len(missing) > 0:
+            raise ValueError(f"Table is missing required column(s): {missing}")
+
+        if not isinstance(pointing_table["time"], Time):
+            raise TypeError("'time' column of pointing table must be astropy.time.Time")
+
+        for col in ("azimuth", "altitude"):
+            unit = pointing_table[col].unit
+            if unit is None or not u.rad.is_equivalent(unit):
+                raise ValueError(f"{col} must have units compatible with 'rad'")
+
+        # sort first, so it's not done twice for each interpolator
+        pointing_table.sort("time")
+        # interpolate in mjd TAI. Float64 mjd is precise enough for pointing
+        # and TAI is contiguous, so no issues with leap seconds.
+        mjd = pointing_table["time"].tai.mjd
+
+        az = pointing_table["azimuth"].quantity.to_value(u.rad)
+        # prepare azimuth for interpolation by "unwrapping": i.e. turning
+        # [359, 1] into [359, 361]. This assumes that if we get values like
+        # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees
+        # the other way around. This should be true for all telescopes given
+        # the sampling speed of pointing values and their maximum movement speed.
+        # No telescope can turn more than 180° in 2 seconds.
+        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 _read_pointing_table(self, tel_id):
+        pointing_table = read_table(
+            self.h5file,
+            f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}",
+        )
+        self.add_table(tel_id, pointing_table)
+
+    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:
+            if self.h5file is not None:
+                self._read_pointing_table(tel_id)
+            else:
+                raise KeyError(f"No pointing table available for tel_id {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

src/ctapipe/io/tableloader.py

@@ -11,11 +11,10 @@
 from astropy.table import Table, hstack, vstack
 from astropy.utils.decorators import lazyproperty
 
-from ctapipe.instrument.optics import FocalLengthKind
-
 from ..core import Component, Provenance, traits
-from ..instrument import SubarrayDescription
+from ..instrument import FocalLengthKind, SubarrayDescription
 from .astropy_helpers import join_allow_empty, read_table
+from .pointing import PointingInterpolator
 
 __all__ = ["TableLoader"]
 
@@ -30,7 +29,8 @@
 SIMULATION_CONFIG_TABLE = "/configuration/simulation/run"
 SHOWER_DISTRIBUTION_TABLE = "/simulation/service/shower_distribution"
 OBSERVATION_TABLE = "/configuration/observation/observation_block"
-FIXED_POINTING_GROUP = "configuration/telescope/pointing"
+FIXED_POINTING_GROUP = "/configuration/telescope/pointing"
+POINTING_GROUP = "/dl0/monitoring/telescope/pointing"
 
 DL2_SUBARRAY_GROUP = "/dl2/event/subarray"
 DL2_TELESCOPE_GROUP = "/dl2/event/telescope"
@@ -204,7 +204,7 @@ class TableLoader(Component):
 
     pointing = traits.Bool(
         True,
-        help="Load pointing information and join to events",
+        help="Load pointing information and interpolate / join to events",
     ).tag(config=True)
 
     focal_length_choice = traits.UseEnum(
@@ -248,6 +248,11 @@ def __init__(self, input_url=None, h5file=None, **kwargs):
 
         Provenance().add_input_file(self.input_url, role="Event data")
 
+        self._pointing_interpolator = PointingInterpolator(
+            h5file=self.h5file,
+            parent=self,
+        )
+
     def close(self):
         """Close the underlying hdf5 file"""
         if self._should_close:
@@ -595,15 +600,21 @@ def _read_telescope_events_for_id(
             )
             table = _join_telescope_events(table, impacts)
 
-        if len(table) > 0 and pointing and FIXED_POINTING_GROUP in self.h5file.root:
-            pointing = read_table(
-                self.h5file, f"{FIXED_POINTING_GROUP}/tel_{tel_id:03d}"
-            )
-            # we know that we only have the tel_id we are looking for, remove optimize joining
-            del pointing["tel_id"]
-            table = join_allow_empty(
-                table, pointing, ["obs_id"], "left", keep_order=True
-            )
+        if len(table) > 0 and pointing:
+            # prefer monitoring pointing
+            if 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
+            elif FIXED_POINTING_GROUP in self.h5file.root:
+                pointing_table = read_table(
+                    self.h5file, f"{FIXED_POINTING_GROUP}/tel_{tel_id:03d}"
+                )
+                # we know that we only have the tel_id we are looking for, remove optimize joining
+                del pointing_table["tel_id"]
+                table = join_allow_empty(
+                    table, pointing_table, ["obs_id"], "left", keep_order=True
+                )
 
         return table
 

src/ctapipe/io/tests/test_hdf5eventsource.py

@@ -251,3 +251,14 @@ def test_dl1_camera_frame(dl1_camera_frame_file):
                 assert sim.true_parameters.hillas.intensity is not None
 
         assert tel_id is not None, "did not test any events"
+
+
+def test_interpolate_pointing(dl1_mon_pointing_file):
+    from ctapipe.io import HDF5EventSource
+
+    with HDF5EventSource(dl1_mon_pointing_file) as source:
+        for e in source:
+            assert set(e.pointing.tel.keys()) == set(e.trigger.tels_with_trigger)
+            for pointing in e.pointing.tel.values():
+                assert not np.isnan(pointing.azimuth)
+                assert not np.isnan(pointing.altitude)