Repository: cta-observatory/ctapipe
Pull request #1980 (2022-06-30 16:11:58)

Add location of subarray reference point to `SubarrayDescription`


Patch:
diff --git a/ctapipe/calib/camera/flatfield.py b/ctapipe/calib/camera/flatfield.py
index 1e3263c6..18411057 100644
--- a/ctapipe/calib/camera/flatfield.py
+++ b/ctapipe/calib/camera/flatfield.py
@@ -115,18 +115,18 @@ class FlatFieldCalculator(Component):
 
 class FlasherFlatFieldCalculator(FlatFieldCalculator):
     """Calculates flat-field parameters from flasher data
-      based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
-      Pixels are defined as outliers on the base of a cut on the pixel charge median
-      over the full sample distribution and the pixel signal time inside the
-      waveform time
+       based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
+       Pixels are defined as outliers on the base of a cut on the pixel charge median
+       over the full sample distribution and the pixel signal time inside the
+       waveform time
 
 
-    Parameters
-    ----------
-    charge_cut_outliers : List[2]
-        Interval of accepted charge values (fraction with respect to camera median value)
-    time_cut_outliers : List[2]
-        Interval (in waveform samples) of accepted time values
+     Parameters
+     ----------
+     charge_cut_outliers : List[2]
+         Interval of accepted charge values (fraction with respect to camera median value)
+     time_cut_outliers : List[2]
+         Interval (in waveform samples) of accepted time values
 
     """
 
@@ -140,18 +140,18 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
     def __init__(self, **kwargs):
         """Calculates flat-field parameters from flasher data
-          based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
-          Pixels are defined as outliers on the base of a cut on the pixel charge median
-          over the full sample distribution and the pixel signal time inside the
-          waveform time
+           based on the best algorithm described by S. Fegan in MST-CAM-TN-0060 (eq. 19)
+           Pixels are defined as outliers on the base of a cut on the pixel charge median
+           over the full sample distribution and the pixel signal time inside the
+           waveform time
 
 
-        Parameters:
-        ----------
-        charge_cut_outliers : List[2]
-            Interval of accepted charge values (fraction with respect to camera median value)
-        time_cut_outliers : List[2]
-            Interval (in waveform samples) of accepted time values
+         Parameters:
+         ----------
+         charge_cut_outliers : List[2]
+             Interval of accepted charge values (fraction with respect to camera median value)
+         time_cut_outliers : List[2]
+             Interval (in waveform samples) of accepted time values
 
         """
         super().__init__(**kwargs)
@@ -190,14 +190,14 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
     def calculate_relative_gain(self, event):
         """
-        calculate the flatfield statistical values
-        and fill mon.tel[tel_id].flatfield container
+         calculate the flatfield statistical values
+         and fill mon.tel[tel_id].flatfield container
 
-        Parameters
-        ----------
-        event : general event container
+         Parameters
+         ----------
+         event : general event container
 
-        """
+         """
 
         # initialize the np array at each cycle
         waveform = event.r1.tel[self.tel_id].waveform
@@ -235,7 +235,7 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
 
         self.collect_sample(dl1.image, pixel_mask, dl1.peak_time)
 
-        sample_age = (trigger_time - self.time_start).to_value(u.s)
+        sample_age = trigger_time - self.time_start
 
         # check if to create a calibration event
         if (
@@ -296,7 +296,7 @@ class FlasherFlatFieldCalculator(FlatFieldCalculator):
     def calculate_time_results(
         self, trace_time, masked_pixels_of_sample, time_start, trigger_time
     ):
-        """Calculate and return the time results"""
+        """Calculate and return the time results """
         masked_trace_time = np.ma.array(trace_time, mask=masked_pixels_of_sample)
 
         # median over the sample per pixel
diff --git a/ctapipe/calib/camera/pedestals.py b/ctapipe/calib/camera/pedestals.py
index ed8c527d..3bdcaceb 100644
--- a/ctapipe/calib/camera/pedestals.py
+++ b/ctapipe/calib/camera/pedestals.py
@@ -72,7 +72,8 @@ class PedestalCalculator(Component):
         Set to None if no configuration to pass.
 
     kwargs
-    """
+
+"""
 
     tel_id = Int(0, help="id of the telescope to calculate the pedestal values").tag(
         config=True
@@ -114,7 +115,7 @@ class PedestalCalculator(Component):
 
         kwargs
 
-        """
+    """
 
         super().__init__(**kwargs)
 
@@ -142,21 +143,21 @@ class PedestalCalculator(Component):
 
 class PedestalIntegrator(PedestalCalculator):
     """Calculates pedestal parameters integrating the charge of pedestal events:
-      the pedestal value corresponds to the charge estimated with the selected
-      charge extractor
-      The pixels are set as outliers on the base of a cut on the pixel charge median
-      over the pedestal sample and the pixel charge standard deviation over
-      the pedestal sample with respect to the camera median values
+       the pedestal value corresponds to the charge estimated with the selected
+       charge extractor
+       The pixels are set as outliers on the base of a cut on the pixel charge median
+       over the pedestal sample and the pixel charge standard deviation over
+       the pedestal sample with respect to the camera median values
 
 
-    Parameters:
-    ----------
-    charge_median_cut_outliers : List[2]
-        Interval (number of std) of accepted charge values around camera median value
-    charge_std_cut_outliers : List[2]
-        Interval (number of std) of accepted charge standard deviation around camera median value
+     Parameters:
+     ----------
+     charge_median_cut_outliers : List[2]
+         Interval (number of std) of accepted charge values around camera median value
+     charge_std_cut_outliers : List[2]
+         Interval (number of std) of accepted charge standard deviation around camera median value
 
