Refactor DL2EventLoader to use FeatureGenerator

docs/changes/2919.bugfix.rst

@@ -0,0 +1 @@
+Removed the option ``output_table_schema`` from ``DL2EventLoader``, which was incorrectly implemented.

src/ctapipe/conftest.py

@@ -985,37 +985,25 @@ def irf_events_table():
 
 
 @pytest.fixture(scope="function")
-def test_config():
-    return {
-        "DL2EventLoader": {"event_reader_function": "read_telescope_events_chunked"},
-        "DL2EventPreprocessor": {
-            "energy_reconstructor": "ExtraTreesRegressor",
-            "gammaness_classifier": "ExtraTreesClassifier",
-            "columns_to_rename": {},
-            "output_table_schema": [
-                Column(
-                    name="obs_id", dtype=np.uint64, description="Observation Block ID"
-                ),
-                Column(name="event_id", dtype=np.uint64, description="Array event ID"),
-                Column(name="tel_id", dtype=np.uint64, description="Telescope ID"),
-                Column(
-                    name="ExtraTreesRegressor_tel_energy",
-                    unit=u.TeV,
-                    description="Reconstructed energy",
-                ),
-                Column(
-                    name="ExtraTreesRegressor_tel_energy_uncert",
-                    unit=u.TeV,
-                    description="Reconstructed energy uncertainty",
-                ),
-            ],
-            "apply_derived_columns": False,
-            # "disable_column_renaming": True,
-            "allow_unsupported_pointing_frames": True,
-        },
-        "DL2EventQualityQuery": {
-            "quality_criteria": [
-                ("valid reco", "ExtraTreesRegressor_tel_is_valid"),
-            ]
-        },
-    }
+def dl2_event_loader_config():
+    """A traitlets Config for the DL2EventLoader class."""
+    from traitlets.config import Config
+
+    return Config(
+        {
+            "DL2EventLoader": {
+                "event_reader_function": "read_telescope_events_chunked"
+            },
+            "DL2EventPreprocessor": {
+                "energy_reconstructor": "ExtraTreesRegressor",
+                "gammaness_classifier": "ExtraTreesClassifier",
+                "apply_derived_columns": False,
+                "allow_unsupported_pointing_frames": True,
+            },
+            "DL2EventQualityQuery": {
+                "quality_criteria": [
+                    ("valid reco", "ExtraTreesRegressor_tel_is_valid"),
+                ]
+            },
+        }
+    )

src/ctapipe/coordinates/__init__.py

@@ -21,7 +21,11 @@
 from .impact_distance import impact_distance, shower_impact_distance
 from .nominal_frame import NominalFrame
 from .telescope_frame import TelescopeFrame
-from .utils import altaz_to_righthanded_cartesian, get_point_on_shower_axis
+from .utils import (
+    altaz_to_fov,
+    altaz_to_righthanded_cartesian,
+    get_point_on_shower_axis,
+)
 
 __all__ = [
     "TelescopeFrame",
@@ -37,6 +41,7 @@
     "impact_distance",
     "shower_impact_distance",
     "get_point_on_shower_axis",
+    "altaz_to_fov",
 ]
 
 

src/ctapipe/coordinates/tests/test_utils.py

@@ -53,3 +53,18 @@ def test_single_telescope(subarray_prod5_paranal):
     # 10 km is around the shower maximum, should be around 1 degree from the source
     with pytest.warns(MissingFrameAttributeWarning):
         assert u.isclose(source.separation(point), 1.0 * u.deg, atol=0.1 * u.deg)
+
+
+def test_altaz_to_fov():
+    from ctapipe.coordinates import altaz_to_fov
+
+    column = altaz_to_fov(
+        az=[220.0, 220.2] * u.deg,
+        alt=[80.0, 79.2] * u.deg,
+        pointing_az=[220.0, 220.0] * u.deg,
+        pointing_alt=[80.0, 80.0] * u.deg,
+    )
+
+    assert column.unit == u.deg
+    assert np.allclose(column[0].value, 0)
+    assert np.allclose(column[1].value, [-0.03747984, -0.79993558])

src/ctapipe/coordinates/utils.py

@@ -1,14 +1,16 @@
 import astropy.units as u
 import numpy as np
-from astropy.coordinates import AltAz
+from astropy.coordinates import AltAz, SkyCoord
 from erfa.ufunc import p2s as cartesian_to_spherical
 from erfa.ufunc import s2p as spherical_to_cartesian
 
 from .ground_frames import _get_xyz
