Repository: cta-observatory/ctapipe
Pull request #2409 (2023-10-15 19:45:42)

Interpolate pointing in TableLoader and HDF5EventSource
This adds the possibility to interpolate pointing information from a not yet used by ctapipe itself `/dl0/monitoring/telescope/pointing/<tel_id>` group.

This is critical for analysis of e.g.  LST observations.

At the moment, I add this table to the DL1 files output by ctapipe-process with a custom script, but eventually we should also add APIs for the datawriter / ctapipe-process / event sources to read this from monitoring files and put it into the outputs.


kosack: This can also be used eventually to fill in the event-wise monitoring info in the event loop as well not just in TableLoader.  

And just to note, this eventually must also interpolate other mis-pointing parameters like camera-z-axis rotation and any other deformation that needs to be applied during the DL1 step.  LST doesn't compute those yet I guess, but it's good not to forget they will be necessary eventually!
codecov[bot]: ## [Codecov](https://app.codecov.io/gh/cta-observatory/ctapipe/pull/2409?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) Report
All modified and coverable lines are covered by tests :white_check_mark:
> Comparison is base [(`49ab82b`)](https://app.codecov.io/gh/cta-observatory/ctapipe/commit/49ab82b84c7eb2b41051827df100fa4489d1065e?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) 60.60% compared to head [(`2f16e3b`)](https://app.codecov.io/gh/cta-observatory/ctapipe/pull/2409?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory) 60.60%.


<details><summary>Additional details and impacted files</summary>


```diff
@@           Coverage Diff           @@
##             main    #2409   +/-   ##
=======================================
  Coverage   60.60%   60.60%           
=======================================
  Files           3        3           
  Lines          33       33           
=======================================
  Hits           20       20           
  Misses         13       13           
```



</details>

[:umbrella: View full report in Codecov by Sentry](https://app.codecov.io/gh/cta-observatory/ctapipe/pull/2409?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).   
:loudspeaker: Have feedback on the report? [Share it here](https://about.codecov.io/codecov-pr-comment-feedback/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=cta-observatory).


Patch:
diff --git a/docs/changes/2409.feature.rst b/docs/changes/2409.feature.rst
new file mode 100644
index 00000000..7d5af104
--- /dev/null
+++ b/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.
diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py
index c270a612..a451bde4 100644
--- a/src/ctapipe/io/hdf5eventsource.py
+++ b/src/ctapipe/io/hdf5eventsource.py
@@ -52,7 +52,8 @@ from .astropy_helpers import read_table
 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 @@ class HDF5EventSource(EventSource):
             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 @@ class HDF5EventSource(EventSource):
                 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 @@ class HDF5EventSource(EventSource):
         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:
diff --git a/src/ctapipe/io/pointing.py b/src/ctapipe/io/pointing.py
new file mode 100644
index 00000000..412c0510
--- /dev/null
+++ b/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
diff --git a/src/ctapipe/io/tableloader.py b/src/ctapipe/io/tableloader.py
index 74a39030..36a954fb 100644
--- a/src/ctapipe/io/tableloader.py
+++ b/src/ctapipe/io/tableloader.py
@@ -11,11 +11,10 @@ import tables
 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 @@ TRUE_IMPACT_GROUP = "/simulation/event/telescope/impact"
 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 @@ class TableLoader(Component):
 
         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 @@ class TableLoader(Component):
             )
             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
 
diff --git a/src/ctapipe/reco/sklearn.py b/src/ctapipe/reco/sklearn.py
index 24c565c4..b9005a46 100644
--- a/src/ctapipe/reco/sklearn.py
+++ b/src/ctapipe/reco/sklearn.py
@@ -823,13 +823,20 @@ class DispReconstructor(Reconstructor):
         fov_lon = table["hillas_fov_lon"].quantity + disp * np.cos(psi)
         fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi)
 
-        # FIXME: Assume constant and parallel pointing for each run
-        self.log.warning("Assuming constant and parallel pointing for each run")
+        # prefer to use pointing interpolated to event
+        if "telescope_pointing_altitude" in table.colnames:
+            pointing_alt = table["telescope_pointing_altitude"]
+            pointing_az = table["telescope_pointing_azimuth"]
+        else:
+            # fallback to fixed pointing of ob
+            pointing_alt = table["subarray_pointing_lat"]
+            pointing_az = table["subarray_pointing_lon"]
+
         alt, az = telescope_to_horizontal(
             lon=fov_lon,
             lat=fov_lat,
-            pointing_alt=table["subarray_pointing_lat"],
-            pointing_az=table["subarray_pointing_lon"],
+            pointing_alt=pointing_alt,
+            pointing_az=pointing_az,
         )
 
         altaz_result = Table(


Test patch:
diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py
index 42d53fc0..4d743e5d 100644
--- a/src/ctapipe/conftest.py
+++ b/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
diff --git a/src/ctapipe/io/tests/test_hdf5eventsource.py b/src/ctapipe/io/tests/test_hdf5eventsource.py
index f426e024..bc8b1aa3 100644
--- a/src/ctapipe/io/tests/test_hdf5eventsource.py
+++ b/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)
diff --git a/src/ctapipe/io/tests/test_pointing.py b/src/ctapipe/io/tests/test_pointing.py
new file mode 100644
index 00000000..10913053
--- /dev/null
+++ b/src/ctapipe/io/tests/test_pointing.py
@@ -0,0 +1,154 @@
+import astropy.units as u
+import numpy as np
+import pytest
+import tables
+from astropy.table import Table
+from astropy.time import Time
+
+
+def test_simple():
+    """Test pointing interpolation"""
+    from ctapipe.io.pointing import PointingInterpolator
+
+    t0 = Time("2022-01-01T00:00:00")
+
+    table = Table(
+        {
+            "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s,
+            "azimuth": np.linspace(0.0, 10.0, 6) * u.deg,
+            "altitude": np.linspace(70.0, 60.0, 6) * u.deg,
+        },
+    )
+
+    interpolator = PointingInterpolator()
+    interpolator.add_table(1, table)
+
+    alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s)
+    assert u.isclose(alt, 69 * u.deg)
+    assert u.isclose(az, 1 * u.deg)
+
+    with pytest.raises(KeyError):
+        interpolator(tel_id=2, time=t0 + 1 * u.s)
+
+
+def test_azimuth_switchover():
+    """Test pointing interpolation"""
+    from ctapipe.io.pointing import PointingInterpolator
+
+    t0 = Time("2022-01-01T00:00:00")
+
+    table = Table(
+        {
+            "time": t0 + [0, 1, 2] * u.s,
+            "azimuth": [359, 1, 3] * u.deg,
+            "altitude": [60, 61, 62] * u.deg,
+        },
+    )
+
+    interpolator = PointingInterpolator()
+    interpolator.add_table(1, table)
+
+    alt, az = interpolator(tel_id=1, time=t0 + 0.5 * u.s)
+    assert u.isclose(az, 360 * u.deg)
+    assert u.isclose(alt, 60.5 * u.deg)
+
+
+def test_invalid_input():
+    """Test invalid pointing tables raise nice errors"""
+    from ctapipe.io.pointing import PointingInterpolator
+
+    wrong_time = Table(
+        {
+            "time": [1, 2, 3] * u.s,
+            "azimuth": [1, 2, 3] * u.deg,
+            "altitude": [1, 2, 3] * u.deg,
+        }
+    )
+
+    interpolator = PointingInterpolator()
+    with pytest.raises(TypeError, match="astropy.time.Time"):
+        interpolator.add_table(1, wrong_time)
+
+    wrong_unit = Table(
+        {
+            "time": Time(1.7e9 + np.arange(3), format="unix"),
+            "azimuth": [1, 2, 3] * u.m,
+            "altitude": [1, 2, 3] * u.deg,
+        }
+    )
+    with pytest.raises(ValueError, match="compatible with 'rad'"):
+        interpolator.add_table(1, wrong_unit)
+
+    wrong_unit = Table(
+        {
+            "time": Time(1.7e9 + np.arange(3), format="unix"),
+            "azimuth": [1, 2, 3] * u.deg,
+            "altitude": [1, 2, 3],
+        }
+    )
+    with pytest.raises(ValueError, match="compatible with 'rad'"):
+        interpolator.add_table(1, wrong_unit)
+
+
+def test_hdf5(tmp_path):
+    from ctapipe.io import write_table
+    from ctapipe.io.pointing import PointingInterpolator
+
+    t0 = Time("2022-01-01T00:00:00")
+
+    table = Table(
+        {
+            "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s,
+            "azimuth": np.linspace(0.0, 10.0, 6) * u.deg,
+            "altitude": np.linspace(70.0, 60.0, 6) * u.deg,
+        },
+    )
+
+    path = tmp_path / "pointing.h5"
+    write_table(table, path, "/dl0/monitoring/telescope/pointing/tel_001")
+    with tables.open_file(path) as h5file:
+        interpolator = PointingInterpolator(h5file)
+        alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s)
+        assert u.isclose(alt, 69 * u.deg)
+        assert u.isclose(az, 1 * u.deg)
+
+
+def test_bounds():
+    """Test invalid pointing tables raise nice errors"""
+    from ctapipe.io.pointing import PointingInterpolator
+
+    t0 = Time("2022-01-01T00:00:00")
+
+    table = Table(
+        {
+            "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s,
+            "azimuth": np.linspace(0.0, 10.0, 6) * u.deg,
+            "altitude": np.linspace(70.0, 60.0, 6) * u.deg,
+        },
+    )
+
+    interpolator = PointingInterpolator()
+    interpolator.add_table(1, table)
+
+    with pytest.raises(ValueError, match="below the interpolation range"):
+        interpolator(tel_id=1, time=t0 - 0.1 * u.s)
+
+    with pytest.raises(ValueError, match="above the interpolation range"):
+        interpolator(tel_id=1, time=t0 + 10.2 * u.s)
+
+    interpolator = PointingInterpolator(bounds_error=False)
+    interpolator.add_table(1, table)
+    for dt in (-0.1, 10.1) * u.s:
+        alt, az = interpolator(tel_id=1, time=t0 + dt)
+        assert np.isnan(alt.value)
+        assert np.isnan(az.value)
+
+    interpolator = PointingInterpolator(bounds_error=False, extrapolate=True)
+    interpolator.add_table(1, table)
+    alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s)
+    assert u.isclose(alt, 71 * u.deg)
+    assert u.isclose(az, -1 * u.deg)
+
+    alt, az = interpolator(tel_id=1, time=t0 + 11 * u.s)
+    assert u.isclose(alt, 59 * u.deg)
+    assert u.isclose(az, 11 * u.deg)
diff --git a/src/ctapipe/io/tests/test_table_loader.py b/src/ctapipe/io/tests/test_table_loader.py
index adfd59b3..398b58d4 100644
--- a/src/ctapipe/io/tests/test_table_loader.py
+++ b/src/ctapipe/io/tests/test_table_loader.py
@@ -416,3 +416,13 @@ 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_mon_pointing_file, tmp_path):
+    from ctapipe.io import TableLoader
+
+    with TableLoader(dl1_mon_pointing_file, 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