-    """
+     """
 
     charge_median_cut_outliers = List(
         [-3, 3],
@@ -169,19 +170,19 @@ class PedestalIntegrator(PedestalCalculator):
 
     def __init__(self, **kwargs):
         """Calculates pedestal parameters integrating the charge of pedestal events:
-          the pedestal value corresponds to the charge estimated with the selected
-          charge extractor
-          The pixels are set as outliers on the base of a cut on the pixel charge median
-          over the pedestal sample and the pixel charge standard deviation over
-          the pedestal sample with respect to the camera median values
-
-
-        Parameters:
-        ----------
-        charge_median_cut_outliers : List[2]
-            Interval (number of std) of accepted charge values around camera median value
-        charge_std_cut_outliers : List[2]
-            Interval (number of std) of accepted charge standard deviation around camera median value
+           the pedestal value corresponds to the charge estimated with the selected
+           charge extractor
+           The pixels are set as outliers on the base of a cut on the pixel charge median
+           over the pedestal sample and the pixel charge standard deviation over
+           the pedestal sample with respect to the camera median values
+
+
+         Parameters:
+         ----------
+         charge_median_cut_outliers : List[2]
+             Interval (number of std) of accepted charge values around camera median value
+         charge_std_cut_outliers : List[2]
+             Interval (number of std) of accepted charge standard deviation around camera median value
         """
 
         super().__init__(**kwargs)
@@ -256,7 +257,7 @@ class PedestalIntegrator(PedestalCalculator):
 
         self.collect_sample(dl1.image, pixel_mask)
 
-        sample_age = (trigger_time - self.time_start).to_value(u.s)
+        sample_age = trigger_time - self.time_start
 
         # check if to create a calibration event
         if (
diff --git a/ctapipe/core/container.py b/ctapipe/core/container.py
index 68322fc7..cb1b536c 100644
--- a/ctapipe/core/container.py
+++ b/ctapipe/core/container.py
@@ -1,11 +1,10 @@
 from collections import defaultdict
 from copy import deepcopy
 from pprint import pformat
-from textwrap import wrap, dedent
+from textwrap import wrap
 import warnings
 import numpy as np
 from astropy.units import UnitConversionError, Quantity, Unit
-from inspect import isclass
 
 import logging
 
@@ -25,25 +24,25 @@ class Field:
 
     Parameters
     ----------
-    default :
-        default value of the item. This will be set when the `Container`
+    default:
+        default value of the item (this will be set when the `Container`
         is constructed, as well as when  ``Container.reset`` is called
-    description : str
+    description: str
         Help text associated with the item
-    unit : str or astropy.units.core.UnitBase
+    unit: str or astropy.units.core.UnitBase
         unit to convert to when writing output, or None for no conversion
-    ucd : str
+    ucd: str
         universal content descriptor (see Virtual Observatory standards)
-    type : type
+    type: type
         expected type of value
-    dtype : str or np.dtype
+    dtype: str or np.dtype
         expected data type of the value, None to ignore in validation.
         Means value is expected to be a numpy array or astropy quantity
-    ndim : int or None
+    ndim: int or None
         expected dimensionality of the data, for arrays, None to ignore
-    allow_none : bool
+    allow_none:
         if the value of None is given to this Field, skip validation
-    max_len : int
+    max_len:
         if type is str, max_len is the maximum number of bytes of the utf-8
         encoded string to be used.
     """
@@ -60,6 +59,7 @@ class Field:
         allow_none=True,
         max_length=None,
     ):
+
         self.default = default
         self.description = description
         self.unit = Unit(unit) if unit is not None else None
@@ -71,34 +71,7 @@ class Field:
         self.max_length = max_length
 
     def __repr__(self):
-        if isinstance(self.default, Container):
-            default = f"{self.default.__class__.__name__}"
-        elif isinstance(self.default, Map):
-            if isclass(self.default.default_factory):
-                cls = self.default.default_factory
-                default = f"Map({cls.__module__}.{cls.__name__})"
-            else:
-                default = f"Map({repr(self.default.default_factory)}"
-        else:
-            default = str(self.default)
-        cmps = [f"Field(default={default}"]
-        if self.unit is not None:
-            cmps.append(f", unit={self.unit}")
-        if self.dtype is not None:
-            cmps.append(f", dtype={self.dtype}")
-        if self.ndim is not None:
-            cmps.append(f", ndim={self.ndim}")
-        if self.type is not None:
-            cmps.append(f", type={self.type.__name__}")
-        if self.allow_none is False:
-            cmps.append(", allow_none=False")
-        if self.max_length is not None:
-            cmps.append(f", max_length={self.max_length}")
-        cmps.append(")")
-        return "".join(cmps)
-
-    def __str__(self):
-        desc = f"{self.description} with default {self.default}"
+        desc = f"{self.description}"
         if self.unit is not None:
             desc += f" [{self.unit}]"
         if self.ndim is not None:
@@ -180,23 +153,6 @@ class DeprecatedField(Field):
         self.reason = reason
 
 
-_doc_template = """{doc}
-
-Attributes
-----------
-{fields}
-meta : dict
-    dict of attached metadata
-prefix : str
-    Prefix attached to column names when saved to a table or file
-"""
-
-
-def _build_docstring(doc, fields):
-    fields = [f"{k} : {f!r}\n    {f.description}" for k, f in fields.items()]
-    return _doc_template.format(doc=dedent(doc), fields="\n".join(fields))
-
-
 class ContainerMeta(type):
     """
     The MetaClass for `Container`
@@ -223,8 +179,6 @@ class ContainerMeta(type):
         for k in field_names:
             dct["fields"][k] = dct.pop(k)
 
-        dct["__doc__"] = _build_docstring(dct.get("__doc__", ""), dct["fields"])
-
         new_cls = type.__new__(cls, name, bases, dct)
 
         # if prefix was not set as a class variable, build a default one
@@ -398,7 +352,7 @@ class Container(metaclass=ContainerMeta):
                 extra = ".*"
             if isinstance(getattr(self, name), Map):
                 extra = "[*]"
-            desc = "{:>30s}: {}".format(name + extra, str(item))
+            desc = "{:>30s}: {}".format(name + extra, repr(item))
             lines = wrap(desc, 80, subsequent_indent=" " * 32)
             text.extend(lines)
         return "\n".join(text)
diff --git a/ctapipe/image/hillas.py b/ctapipe/image/hillas.py
index fae58ccb..8fdd405a 100644
--- a/ctapipe/image/hillas.py
+++ b/ctapipe/image/hillas.py
@@ -7,6 +7,7 @@ Hillas-style moment-based shower image parametrization.
 import astropy.units as u
 import numpy as np
 from astropy.coordinates import Angle
+from astropy.units import Quantity
 from ..containers import CameraHillasParametersContainer, HillasParametersContainer
 
 
@@ -104,15 +105,12 @@ def hillas_parameters(geom, image):
         container of hillas parametesr
     """
     unit = geom.pix_x.unit
-    pix_x = geom.pix_x.to_value(unit)
-    pix_y = geom.pix_y.to_value(unit)
+    pix_x = Quantity(np.asanyarray(geom.pix_x, dtype=np.float64)).value
+    pix_y = Quantity(np.asanyarray(geom.pix_y, dtype=np.float64)).value
     image = np.asanyarray(image, dtype=np.float64)
-
-    if isinstance(image, np.ma.masked_array):
-        image = np.ma.filled(image, 0)
-
-    if not (pix_x.shape == pix_y.shape == image.shape):
-        raise ValueError("Image and pixel shape do not match")
+    image = np.ma.filled(image, 0)
+    msg = "Image and pixel shape do not match"
+    assert pix_x.shape == pix_y.shape == image.shape, msg
 
     size = np.sum(image)
 
@@ -160,11 +158,11 @@ def hillas_parameters(geom, image):
         # calculate higher order moments along shower axes
         longitudinal = delta_x * np.cos(psi) + delta_y * np.sin(psi)
 
-        m3_long = np.average(longitudinal**3, weights=image)
-        skewness_long = m3_long / length**3
+        m3_long = np.average(longitudinal ** 3, weights=image)
+        skewness_long = m3_long / length ** 3
 
-        m4_long = np.average(longitudinal**4, weights=image)
-        kurtosis_long = m4_long / length**4
+        m4_long = np.average(longitudinal ** 4, weights=image)
+        kurtosis_long = m4_long / length ** 4
 
     # Compute of the Hillas parameters uncertainties.
     # Implementation described in [hillas_uncertainties]_ This is an internal MAGIC document
@@ -176,8 +174,8 @@ def hillas_parameters(geom, image):
     b = (1 - cos_2psi) / 2
     c = np.sin(2 * psi)
 
-    A = ((delta_x**2.0) - cov[0][0]) / size
-    B = ((delta_y**2.0) - cov[1][1]) / size
+    A = ((delta_x ** 2.0) - cov[0][0]) / size
+    B = ((delta_y ** 2.0) - cov[1][1]) / size
     C = ((delta_x * delta_y) - cov[0][1]) / size
 
     # Hillas's uncertainties
diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py
index 853ddfb6..5e4f8720 100644
--- a/ctapipe/instrument/subarray.py
+++ b/ctapipe/instrument/subarray.py
@@ -6,6 +6,7 @@ from contextlib import ExitStack
 from copy import copy
 from itertools import groupby
 from typing import Dict, List, Union
+from astropy.coordinates.earth import EarthLocation
 
 import numpy as np
 import tables
@@ -51,16 +52,6 @@ class SubarrayDescription:
     Collects the `~ctapipe.instrument.TelescopeDescription` of all telescopes
     along with their positions on the ground.
 
-    Parameters
-    ----------
-    name: str
-        name of this subarray
-    tel_positions: Dict[Array]
-        dict of x,y,z telescope positions on the ground by tel_id. These are
-        converted internally to a coordinate in the `~ctapipe.coordinates.GroundFrame`
-    tel_descriptions: Dict[TelescopeDescription]
-        dict of TelescopeDescriptions by tel_id
-
     Attributes
     ----------
     name: str
@@ -69,16 +60,35 @@ class SubarrayDescription:
        coordinates of all telescopes
     tels:
        dict of TelescopeDescription for each telescope in the subarray
-    tel_ids: np.ndarray
-        array of tel_ids
-    tel_indices: dict
-        dict mapping tel_id to index in array attributes
     """
 
-    def __init__(self, name, tel_positions=None, tel_descriptions=None):
+    def __init__(
+        self,
+        name,
+        tel_positions=None,
+        tel_descriptions=None,
+        reference_location=None,
+    ):
+        """
+        Initialize a new SubarrayDescription
+
+        Parameters
+        ----------
+        name : str
+            name of this subarray
+        tel_positions : Dict[int, np.ndarray]
+            dict of x,y,z telescope positions on the ground by tel_id. These are
+            converted internally to a coordinate in the `~ctapipe.coordinates.GroundFrame`
+        tel_descriptions : Dict[TelescopeDescription]
+            dict of TelescopeDescriptions by tel_id
+        reference_location : `astropy.coordinates.EarthLocation`
+            EarthLocation of the array reference position, (0, 0, 0) of the
+            coordinate system used for `tel_positions`.
+        """
         self.name = name
         self.positions = tel_positions or dict()
         self.tels: Dict[int, TelescopeDescription] = tel_descriptions or dict()
+        self.reference_location = reference_location
 
         if self.positions.keys() != self.tels.keys():
             raise ValueError("Telescope ids in positions and descriptions do not match")
@@ -92,8 +102,8 @@ class SubarrayDescription:
         )
 
     @property
-    def tel(self):
-        """for backward compatibility"""
+    def tel(self) -> Dict[int, TelescopeDescription]:
+        """Dictionary mapping tel_ids to TelescopeDescriptions"""
         return self.tels
 
     @property
@@ -238,7 +248,7 @@ class SubarrayDescription:
             which table to generate (subarray or optics)
         """
 
-        if kind == 'joined':
+        if kind == "joined":
             table = self.to_table()
             optics = self.to_table(kind="optics")
             optics["optics_index"] = np.arange(len(optics))
@@ -252,10 +262,9 @@ class SubarrayDescription:
             )
 
             table.remove_columns(["optics_index", "camera_index"])
-            table.add_index('tel_id')
+            table.add_index("tel_id")
             return table
 
-
         meta = {
             "ORIGIN": "ctapipe.instrument.SubarrayDescription",
             "SUBARRAY": self.name,
@@ -263,6 +272,12 @@ class SubarrayDescription:
             "TAB_TYPE": kind,
         }
 
+        if self.reference_location is not None:
+            itrs = self.reference_location.itrs
+            meta["OBSGEO-X"] = itrs.x.to_value(u.m)
+            meta["OBSGEO-Y"] = itrs.y.to_value(u.m)
+            meta["OBSGEO-Z"] = itrs.z.to_value(u.m)
+
         if kind == "subarray":
 
             unique_types = self.telescope_types
@@ -447,6 +462,9 @@ class SubarrayDescription:
         if self.positions.keys() != other.positions.keys():
             return False
 
+        if self.reference_location != other.reference_location:
+            return False
+
         for tel_id in self.tels.keys():
             if self.tels[tel_id] != other.tels[tel_id]:
                 return False
@@ -454,6 +472,7 @@ class SubarrayDescription:
         for tel_id in self.tels.keys():
             if np.any(self.positions[tel_id] != other.positions[tel_id]):
                 return False
+
         return True
 
     def to_hdf(self, h5file, overwrite=False, mode="a"):
@@ -510,7 +529,13 @@ class SubarrayDescription:
                     overwrite=overwrite,
                 )
 
-            h5file.root.configuration.instrument.subarray._v_attrs.name = self.name
+            meta = h5file.root.configuration.instrument.subarray._v_attrs
+            meta["name"] = self.name
+            if self.reference_location is not None:
+                itrs = self.reference_location.itrs
+                meta["reference_itrs_x"] = itrs.x.to_value(u.m)
+                meta["reference_itrs_y"] = itrs.y.to_value(u.m)
+                meta["reference_itrs_z"] = itrs.z.to_value(u.m)
 
     @classmethod
     def from_hdf(cls, path):
@@ -588,6 +613,7 @@ class SubarrayDescription:
         positions = np.column_stack([layout[f"pos_{c}"] for c in "xyz"])
 
         name = "Unknown"
+        reference_location = None
         with ExitStack() as stack:
             if not isinstance(path, tables.File):
                 path = stack.enter_context(tables.open_file(path, mode="r"))
@@ -596,10 +622,18 @@ class SubarrayDescription:
             if "name" in attrs:
                 name = str(attrs.name)
 
+            if "reference_itrs_x" in attrs:
+                reference_location = EarthLocation(
+                    x=attrs["reference_itrs_x"] * u.m,
+                    y=attrs["reference_itrs_y"] * u.m,
+                    z=attrs["reference_itrs_z"] * u.m,
+                )
+
         return cls(
             name=name,
             tel_positions={
                 tel_id: pos for tel_id, pos in zip(layout["tel_id"], positions)
             },
             tel_descriptions=telescope_descriptions,
+            reference_location=reference_location,
         )
diff --git a/ctapipe/io/simteleventsource.py b/ctapipe/io/simteleventsource.py
index d845feec..c6e0d3e7 100644
--- a/ctapipe/io/simteleventsource.py
+++ b/ctapipe/io/simteleventsource.py
@@ -1,3 +1,4 @@
+from collections import defaultdict
 import warnings
 from gzip import GzipFile
 from io import BufferedReader
@@ -28,7 +29,6 @@ from ..containers import (
     TriggerContainer,
 )
 from ..coordinates import CameraFrame
-from ..core import Map
 from ..core.traits import (
     Bool,
     CaselessStrEnum,
@@ -567,7 +567,7 @@ class SimTelEventSource(EventSource):
 
         central_time = parse_simtel_time(trigger["gps_time"])
 
-        tel = Map(TelescopeTriggerContainer)
+        tel = defaultdict(TelescopeTriggerContainer)
         for tel_id, time in zip(
             trigger["triggered_telescopes"], trigger["trigger_times"]
         ):
diff --git a/ctapipe/io/tableloader.py b/ctapipe/io/tableloader.py
index 39fe12c6..ba83e5ae 100644
--- a/ctapipe/io/tableloader.py
+++ b/ctapipe/io/tableloader.py
@@ -97,6 +97,14 @@ def _empty_telescope_events_table():
     return Table(names=TELESCOPE_EVENT_KEYS, dtype=[np.int32, np.int64, np.int16])
 
 
+def _empty_subarray_events_table():
+    """
+    Create a new astropy table with correct column names and dtypes
+    for subarray event based data.
+    """
+    return Table(names=SUBARRAY_EVENT_KEYS, dtype=[np.int32, np.int64])
+
+
 def _get_structure(h5file):
     """
     Check if data is stored by telescope id or by telescope type in h5file.
@@ -127,18 +135,12 @@ def _add_column_prefix(table, prefix, ignore=()):
 
 def _join_subarray_events(table1, table2):
     """Outer join two tables on the telescope subarray keys"""
-    return join_allow_empty(table1, table2, SUBARRAY_EVENT_KEYS, "left")
+    return join_allow_empty(table1, table2, SUBARRAY_EVENT_KEYS, "outer")
 
 
 def _join_telescope_events(table1, table2):
     """Outer join two tables on the telescope event keys"""
-    # we start with an empty table, but after the first non-empty, we perform
-    # left joins
-    if len(table1) == 0:
-        how = "right"
-    else:
-        how = "left"
-    return join_allow_empty(table1, table2, TELESCOPE_EVENT_KEYS, how)
+    return join_allow_empty(table1, table2, TELESCOPE_EVENT_KEYS, "outer")
 
 
 class TableLoader(Component):
@@ -186,7 +188,9 @@ class TableLoader(Component):
     load_true_parameters = traits.Bool(
         False, help="load image parameters obtained from true images"
     ).tag(config=True)
-
+    load_trigger = traits.Bool(True, help="load subarray trigger information").tag(
+        config=True
+    )
     load_instrument = traits.Bool(
         False, help="join subarray instrument information to each event"
     ).tag(config=True)
@@ -222,7 +226,7 @@ class TableLoader(Component):
 
         self.instrument_table = None
         if self.load_instrument:
-            self.instrument_table = self.subarray.to_table("joined")
+            self.instrument_table = self.subarray.to_table('joined')
 
     def close(self):
         """Close the underlying hdf5 file"""
@@ -238,10 +242,6 @@ class TableLoader(Component):
     def __exit__(self, exc_type, exc_value, traceback):
         self.close()
 
-    def __len__(self):
-        """Number of subarray events in input file"""
-        return self.h5file.root[TRIGGER_TABLE].shape[0]
-
     def _read_telescope_table(self, group, tel_id, start=None, stop=None):
         if self.structure == "by_id":
             key = f"{group}/tel_{tel_id:03d}"
@@ -259,37 +259,7 @@ class TableLoader(Component):
 
         return table
 
-    def _get_sort_index(self, start=None, stop=None):
-        """
-        Get plain index of increasing integers in the order in the file.
-
-        Astropy.table.join orders by the join keys, we want to keep
-        the original order. Adding an index before the first join and then
-        using that in the end to sort back to the original order solves this.
-        """
-        table = read_table(
-            self.h5file,
-            TRIGGER_TABLE,
-            start=start,
-            stop=stop,
-        )[["obs_id", "event_id"]]
-        self._add_index_if_needed(table)
-        return table
-
-    @staticmethod
-    def _sort_to_original_order(table, include_tel_id=False):
-        if include_tel_id:
-            table.sort(("__index__", "tel_id"))
-        else:
-            table.sort("__index__")
-        table.remove_column("__index__")
-
-    @staticmethod
-    def _add_index_if_needed(table):
-        if "__index__" not in table.colnames:
-            table["__index__"] = np.arange(len(table))
-
-    def read_subarray_events(self, start=None, stop=None, keep_order=True):
+    def read_subarray_events(self, start=None, stop=None):
         """Read subarray-based event information.
 
         Returns
@@ -297,9 +267,11 @@ class TableLoader(Component):
         table: astropy.io.Table
             Table with primary index columns "obs_id" and "event_id".
         """
-        table = read_table(self.h5file, TRIGGER_TABLE, start=start, stop=stop)
-        if keep_order:
-            self._add_index_if_needed(table)
+        table = _empty_subarray_events_table()
+
+        if self.load_trigger:
+            trigger = read_table(self.h5file, TRIGGER_TABLE, start=start, stop=stop)
+            table = _join_subarray_events(table, trigger)
 
         if self.load_simulated and SHOWER_TABLE in self.h5file:
             showers = read_table(self.h5file, SHOWER_TABLE, start=start, stop=stop)
@@ -325,9 +297,6 @@ class TableLoader(Component):
                             dl2, prefix=algorithm, ignore=SUBARRAY_EVENT_KEYS
                         )
                         table = _join_subarray_events(table, dl2)
-
-        if keep_order:
-            self._sort_to_original_order(table)
         return table
 
     def read_subarray_events_chunked(self, chunk_size):
@@ -339,9 +308,11 @@ class TableLoader(Component):
         chunk_size: int
             Number of subarray events to load per chunk
         """
+        n_events = self.h5file.root[TRIGGER_TABLE].shape[0]
+
         return ChunkIterator(
             self.read_subarray_events,
-            n_total=len(self),
+            n_total=n_events,
             chunk_size=chunk_size,
         )
 
@@ -420,62 +391,34 @@ class TableLoader(Component):
             )
 
         if self.load_simulated and TRUE_IMPACT_GROUP in self.h5file.root:
-            impacts = self._read_telescope_table(
-                TRUE_IMPACT_GROUP,
-                tel_id,
-                start=start,
-                stop=stop,
-            )
+            impacts = self._read_telescope_table(TRUE_IMPACT_GROUP, tel_id)
             table = _join_telescope_events(table, impacts)
 
         return table
 
     def _read_telescope_events_for_ids(self, tel_ids, tel_start=None, tel_stop=None):
-        tel_start = tel_start if tel_start is not None else [None] * len(tel_ids)
-        tel_stop = tel_stop if tel_stop is not None else [None] * len(tel_ids)
+        if tel_start is None:
+            tel_start = [None] * len(tel_ids)
 
-        tables = []
-        for tel_id, start, stop in zip(tel_ids, tel_start, tel_stop):
-            # no events for this telescope in chunk
-            if start is not None and stop is not None and (stop - start) == 0:
-                continue
+        if tel_stop is None:
+            tel_stop = [None] * len(tel_ids)
 
-            tables.append(
+        table = vstack(
+            [
                 self._read_telescope_events_for_id(tel_id, start=start, stop=stop)
-            )
+                for tel_id, start, stop in zip(tel_ids, tel_start, tel_stop)
+            ]
+        )
 
-        return vstack(tables)
+        return table
 
     def _join_subarray_info(self, table, start=None, stop=None):
-        subarray_events = self.read_subarray_events(
-            start=start,
-            stop=stop,
-            keep_order=False,
-        )
+        subarray_events = self.read_subarray_events(start=start, stop=stop)
         table = join_allow_empty(
             table, subarray_events, keys=SUBARRAY_EVENT_KEYS, join_type="left"
         )
         return table
 
-    def _get_tel_start_stop(self, tel_ids, start, stop):
-        tel_start = None
-        tel_stop = None
-        if start is not None or stop is not None:
-
-            indices = self.subarray.tel_ids_to_indices(tel_ids)
-
-            # find first/last row for each telescope
-            if start is not None:
-                tel_start = self._n_telescope_events[start][indices]
-
-            if stop is not None:
-                if stop >= len(self._n_telescope_events):
-                    tel_stop = None
-                else:
-                    tel_stop = self._n_telescope_events[stop][indices]
-
-        return tel_start, tel_stop
-
     def read_telescope_events(self, telescopes=None, start=None, stop=None):
         """
         Read telescope-based event information.
@@ -510,15 +453,26 @@ class TableLoader(Component):
         else:
             tel_ids = self.subarray.get_tel_ids(telescopes)
 
-        tel_start, tel_stop = self._get_tel_start_stop(tel_ids, start, stop)
+        tel_start = None
+        tel_stop = None
+        if start is not None or stop is not None:
+
+            indices = self.subarray.tel_ids_to_indices(tel_ids)
+
+            # find first/last row for each telescope
+            if start is not None:
+                tel_start = self._n_telescope_events[start][indices]
+
+            if stop is not None:
+                if stop >= len(self._n_telescope_events):
+                    tel_stop = None
+                else:
+                    tel_stop = self._n_telescope_events[stop][indices]
+
         table = self._read_telescope_events_for_ids(tel_ids, tel_start, tel_stop)
-        table = self._join_subarray_info(table, start=start, stop=stop)
 
-        # sort back to order in the file
-        table = _join_subarray_events(
-            table, self._get_sort_index(start=start, stop=stop)
-        )
-        self._sort_to_original_order(table, include_tel_id=True)
+        if any([self.load_trigger, self.load_simulated, self.load_dl2]):
+            table = self._join_subarray_info(table, start=start, stop=stop)
 
         return table
 
@@ -535,9 +489,11 @@ class TableLoader(Component):
 
         **kwargs are passed to `read_telescope_events`
         """
+        n_events = self.h5file.root[TRIGGER_TABLE].shape[0]
+
         return ChunkIterator(
             self.read_telescope_events,
-            n_total=len(self),
+            n_total=n_events,
             chunk_size=chunk_size,
             **kwargs,
         )
@@ -561,9 +517,7 @@ class TableLoader(Component):
         np.cumsum(tels_with_trigger, out=tels_with_trigger, axis=0)
         return tels_with_trigger
 
-    def read_telescope_events_by_type(
-        self, telescopes=None, start=None, stop=None
-    ) -> Dict[str, Table]:
+    def read_telescope_events_by_type(self, telescopes=None) -> Dict[str, Table]:
         """Read telescope-based event information.
 
         Parameters
@@ -579,54 +533,20 @@ class TableLoader(Component):
         """
 
         if telescopes is None:
-            tel_ids = tuple(self.subarray.tel.keys())
+            tel_ids = self.subarray.tel.keys()
         else:
             tel_ids = self.subarray.get_tel_ids(telescopes)
 
-        tel_start, tel_stop = self._get_tel_start_stop(tel_ids, start, stop)
-        tel_start = tel_start if tel_start is not None else [None] * len(tel_ids)
-        tel_stop = tel_stop if tel_stop is not None else [None] * len(tel_ids)
-
         by_type = defaultdict(list)
-        sort_index = self._get_sort_index(start=start, stop=stop)
-
-        for tel_id, start, stop in zip(tel_ids, tel_start, tel_stop):
-            # no events for this telescope in range start/stop
-            if start is not None and stop is not None and (stop - start) == 0:
-                continue
 
+        for tel_id in tel_ids:
             key = str(self.subarray.tel[tel_id])
-            by_type[key].append(
-                self._read_telescope_events_for_id(tel_id, start=start, stop=stop)
-            )
+            by_type[key].append(self._read_telescope_events_for_id(tel_id))
 
         by_type = {k: vstack(ts) for k, ts in by_type.items()}
 
-        for key in by_type.keys():
-            by_type[key] = self._join_subarray_info(
-                by_type[key], start=start, stop=stop
-            )
-            by_type[key] = _join_subarray_events(by_type[key], sort_index)
-            self._sort_to_original_order(by_type[key], include_tel_id=True)
+        if any([self.load_trigger, self.load_simulated, self.load_dl2]):
+            for key, table in by_type.items():
+                by_type[key] = self._join_subarray_info(table)
 
         return by_type
-
-    def read_telescope_events_by_type_chunked(self, chunk_size, **kwargs):
-        """
-        Iterate over chunks of telescope events.
-
-        Parameters
-        ----------
-        chunk_size: int
-            Number of subarray events to load per chunk.
-            The telescope tables might be larger or smaller than chunk_size
-            depending on the selected telescopes.
-
-        **kwargs are passed to `read_telescope_events`
-        """
-        return ChunkIterator(
-            self.read_telescope_events_by_type,
-            n_total=len(self),
-            chunk_size=chunk_size,
-            **kwargs,
-        )


Test patch:
diff --git a/ctapipe/conftest.py b/ctapipe/conftest.py
index 8742faed..4676ca16 100644
--- a/ctapipe/conftest.py
+++ b/ctapipe/conftest.py
@@ -233,30 +233,6 @@ def dl2_shower_geometry_file_type(dl2_tmp_path, prod5_gamma_simtel_path):
         return output
 
 
-@pytest.fixture(scope="session")
-def dl2_merged_file(dl2_tmp_path, dl2_shower_geometry_file, dl2_proton_geometry_file):
-    """
-    File containing both parameters and shower geometry from a gamma simulation set.
-    """
-    from ctapipe.core import run_tool
-    from ctapipe.tools.merge import MergeTool
-
-    output = dl2_tmp_path / "merged.training.h5"
-
-    # prevent running process multiple times in case of parallel tests
-    with FileLock(output.with_suffix(output.suffix + ".lock")):
-        if output.is_file():
-            return output
-
-        argv = [
-            f"--output={output}",
-            str(dl2_proton_geometry_file),
-            str(dl2_shower_geometry_file),
-        ]
-        assert run_tool(MergeTool(), argv=argv, cwd=dl2_tmp_path) == 0
-        return output
-
-
 @pytest.fixture(scope="session")
 def dl1_file(dl1_tmp_path, prod5_gamma_simtel_path):
     """
diff --git a/ctapipe/instrument/tests/test_subarray.py b/ctapipe/instrument/tests/test_subarray.py
index c594730c..1124b5f7 100644
--- a/ctapipe/instrument/tests/test_subarray.py
+++ b/ctapipe/instrument/tests/test_subarray.py
@@ -1,5 +1,6 @@
 """ Tests for SubarrayDescriptions """
 from copy import deepcopy
+from astropy.coordinates.earth import EarthLocation
 
 import numpy as np
 from astropy import units as u
@@ -28,7 +29,16 @@ def example_subarray(n_tels=10):
         )
         pos[tel_id] = rng.uniform(-100, 100, size=3) * u.m
 
