Repository: cta-observatory/ctapipe
Pull request #2785 (2025-06-23 21:14:04)

Raise clear error in case of non-altaz fixed pointing in disp reconstructor


Patch:
diff --git a/docs/changes/2776.feature.rst b/docs/changes/2776.feature.rst
deleted file mode 100644
index 5d03ca49..00000000
--- a/docs/changes/2776.feature.rst
+++ /dev/null
@@ -1,2 +0,0 @@
-Use the HDF5Merger to add the observation block and scheduling block info
-to the DL1 files produced by ctapipe-calculate-pixel-statistics tool.
diff --git a/docs/changes/2785.bugfix.rst b/docs/changes/2785.bugfix.rst
new file mode 100644
index 00000000..31a01a36
--- /dev/null
+++ b/docs/changes/2785.bugfix.rst
@@ -0,0 +1,3 @@
+Raise a clear error message in the ``DispReconstructor``
+in case of non-supported coordinate frames for the pointing
+directions instead of silently assuming AltAz.
diff --git a/src/ctapipe/reco/sklearn.py b/src/ctapipe/reco/sklearn.py
index acc77e96..49425dd0 100644
--- a/src/ctapipe/reco/sklearn.py
+++ b/src/ctapipe/reco/sklearn.py
@@ -23,6 +23,7 @@ from ctapipe.exceptions import TooFewEvents
 
 from ..containers import (
     ArrayEventContainer,
+    CoordinateFrameType,
     DispContainer,
     ParticleClassificationContainer,
     ReconstructedEnergyContainer,
@@ -796,6 +797,7 @@ class DispReconstructor(Reconstructor):
         altaz_table : `~astropy.table.Table`
             Table with resulting predictions of horizontal coordinates
         """
+
         table = self.feature_generator(table, subarray=self.subarray)
 
         n_rows = len(table)
@@ -825,15 +827,7 @@ 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)
 
-        # 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"]
-
+        pointing_alt, pointing_az = self._get_pointing(table)
         alt, az = telescope_to_horizontal(
             lon=fov_lon,
             lat=fov_lat,
@@ -870,6 +864,31 @@ class DispReconstructor(Reconstructor):
                 disp.n_jobs = n_jobs.new
                 sign.n_jobs = n_jobs.new
 
+    def _get_pointing(self, table):
+        # prefer to use pointing interpolated to event
+        if "telescope_pointing_altitude" in table.colnames:
+            return (
+                table["telescope_pointing_altitude"].quantity,
+                table["telescope_pointing_azimuth"].quantity,
+            )
+
+        # fallback to fixed pointing of ob
+        if len(np.unique(table["subarray_pointing_frame"])) > 1:
+            msg = "Subarray pointing frame must be the same for all events"
+            raise NotImplementedError(msg)
+
+        # for now only allow fixed altaz, real data should have telescope monitoring
+        # pointing and simulations have fixed pointing in alt az
+        frame_type = CoordinateFrameType(table["subarray_pointing_frame"][0])
+        if frame_type is CoordinateFrameType.ALTAZ:
+            return (
+                table["subarray_pointing_lat"].quantity,
+                table["subarray_pointing_lon"].quantity,
+            )
+
+        msg = f"Only AltAz frame supported for fixed subarray pointing, got {frame_type.name}"
+        raise NotImplementedError(msg)
+
 
 class CrossValidator(Component):
     """Class to train sklearn based reconstructors in a cross validation."""
diff --git a/src/ctapipe/tools/calculate_pixel_stats.py b/src/ctapipe/tools/calculate_pixel_stats.py
index 5a5980b4..8b8c4f65 100644
--- a/src/ctapipe/tools/calculate_pixel_stats.py
+++ b/src/ctapipe/tools/calculate_pixel_stats.py
@@ -10,13 +10,15 @@ from astropy.table import vstack
 from ctapipe.core import Tool
 from ctapipe.core.tool import ToolConfigurationError
 from ctapipe.core.traits import (
+    Bool,
     CInt,
     Path,
     Set,
     Unicode,
     classes_with_traits,
 )
-from ctapipe.io import HDF5Merger, write_table
+from ctapipe.instrument import SubarrayDescription
+from ctapipe.io import write_table
 from ctapipe.io.tableloader import TableLoader
 from ctapipe.monitoring.calculator import PixelStatisticsCalculator
 
@@ -64,6 +66,8 @@ class PixelStatisticsCalculatorTool(Tool):
         help="Output filename", default_value=pathlib.Path("monitoring.h5")
     ).tag(config=True)
 
+    overwrite = Bool(help="Overwrite output file if it exists").tag(config=True)
+
     aliases = {
         ("i", "input_url"): "TableLoader.input_url",
         ("o", "output_path"): "PixelStatisticsCalculatorTool.output_path",
@@ -71,13 +75,9 @@ class PixelStatisticsCalculatorTool(Tool):
 
     flags = {
         "overwrite": (
-            {"HDF5Merger": {"overwrite": True}},
+            {"PixelStatisticsCalculatorTool": {"overwrite": True}},
             "Overwrite existing files",
         ),
-        "append": (
-            {"HDF5Merger": {"append": True}},
-            "Append to existing files",
-        ),
     }
 
     classes = [
@@ -91,32 +91,26 @@ class PixelStatisticsCalculatorTool(Tool):
         self.input_data = self.enter_context(
             TableLoader(
                 parent=self,
-                dl1_images=True,  # Ensure that dl1 images are read
             )
         )
-        # Copy selected tables from the input file to the output file
-        self.log.info(
-            "Copying selected data and metadata to output destination using the HDF5Merger component."
-        )
-        # Disable the copy of waveforms and images in the HDF5Merger
-        self.log.debug(
-            "Overwriting the default configuration of the HDF5Merger "
-            "component to disable the copy of waveforms and images."
-        )
-        with HDF5Merger(
-            parent=self,
-            output_path=self.output_path,
-            r0_waveforms=False,
-            r1_waveforms=False,
-            dl1_images=False,
-            true_images=False,
-        ) as merger:
-            merger(self.input_data.input_url)
+        # Check that the input and output files are not the same
+        if self.input_data.input_url == self.output_path:
+            raise ToolConfigurationError(
+                "Input and output files are same. Fix your configuration / cli arguments."
+            )
+        if "dl1_images" in self.input_data.config.TableLoader:
+            if not self.input_data.dl1_images:
+                raise ToolConfigurationError(
+                    "The TableLoader must read dl1 images. Set 'dl1_images' to True."
+                )
+        self.input_data.dl1_images = True
+        # Load the subarray description from the input file
+        subarray = SubarrayDescription.from_hdf(self.input_data.input_url)
         # Select a new subarray if the allowed_tels configuration is used
         self.subarray = (
-            self.input_data.subarray
+            subarray
             if self.allowed_tels is None
-            else self.input_data.subarray.select_subarray(self.allowed_tels)
+            else subarray.select_subarray(self.allowed_tels)
         )
         # Initialization of the statistics calculator
         self.stats_calculator = PixelStatisticsCalculator(
@@ -127,18 +121,11 @@ class PixelStatisticsCalculatorTool(Tool):
         # Iterate over the telescope ids and calculate the statistics
         for tel_id in self.subarray.tel_ids:
             # Read the whole dl1 images for one particular telescope
-            dl1_table = self.input_data.read_telescope_events(
+            dl1_table = self.input_data.read_telescope_events_by_id(
                 telescopes=[
                     tel_id,
                 ],
-            )
-            # Check if the dl1 table is empty and skip the telescope if so
-            if len(dl1_table) == 0:
-                self.log.warning(
-                    "No dl1 images found for telescope 'tel_id=%d'. Skipping.",
-                    tel_id,
-                )
-                continue
+            )[tel_id]
             # Check if the chunk size does not exceed the table length of the input data
             if self.stats_calculator.stats_aggregators[
                 self.stats_calculator.stats_aggregator_type.tel[tel_id]
@@ -203,6 +190,12 @@ class PixelStatisticsCalculatorTool(Tool):
         )
 
     def finish(self):
+        # Store the subarray description in the output file
+        self.subarray.to_hdf(self.output_path, overwrite=self.overwrite)
+        self.log.info(
+            "Subarray description was stored in '%s'",
+            self.output_path,
+        )
         self.log.info("Tool is shutting down")
 
 


Test patch:
diff --git a/src/ctapipe/reco/tests/test_sklearn.py b/src/ctapipe/reco/tests/test_sklearn.py
index d2d715d7..fd0e1cbe 100644
--- a/src/ctapipe/reco/tests/test_sklearn.py
+++ b/src/ctapipe/reco/tests/test_sklearn.py
@@ -6,10 +6,12 @@ from numpy.testing import assert_array_equal
 from traitlets import TraitError
 from traitlets.config import Config
 
+from ctapipe.containers import CoordinateFrameType
 from ctapipe.core import Component
 from ctapipe.reco import EnergyRegressor, ParticleClassifier
 from ctapipe.reco.reconstructor import ReconstructionProperty
 from ctapipe.reco.sklearn import DispReconstructor
+from ctapipe.utils import get_dataset_path
 
 KEY = "LST_LST_LSTCam"
 
@@ -290,3 +292,23 @@ def test_io_with_parent(example_table, tmp_path, example_subarray):
 
     with pytest.raises(TypeError):
         EnergyRegressor.read(path)
+
+
+def test_disp_fixed_icrs_pointing(disp_reconstructor_path):
+    from ctapipe.io import TableLoader
+    from ctapipe.reco.sklearn import DispReconstructor
+
+    path = get_dataset_path("gamma_diffuse_dl2_train_small.dl2.h5")
+    with TableLoader(path) as loader:
+        events = loader.read_telescope_events([1, 2, 3, 4], stop=20)
+
+    events["subarray_pointing_frame"] = CoordinateFrameType.ICRS.value
+    rng = np.random.default_rng(0)
+    events["ExtraTreesRegressor_energy"] = rng.uniform(1, 10, len(events)) * u.TeV
+    events["ExtraTreesRegressor_tel_energy"] = rng.uniform(1, 10, len(events)) * u.TeV
+
+    disp_reconstructor_path, _ = disp_reconstructor_path
+    disp_model = DispReconstructor.read(disp_reconstructor_path)
+
+    with pytest.raises(NotImplementedError, match="Only AltAz frame supported"):
+        disp_model.predict_table(disp_model.subarray.tel[1], events)
diff --git a/src/ctapipe/tools/tests/test_calculate_pixel_stats.py b/src/ctapipe/tools/tests/test_calculate_pixel_stats.py
index cd56f314..42c842ed 100644
--- a/src/ctapipe/tools/tests/test_calculate_pixel_stats.py
+++ b/src/ctapipe/tools/tests/test_calculate_pixel_stats.py
@@ -8,6 +8,7 @@ from traitlets.config.loader import Config
 
 from ctapipe.core import run_tool
 from ctapipe.core.tool import ToolConfigurationError
+from ctapipe.instrument import SubarrayDescription
 from ctapipe.io import read_table
 from ctapipe.tools.calculate_pixel_stats import PixelStatisticsCalculatorTool
 
@@ -20,13 +21,13 @@ def test_calculate_pixel_stats_tool(tmp_path, dl1_image_file):
     config = Config(
         {
             "PixelStatisticsCalculatorTool": {
-                "allowed_tels": [3],
+                "allowed_tels": [tel_id],
                 "input_column_name": "image",
                 "output_table_name": "statistics",
             },
             "PixelStatisticsCalculator": {
                 "stats_aggregator_type": [
-                    ("type", "*", "PlainAggregator"),
+                    ("id", tel_id, "PlainAggregator"),
                 ],
                 "outlier_detector_list": [
                     {
@@ -64,6 +65,10 @@ def test_calculate_pixel_stats_tool(tmp_path, dl1_image_file):
         )["mean"].ndim
         == 3
     )
+    # Read subarray description from the created monitoring file
+    subarray = SubarrayDescription.from_hdf(monitoring_file)
+    # Check for the selected telescope
+    assert subarray.tel_ids[0] == tel_id
 
 
 def test_tool_config_error(tmp_path, dl1_image_file):
@@ -73,13 +78,14 @@ def test_tool_config_error(tmp_path, dl1_image_file):
     config = Config(
         {
             "PixelStatisticsCalculatorTool": {
+                "allowed_tels": [3],
                 "input_column_name": "image_charges",
                 "output_table_name": "statistics",
             }
         }
     )
     # Set the output file path
-    monitoring_failure_colname_file = tmp_path / "monitoring_failure_colname.dl1.h5"
+    monitoring_file = tmp_path / "monitoring.dl1.h5"
     # Check if ToolConfigurationError is raised
     # when the column name of the pixel-wise image data is not correct
     with pytest.raises(
@@ -89,7 +95,7 @@ def test_tool_config_error(tmp_path, dl1_image_file):
             PixelStatisticsCalculatorTool(config=config),
             argv=[
                 f"--input_url={dl1_image_file}",
-                f"--output_path={monitoring_failure_colname_file}",
+                f"--output_path={monitoring_file}",
                 "--StatisticsAggregator.chunk_size=1",
                 "--overwrite",
             ],
@@ -97,11 +103,22 @@ def test_tool_config_error(tmp_path, dl1_image_file):
             raises=True,
         )
     # Check if ToolConfigurationError is raised
+    # when the input and output files are the same
+    with pytest.raises(
+        ToolConfigurationError, match="Input and output files are same."
+    ):
+        run_tool(
+            PixelStatisticsCalculatorTool(),
+            argv=[
+                f"--input_url={dl1_image_file}",
+                f"--output_path={dl1_image_file}",
+                "--overwrite",
+            ],
+            cwd=tmp_path,
+            raises=True,
+        )
+    # Check if ToolConfigurationError is raised
     # when the chunk size is larger than the number of events in the input file
-    monitoring_failure_chunk_size_file = (
-        tmp_path / "monitoring_failure_chunk_size.dl1.h5"
-    )
-
     with pytest.raises(
         ToolConfigurationError, match="Change --StatisticsAggregator.chunk_size"
     ):
@@ -109,8 +126,10 @@ def test_tool_config_error(tmp_path, dl1_image_file):
             PixelStatisticsCalculatorTool(),
             argv=[
                 f"--input_url={dl1_image_file}",
-                f"--output_path={monitoring_failure_chunk_size_file}",
+                f"--output_path={monitoring_file}",
+                "--PixelStatisticsCalculatorTool.allowed_tels=3",
                 "--StatisticsAggregator.chunk_size=2500",
+                "--overwrite",
             ],
             cwd=tmp_path,
             raises=True,