Merge telescope-wise data from the same OB

docs/changes/2916.feature.rst

@@ -0,0 +1 @@
+Add a new merge_strategy option 'combine-telescope-events' to HDF5Merger component to support the merging of telescope-wise data from different files for the same observation block.

src/ctapipe/instrument/subarray.py

@@ -614,12 +614,18 @@ def to_hdf(self, h5file, overwrite=False, mode="a"):
         """
         # here to prevent circular import
         from ..io import write_table
+        from ..io.hdf5dataformat import (
+            CONFIG_INSTRUMENT_SUBARRAY,
+            CONFIG_INSTRUMENT_SUBARRAY_LAYOUT,
+            CONFIG_INSTRUMENT_TEL_CAMERA,
+            CONFIG_INSTRUMENT_TEL_OPTICS,
+        )
 
         with ExitStack() as stack:
             if not isinstance(h5file, tables.File):
                 h5file = stack.enter_context(tables.open_file(h5file, mode=mode))
 
-            if "/configuration/instrument/subarray" in h5file.root and not overwrite:
+            if CONFIG_INSTRUMENT_SUBARRAY in h5file.root and not overwrite:
                 raise OSError(
                     "File already contains a SubarrayDescription and overwrite=False"
                 )
@@ -630,40 +636,43 @@ def to_hdf(self, h5file, overwrite=False, mode="a"):
             write_table(
                 subarray_table,
                 h5file,
-                path="/configuration/instrument/subarray/layout",
+                path=CONFIG_INSTRUMENT_SUBARRAY_LAYOUT,
                 overwrite=overwrite,
             )
             write_table(
                 self.to_table(kind="optics"),
                 h5file,
-                path="/configuration/instrument/telescope/optics",
+                path=CONFIG_INSTRUMENT_TEL_OPTICS,
                 overwrite=overwrite,
             )
             for i, camera in enumerate(self.camera_types):
                 write_table(
                     camera.geometry.to_table(),
                     h5file,
-                    path=f"/configuration/instrument/telescope/camera/geometry_{i}",
+                    path=f"{CONFIG_INSTRUMENT_TEL_CAMERA}/geometry_{i}",
                     overwrite=overwrite,
                 )
                 write_table(
                     camera.readout.to_table(),
                     h5file,
-                    path=f"/configuration/instrument/telescope/camera/readout_{i}",
+                    path=f"{CONFIG_INSTRUMENT_TEL_CAMERA}/readout_{i}",
                     overwrite=overwrite,
                 )
 
     @classmethod
     def from_hdf(cls, path, focal_length_choice=FocalLengthKind.EFFECTIVE):
         # here to prevent circular import
         from ..io import read_table
+        from ..io.hdf5dataformat import (
+            CONFIG_INSTRUMENT_SUBARRAY_LAYOUT,
+            CONFIG_INSTRUMENT_TEL_CAMERA,
+            CONFIG_INSTRUMENT_TEL_OPTICS,
+        )
 
         if isinstance(focal_length_choice, str):
             focal_length_choice = FocalLengthKind[focal_length_choice.upper()]
 
-        layout = read_table(
-            path, "/configuration/instrument/subarray/layout", table_cls=QTable
-        )
+        layout = read_table(path, CONFIG_INSTRUMENT_SUBARRAY_LAYOUT, table_cls=QTable)
 
         version = layout.meta.get("TAB_VER")
         if version not in cls.COMPATIBLE_VERSIONS:
@@ -677,22 +686,16 @@ def from_hdf(cls, path, focal_length_choice=FocalLengthKind.EFFECTIVE):
 
         for idx in set(layout["camera_index"]):
             geometry = CameraGeometry.from_table(
-                read_table(
-                    path, f"/configuration/instrument/telescope/camera/geometry_{idx}"
-                )
+                read_table(path, f"{CONFIG_INSTRUMENT_TEL_CAMERA}/geometry_{idx}")
             )
             readout = CameraReadout.from_table(
-                read_table(
-                    path, f"/configuration/instrument/telescope/camera/readout_{idx}"
-                )
+                read_table(path, f"{CONFIG_INSTRUMENT_TEL_CAMERA}/readout_{idx}")
             )
             cameras[idx] = CameraDescription(
                 name=geometry.name, readout=readout, geometry=geometry
             )
 
-        optics_table = read_table(
-            path, "/configuration/instrument/telescope/optics", table_cls=QTable
-        )
+        optics_table = read_table(path, CONFIG_INSTRUMENT_TEL_OPTICS, table_cls=QTable)
 
         optics_version = optics_table.meta.get("TAB_VER")
         if optics_version not in OpticsDescription.COMPATIBLE_VERSIONS:
@@ -808,3 +811,81 @@ def check_matching_subarrays(subarray_list: list) -> bool:
             set(subarray.tel_ids) == set(subarray_list[0].tel_ids)
             for subarray in subarray_list
         )
+
+    @staticmethod
+    def merge_subarrays(
+        subarray_list: list,
+        name=None,
+        overwrite_tel_ids: bool = False,
+    ) -> "SubarrayDescription":
+        """Merge multiple subarrays into one
+
+        Parameters
+        ----------
+        subarray_list: list(SubarrayDescription)
+            list of subarrays to merge
+        name: str
+            name of new merged subarray
+        overwrite_tel_ids: bool
+            if True, telescope entries from later subarrays replace earlier ones;
+            if False (default), encountering a duplicate telescope id raises a
+            ValueError.
+        Returns
+        -------
+        SubarrayDescription
+
+        """
+
+        tel_positions, tel_descriptions = {}, {}
+        tel_ids, tid_to_subarray = set(), {}
+        reference_location = subarray_list[
+            0
+        ].reference_location  # Get the reference location from the first subarray
+        for s, subarray in enumerate(subarray_list):
+            if not isinstance(subarray, SubarrayDescription):
+                raise TypeError(
+                    "All elements of subarray_list must be 'SubarrayDescription' "
+                    f"instances, got '{type(subarray)}' for element '{s}'."
+                )
+            for tid in subarray.tel_ids:
+                if tid in tel_ids:
+                    if overwrite_tel_ids:
+                        # Warn about overwriting telescope entry
+                        msg = (
+                            f"Overwriting telescope id '{tid}' from subarray "
+                            f"'{subarray.name}' into merged subarray."
+                        )
+                        warnings.warn(msg, UserWarning)
+                    else:
+                        raise ValueError(
+                            "Duplicate telescope id encountered while merging subarrays. "
+                            f"Telescope '{tid}' already defined; set overwrite_tel_ids=True to "
+                            "allow later subarrays to replace earlier entries."
+                        )
+                tid_to_subarray[tid] = subarray
+                tel_ids.add(tid)
+            if subarray.reference_location != reference_location:
+                raise ValueError(
+                    "All subarrays must have the same reference_location to be merged. "
+                    f"Subarray '{subarray.name}' ({subarray.reference_location}) does not match "
+                    f"the reference location of the first subarray '{subarray_list[0].name}' "
+                    f"({reference_location})."
+                )
+
+        # Merge telescope positions and descriptions, optionally allowing later entries to overwrite
+        for tid in sorted(tel_ids):
+            # Copy/overwrite telescope position and description from the current subarray
+            subarray = tid_to_subarray[tid]
+            tel_positions[tid] = subarray.positions[tid]
+            tel_descriptions[tid] = subarray.tels[tid]
+
+        if name is None:
+            name = "Merged_" + _range_extraction(sorted(tel_ids))
+
+        newsub = SubarrayDescription(
+            name,
+            tel_positions=tel_positions,
+            tel_descriptions=tel_descriptions,
+            reference_location=reference_location,
+        )
+        return newsub

src/ctapipe/instrument/tests/test_subarray.py

@@ -313,3 +313,71 @@ def test_check_matchings_subarray(example_subarray, subarray_prod5_paranal):
     assert not SubarrayDescription.check_matching_subarrays(
         [example_subarray, subarray_prod5_paranal]
     )
+
+
+def test_merge_subarrays(example_subarray):
+    """Test SubarrayDescription.merge_subarrays static method"""
+
+    sub1 = example_subarray.select_subarray([1, 2], name="s1")
+    sub2 = example_subarray.select_subarray([3, 4], name="s2")
+    expected_sub = example_subarray.select_subarray([1, 2, 3, 4], name="Merged_1-4")
+    merged_sub = SubarrayDescription.merge_subarrays([sub1, sub2])
+
+    assert expected_sub == merged_sub
+
+
+def test_merge_subarrays_exceptions(example_subarray):
+    """Merging subarrays with invalid parameters should raise exceptions."""
+
+    sub1 = example_subarray.select_subarray([1, 2], name="s1")
+    sub2 = example_subarray.select_subarray([3, 4], name="s2")
+
+    # Check that invalid inputs raise exceptions
+    with pytest.raises(
+        TypeError, match="All elements of subarray_list must be 'SubarrayDescription'"
+    ):
+        SubarrayDescription.merge_subarrays([sub1, int(67)])
+
+    # Check that duplicate telescope ids without overwrite raises exception
+    with pytest.raises(ValueError, match="Duplicate telescope id encountered"):
+        SubarrayDescription.merge_subarrays([sub1, sub1], overwrite_tel_ids=False)
+
+    # Check that different reference locations raises exception
+    shifted_location = EarthLocation(
+        lon=sub2.reference_location.lon,
+        lat=sub2.reference_location.lat,
+        height=sub2.reference_location.height + 1 * u.m,
+    )
+    sub2_shifted = SubarrayDescription(
+        name="shifted",
+        tel_positions=sub2.positions,
+        tel_descriptions=sub2.tel,
+        reference_location=shifted_location,
+    )
+    with pytest.raises(
+        ValueError, match="All subarrays must have the same reference_location"
+    ):
+        SubarrayDescription.merge_subarrays(
+            [sub1, sub2_shifted], overwrite_tel_ids=True
+        )
+
+
+def test_merge_subarrays_overwrite_tel_ids(example_subarray):
+    """Later subarrays overwrite earlier telescope entries when enabled."""
+
+    sub1 = example_subarray.select_subarray([1, 2], name="s1")
+    sub2 = example_subarray.select_subarray([1, 2], name="s2")
+
+    # Overwrite positions of the telescopes in sub2
+    sub2.positions[1] = np.array([10, 0, 0]) * u.m
+    sub2.positions[2] = np.array([-10, 0, 0]) * u.m
+
+    with pytest.warns(UserWarning, match="Overwriting telescope id"):
+        merged_sub = SubarrayDescription.merge_subarrays(
+            [sub1, sub2], overwrite_tel_ids=True
+        )
+
+    assert merged_sub.name == "Merged_1,2"
+    assert merged_sub.n_tels == 2
+    np.testing.assert_allclose(merged_sub.positions[1], [10, 0, 0] * u.m)
+    np.testing.assert_allclose(merged_sub.positions[2], [-10, 0, 0] * u.m)

src/ctapipe/io/hdf5dataformat.py

@@ -9,6 +9,11 @@
     "SIMULATION_GROUP",
     "SIMULATION_TEL_TABLE",
     "CONFIG_GROUP",
+    "CONFIG_INSTRUMENT_SUBARRAY",
+    "CONFIG_INSTRUMENT_SUBARRAY_LAYOUT",
+    "CONFIG_INSTRUMENT_TEL",
+    "CONFIG_INSTRUMENT_TEL_OPTICS",
+    "CONFIG_INSTRUMENT_TEL_CAMERA",
     "SCHEDULING_BLOCK_TABLE",
     "OBSERVATION_BLOCK_TABLE",
     "SIMULATION_RUN_TABLE",
@@ -38,6 +43,7 @@
     "DL2_TEL_GEOMETRY_GROUP",
     "DL2_TEL_ENERGY_GROUP",
     "DL2_TEL_PARTICLETYPE_GROUP",
+    "DL2_TEL_IMPACT_GROUP",
     "DL0_TEL_POINTING_GROUP",
     "DL1_SUBARRAY_POINTING_GROUP",
     "DL1_TEL_POINTING_GROUP",
@@ -60,10 +66,16 @@
 
 # Configuration, service, and simulation group
 CONFIG_GROUP = "/configuration"
+CONFIG_INSTRUMENT_SUBARRAY = "/configuration/instrument/subarray"
+CONFIG_INSTRUMENT_SUBARRAY_LAYOUT = "/configuration/instrument/subarray/layout"
+CONFIG_INSTRUMENT_TEL = "/configuration/instrument/telescope"
+CONFIG_INSTRUMENT_TEL_OPTICS = "/configuration/instrument/telescope/optics"
+CONFIG_INSTRUMENT_TEL_CAMERA = "/configuration/instrument/telescope/camera"
 SCHEDULING_BLOCK_TABLE = "/configuration/observation/scheduling_block"
 OBSERVATION_BLOCK_TABLE = "/configuration/observation/observation_block"
 SIMULATION_RUN_TABLE = "/configuration/simulation/run"
 FIXED_POINTING_GROUP = "/configuration/telescope/pointing"
+
 DL1_IMAGE_STATISTICS_TABLE = "/dl1/service/image_statistics"
 DL2_EVENT_STATISTICS_GROUP = "/dl2/service/tel_event_statistics"
 SIMULATION_GROUP = "/simulation"
@@ -93,6 +105,7 @@
 DL2_TEL_GEOMETRY_GROUP = "/dl2/event/telescope/geometry"
 DL2_TEL_ENERGY_GROUP = "/dl2/event/telescope/energy"
 DL2_TEL_PARTICLETYPE_GROUP = "/dl2/event/telescope/particle_type"
+DL2_TEL_IMPACT_GROUP = "/dl2/event/telescope/impact"
 
 DL2_GROUP = "/dl2"
 DL2_SUBARRAY_GROUP = "/dl2/event/subarray"