-    return SubarrayDescription("test array", tel_positions=pos, tel_descriptions=tel)
+    return SubarrayDescription(
+        "test array",
+        tel_positions=pos,
+        tel_descriptions=tel,
+        reference_location=EarthLocation(
+            lon=-17 * u.deg,
+            lat=28 * u.deg,
+            height=2200 * u.m,
+        ),
+    )
 
 
 def test_subarray_description():
diff --git a/ctapipe/io/tests/test_table_loader.py b/ctapipe/io/tests/test_table_loader.py
index b8a63ff8..b4701be2 100644
--- a/ctapipe/io/tests/test_table_loader.py
+++ b/ctapipe/io/tests/test_table_loader.py
@@ -1,67 +1,6 @@
 import pytest
 import tables
 import numpy as np
-from astropy.table import Table
-
-from ctapipe.io.astropy_helpers import read_table
-
-
-def check_equal_array_event_order(table1, table2):
-    """
-    Check that values and order of array events is consistent in two tables.
-
-    Works for both tables of subarray events (obs_id, tel_id) and
-    tables of telescope events (obs_id, tel_id, event_id) and combinations
-    of the two.
-    """
-
-    def unique_events(table):
-        return np.unique(np.array(table[["obs_id", "event_id"]]), return_index=True)
-
-    unique_events1, indicies1 = unique_events(table1)
-    unique_events2, indicies2 = unique_events(table2)
-
-    if len(unique_events1) != len(unique_events2):
-        raise ValueError("Tables have different numbers of events")
-
-    if np.any(unique_events1 != unique_events2):
-        raise ValueError("Tables have different subarray events")
-
-    # we expect the rank in which the events appear in indicies to be the
-    # same, that means argsort should produce the same result on both:
-    order1 = np.argsort(indicies1)
-    order2 = np.argsort(indicies2)
-
-    if np.any(order1 != order2):
-        raise ValueError("Tables have subarray events in different order")
-
-
-def test_check_order():
-    with pytest.raises(ValueError, match="Tables have different numbers"):
-        check_equal_array_event_order(
-            Table({"obs_id": [1, 1, 2, 2, 3], "event_id": [1, 2, 1, 2, 1]}),
-            Table({"obs_id": [1, 1, 2, 2], "event_id": [1, 2, 1, 2]}),
-        )
-
-    with pytest.raises(ValueError, match="Tables have different subarray events"):
-        check_equal_array_event_order(
-            Table({"obs_id": [1, 1, 2, 2, 3], "event_id": [1, 2, 1, 2, 1]}),
-            Table({"obs_id": [1, 1, 2, 2, 4], "event_id": [1, 2, 1, 2, 1]}),
-        )
-
-    with pytest.raises(
-        ValueError, match="Tables have subarray events in different order"
-    ):
-        check_equal_array_event_order(
-            Table({"obs_id": [1, 1, 3, 2, 2], "event_id": [1, 2, 1, 1, 2]}),
-            Table({"obs_id": [1, 1, 2, 2, 3], "event_id": [1, 2, 1, 2, 1]}),
-        )
-
-    check_equal_array_event_order(
-        Table({"obs_id": [1, 1, 3, 2, 2], "event_id": [1, 2, 1, 1, 2]}),
-        Table({"obs_id": [1, 1, 3, 2, 2], "event_id": [1, 2, 1, 1, 2]}),
-    )
-
 
 from ctapipe.instrument.subarray import SubarrayDescription
 
