Restore reading of pointing for files written with ctapipe<0.21, fixes …

…#2651

src/ctapipe/io/hdf5eventsource.py

@@ -49,6 +49,7 @@
 from ..instrument import SubarrayDescription
 from ..instrument.optics import FocalLengthKind
 from ..monitoring.interpolation import PointingInterpolator
+from ..utils import IndexFinder
 from .astropy_helpers import read_table
 from .datalevels import DataLevel
 from .eventsource import EventSource
@@ -230,13 +231,29 @@ def __init__(self, input_url=None, config=None, parent=None, **kwargs):
             self.file_.root.configuration.observation.observation_block.col("obs_id")
         )
         pointing_key = "/configuration/telescope/pointing"
+        # for ctapipe <0.21
+        legacy_pointing_key = "/dl1/monitoring/telescope/pointing"
+        self._legacy_tel_pointing_finders = {}
+        self._legacy_tel_pointing_tables = {}
+
         self._constant_telescope_pointing = {}
         if pointing_key in self.file_.root:
             for h5table in self.file_.root[pointing_key]._f_iter_nodes("Table"):
                 tel_id = int(h5table._v_name.partition("tel_")[-1])
                 table = QTable(read_table(self.file_, h5table._v_pathname), copy=False)
                 table.add_index("obs_id")
                 self._constant_telescope_pointing[tel_id] = table
+        elif legacy_pointing_key in self.file_.root:
+            self.log.info(
+                "Found file written using ctapipe<0.21, using legacy pointing information"
+            )
+            for node in self.file_.root[legacy_pointing_key]:
+                tel_id = int(node.name.removeprefix("tel_"))
+                table = QTable(read_table(self.file_, node._v_pathname), copy=False)
+                self._legacy_tel_pointing_tables[tel_id] = table
+                self._legacy_tel_pointing_finders[tel_id] = IndexFinder(
+                    table["time"].mjd
+                )
 
         self._simulated_shower_distributions = (
             self._read_simulated_shower_distributions()
@@ -767,12 +784,18 @@ def _fill_telescope_pointing(self, event, tel_pointing_interpolator=None):
             return
 
         for tel_id in event.trigger.tels_with_trigger:
-            tel_pointing = self._constant_telescope_pointing.get(tel_id)
-            if tel_pointing is None:
-                continue
-
-            current = tel_pointing.loc[event.index.obs_id]
-            event.pointing.tel[tel_id] = TelescopePointingContainer(
-                altitude=current["telescope_pointing_altitude"],
-                azimuth=current["telescope_pointing_azimuth"],
-            )
+            if (
+                tel_pointing := self._constant_telescope_pointing.get(tel_id)
+            ) is not None:
+                current = tel_pointing.loc[event.index.obs_id]
+                event.pointing.tel[tel_id] = TelescopePointingContainer(
+                    altitude=current["telescope_pointing_altitude"],
+                    azimuth=current["telescope_pointing_azimuth"],
+                )
+            elif (finder := self._legacy_tel_pointing_finders.get(tel_id)) is not None:
+                index = finder.closest(event.trigger.time.mjd)
+                row = self._legacy_tel_pointing_tables[tel_id][index]
+                event.pointing.tel[tel_id] = TelescopePointingContainer(
+                    altitude=row["altitude"],
+                    azimuth=row["azimuth"],
+                )

src/ctapipe/io/tests/test_hdf5eventsource.py

@@ -285,3 +285,17 @@ def test_provenance(dl1_file, provenance):
     assert inputs[0]["url"] == str(dl1_file)
     meta = _read_reference_metadata_hdf5(dl1_file)
     assert inputs[0]["reference_meta"].product.id_ == meta.product.id_
+
+
+def test_pointing_old_file():
+    input_url = "dataset://gamma_diffuse_dl2_train_small.dl2.h5"
+
+    n_read = 0
+    with HDF5EventSource(input_url, max_events=5) as source:
+        for e in source:
+            assert e.pointing.tel.keys() == set(e.trigger.tels_with_trigger)
+            for pointing in e.pointing.tel.values():
+                assert u.isclose(pointing.altitude, 70 * u.deg)
+                assert u.isclose(pointing.azimuth, 0 * u.deg)
+            n_read += 1
+    assert n_read == 5

src/ctapipe/tools/tests/test_process.py

@@ -499,7 +499,7 @@ def test_only_trigger_and_simulation(tmp_path):
     [
         pytest.param(
             "dataset://gamma_diffuse_dl2_train_small.dl2.h5",
-            ["--no-write-images", "--max-events=5"],
+            ["--no-write-images", "--max-events=20"],
             id="0.17",
         )
     ],
@@ -514,9 +514,24 @@ def test_on_old_file(input_url, args, tmp_path):
             f"--config={config}",
             f"--input={input_url}",
             f"--output={output_path}",
+            "--write-showers",
             "--overwrite",
             *args,
         ],
         cwd=tmp_path,
         raises=True,
     )
+
+    with tables.open_file(output_path) as f:
+        assert "/configuration/telescope/pointing" in f.root
+
+    with TableLoader(output_path) as loader:
+        events = loader.read_subarray_events()
+
+        # check that we have valid reconstructions and that in case
+        # we don't, is_valid is False, regression test for #2651
+        finite_reco = np.isfinite(events["HillasReconstructor_alt"])
+        assert np.any(finite_reco)
+        np.testing.assert_array_equal(
+            finite_reco, events["HillasReconstructor_is_valid"]
+        )