+from .nominal_frame import NominalFrame
 
 __all__ = [
     "altaz_to_righthanded_cartesian",
     "get_point_on_shower_axis",
+    "altaz_to_fov",
 ]
 
 
@@ -80,3 +82,25 @@ def get_point_on_shower_axis(core_x, core_y, alt, az, telescope_position, distan
     cartesian = point[np.newaxis, :] - _get_xyz(telescope_position).T
     lon, lat, _ = cartesian_to_spherical(cartesian)
     return AltAz(alt=lat, az=-lon, copy=False)
+
+
+def altaz_to_fov(az, alt, pointing_az, pointing_alt) -> u.Quantity[2]:
+    """
+    Compute FOV coordinates from alt/az coordinates.
+
+    This can be used in an FeatureGenerator or ExpressionEngine to get a single
+    column with fov_lon, fov_lat coordinates.
+
+    Returns
+    -------
+    u.Quantity[2]:
+       2D array of coordinates with 2 columns: fov_lon, fov_lat
+    """
+    pointing_coord = SkyCoord(az=pointing_az, alt=pointing_alt, frame="altaz")
+    event_coord = SkyCoord(az=az, alt=alt, frame="altaz", origin=pointing_coord)
+    nominal_coord = event_coord.transform_to(NominalFrame)
+
+    # note the minus sign for the fov_lon, which is to match GADF conventions
+    return u.Quantity(
+        np.column_stack((-nominal_coord.fov_lon.deg, nominal_coord.fov_lat.deg)), u.deg
+    )

src/ctapipe/core/feature_generator.py

@@ -23,7 +23,13 @@ def _shallow_copy_table(table):
     the original table.
     """
     # automatically return Table or QTable depending on input
-    return table.__class__({col: table[col] for col in table.colnames}, copy=False)
+    table_copy = table.__class__(
+        {col: table[col] for col in table.colnames}, copy=False
+    )
+    table_copy.meta = (
+        table.meta.copy()
+    )  # here we don't want a shallow copy, as we might add metadata
+    return table_copy
 
 
 class FeatureGeneratorException(TypeError):

src/ctapipe/core/tests/test_feature_generator.py

@@ -19,19 +19,26 @@ def test_feature_generator():
     generator = FeatureGenerator(features=expressions)
 
     table = Table({"intensity": [1, 10, 100], "length": [2, 4, 8], "width": [1, 2, 4]})
+    table.meta["VERSION"] = 1.0
+
     log_intensity = [0, 1, 2]
     area = [2, 8, 32]
     eccentricity = np.sqrt(0.75)
 
-    table = generator(table)
+    processed_table = generator(table)
+    processed_table.meta["GEN"] = True
+
+    assert "log_intensity" in processed_table.colnames
+    assert "area" in processed_table.colnames
+    assert "eccentricity" in processed_table.colnames
 
-    assert "log_intensity" in table.colnames
-    assert "area" in table.colnames
-    assert "eccentricity" in table.colnames
+    assert np.all(processed_table["log_intensity"] == log_intensity)
+    assert np.all(processed_table["area"] == area)
+    assert np.all(processed_table["eccentricity"] == eccentricity)
 
-    assert np.all(table["log_intensity"] == log_intensity)
-    assert np.all(table["area"] == area)
-    assert np.all(table["eccentricity"] == eccentricity)
+    # also check that metadata was propagated and doesn't overwrite old table:
+    assert "VERSION" in processed_table.meta
+    assert "GEN" not in table.meta
 
 
 def test_existing_feature():

src/ctapipe/io/__init__.py

@@ -7,7 +7,7 @@
 
 from .astropy_helpers import read_table, write_table  # noqa: I001
 from .datalevels import DataLevel
-from .dl2_tables_preprocessing import DL2EventPreprocessor, DL2EventLoader
+from .event_preprocessor import EventPreprocessor, PreprocessorFeatureSet
 from .eventsource import EventSource
 from .eventseeker import EventSeeker
 from .tableio import TableReader, TableWriter
@@ -43,7 +43,7 @@
     "DataWriter",
     "DATA_MODEL_VERSION",
     "get_hdf5_datalevels",
-    "DL2EventPreprocessor",
-    "DL2EventLoader",
+    "EventPreprocessor",
     "get_hdf5_monitoring_types",
+    "PreprocessorFeatureSet",
 ]