@@ -104,7 +43,7 @@ def test_telescope_events_for_tel_id(test_file):
 
     _, dl1_file = test_file
 
-    loader = TableLoader(dl1_file, load_dl1_parameters=True)
+    loader = TableLoader(dl1_file, load_dl1_parameters=True, load_trigger=True)
 
     with loader as table_loader:
         table = table_loader.read_telescope_events([8])
@@ -179,36 +118,20 @@ def test_true_parameters(test_file):
 
 def test_read_subarray_events(test_file_dl2):
     """Test reading subarray events"""
+
     from ctapipe.io.tableloader import TableLoader
 
     _, dl2_file = test_file_dl2
 
-    with TableLoader(dl2_file, load_dl2=True, load_simulated=True) as table_loader:
+    with TableLoader(
+        dl2_file, load_dl2=True, load_simulated=True, load_trigger=True
+    ) as table_loader:
         table = table_loader.read_subarray_events()
         assert "HillasReconstructor_alt" in table.colnames
         assert "true_energy" in table.colnames
         assert "time" in table.colnames
 
 
-def test_table_loader_keeps_original_order(dl2_merged_file):
-    """Test reading subarray events keeps order in file"""
-    from ctapipe.io.tableloader import TableLoader
-
-    # check that the order is the same as in the file itself
-    trigger = read_table(dl2_merged_file, "/dl1/event/subarray/trigger")
-    # check we actually have unsorted input
-    assert not np.all(np.diff(trigger["obs_id"]) >= 0)
-
-    with TableLoader(
-        dl2_merged_file, load_dl2=True, load_simulated=True
-    ) as table_loader:
-        events = table_loader.read_subarray_events()
-        tel_events = table_loader.read_telescope_events()
-
-    check_equal_array_event_order(events, trigger)
-    check_equal_array_event_order(events, tel_events)
-
-
 def test_read_telescope_events_type(test_file_dl2):
     """Test reading telescope events for a given telescope type"""
 
@@ -225,6 +148,7 @@ def test_read_telescope_events_type(test_file_dl2):
         load_dl2=True,
         load_simulated=True,
         load_true_images=True,
+        load_trigger=False,
         load_instrument=True,
     ) as table_loader:
 
@@ -254,6 +178,7 @@ def test_read_telescope_events_by_type(test_file_dl2):
         load_dl2=True,
         load_simulated=True,
         load_true_images=True,
+        load_trigger=False,
         load_instrument=True,
     ) as table_loader:
 
@@ -300,8 +225,6 @@ def test_chunked(dl2_shower_geometry_file):
 
     n_chunks = 2
     chunk_size = int(np.ceil(n_events / n_chunks))
-    start = 0
-    stop = chunk_size
 
     with TableLoader(
         dl2_shower_geometry_file,
@@ -310,18 +233,14 @@ def test_chunked(dl2_shower_geometry_file):
         load_dl1_parameters=True,
         load_dl2=True,
         load_simulated=True,
+        load_trigger=True,
     ) as table_loader:
 
         tel_event_it = table_loader.read_telescope_events_chunked(chunk_size)
         event_it = table_loader.read_subarray_events_chunked(chunk_size)
-        by_type_it = table_loader.read_telescope_events_by_type_chunked(chunk_size)
-
-        iters = (event_it, tel_event_it, by_type_it)
 
-        for chunk, (events, tel_events, by_type) in enumerate(zip(*iters)):
+        for chunk, (events, tel_events) in enumerate(zip(event_it, tel_event_it)):
             n_read += len(events)
-            start = chunk * chunk_size
-            stop = min(n_events, (chunk + 1) * chunk_size)
 
             # last chunk might be smaller
             if chunk == (n_chunks - 1):
@@ -329,17 +248,13 @@ def test_chunked(dl2_shower_geometry_file):
             else:
                 assert len(events) == chunk_size
 
-            # check events are in compatible order
-            check_equal_array_event_order(events, tel_events)
-            check_equal_array_event_order(trigger[start:stop], events)
-
-            # check number of telescope events is correct
-            assert len(tel_events) == np.count_nonzero(events["tels_with_trigger"])
+            assert len(tel_events) == np.sum(events["tels_with_trigger"])
 
-            n_events_by_type = 0
-            for table in by_type.values():
-                n_events_by_type += len(table)
+            unique_events_from_tel_table = np.unique(tel_events[["obs_id", "event_id"]])
 
-            assert n_events_by_type == len(tel_events)
+            # test we get matching obs_ids / event_ids
+            idx = events[["obs_id", "event_id"]].argsort()
+            sorted_index = np.array(events[["obs_id", "event_id"]][idx])
+            assert (unique_events_from_tel_table == sorted_index).all()
 
     assert n_read == n_events