diff --git a/docs/changes/2409.feature.rst b/docs/changes/2409.feature.rst new file mode 100644 index 00000000000..7d5af104a79 --- /dev/null +++ b/docs/changes/2409.feature.rst @@ -0,0 +1,5 @@ +Add support for interpolating a monitoring pointing table +in ``TableLoader``. The corresponding table is not yet written by ``ctapipe``, +but can be written by external tools. +This is to enable analysis of real observations, where the pointing changes over time in +alt/az. diff --git a/src/ctapipe/conftest.py b/src/ctapipe/conftest.py index 42d53fc0e39..4d743e5d217 100644 --- a/src/ctapipe/conftest.py +++ b/src/ctapipe/conftest.py @@ -1,12 +1,15 @@ """ common pytest fixtures for tests in ctapipe """ - +import shutil from copy import deepcopy import astropy.units as u +import numpy as np import pytest +import tables from astropy.coordinates import EarthLocation +from astropy.table import Table from pytest_astropy_header.display import PYTEST_HEADER_MODULES from ctapipe.core import run_tool @@ -646,3 +649,38 @@ def disp_reconstructor_path(model_tmp_path, gamma_train_clf): def reference_location(): """a dummy EarthLocation to use for SubarrayDescriptions""" return EarthLocation(lon=-17 * u.deg, lat=28 * u.deg, height=2200 * u.m) + + +@pytest.fixture(scope="session") +def dl1_mon_pointing_file(dl1_file, dl1_tmp_path): + from ctapipe.instrument import SubarrayDescription + from ctapipe.io import read_table, write_table + + path = dl1_tmp_path / "dl1_mon_ponting.dl1.h5" + shutil.copy(dl1_file, path) + + events = read_table(path, "/dl1/event/subarray/trigger") + subarray = SubarrayDescription.from_hdf(path) + + # create some dummy monitoring data + time = events["time"] + start, stop = time[[0, -1]] + duration = (stop - start).to_value(u.s) + + # start a bit before, go a bit longer + dt = np.arange(-1, duration + 2, 1) * u.s + time_mon = start + dt + + alt = (69 + 2 * dt / dt[-1]) * u.deg + az = (180 + 5 * dt / dt[-1]) * u.deg + + table = Table({"time": time_mon, "azimuth": az, "altitude": alt}) + + for tel_id in subarray.tel: + write_table(table, path, f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}") + + # remove static pointing table + with tables.open_file(path, "r+") as f: + f.remove_node("/configuration/telescope/pointing", recursive=True) + + return path diff --git a/src/ctapipe/io/hdf5eventsource.py b/src/ctapipe/io/hdf5eventsource.py index c270a612f12..a451bde4116 100644 --- a/src/ctapipe/io/hdf5eventsource.py +++ b/src/ctapipe/io/hdf5eventsource.py @@ -52,7 +52,8 @@ from .datalevels import DataLevel from .eventsource import EventSource from .hdf5tableio import HDF5TableReader -from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP +from .pointing import PointingInterpolator +from .tableloader import DL2_SUBARRAY_GROUP, DL2_TELESCOPE_GROUP, POINTING_GROUP __all__ = ["HDF5EventSource"] @@ -581,6 +582,13 @@ def _generator(self): ignore_columns={"trigger_pixels"}, ) + pointing_interpolator = None + if POINTING_GROUP in self.file_.root: + pointing_interpolator = PointingInterpolator( + h5file=self.file_, + parent=self, + ) + counter = 0 for trigger, index in events: data = ArrayEventContainer( @@ -630,7 +638,7 @@ def _generator(self): continue self._fill_array_pointing(data) - self._fill_telescope_pointing(data) + self._fill_telescope_pointing(data, pointing_interpolator) for tel_id in data.trigger.tel.keys(): key = f"tel_{tel_id:03d}" @@ -715,10 +723,19 @@ def _fill_array_pointing(self, event): else: raise ValueError(f"Unsupported pointing frame: {frame}") - def _fill_telescope_pointing(self, event): + def _fill_telescope_pointing(self, event, tel_pointing_interpolator=None): """ Fill the telescope pointing information of a given event """ + if tel_pointing_interpolator is not None: + for tel_id, trigger in event.trigger.tel.items(): + alt, az = tel_pointing_interpolator(tel_id, trigger.time) + event.pointing.tel[tel_id] = TelescopePointingContainer( + altitude=alt, + azimuth=az, + ) + return + for tel_id in event.trigger.tels_with_trigger: tel_pointing = self._constant_telescope_pointing.get(tel_id) if tel_pointing is None: diff --git a/src/ctapipe/io/pointing.py b/src/ctapipe/io/pointing.py new file mode 100644 index 00000000000..412c0510a98 --- /dev/null +++ b/src/ctapipe/io/pointing.py @@ -0,0 +1,137 @@ +from typing import Any + +import astropy.units as u +import numpy as np +import tables +from astropy.time import Time +from scipy.interpolate import interp1d + +from ctapipe.core import Component, traits + +from .astropy_helpers import read_table + + +class PointingInterpolator(Component): + """ + Interpolate pointing from a monitoring table to a given timestamp. + + Parameters + ---------- + h5file : None | tables.File + A open hdf5 file with read access. + The monitoring table is expected to be stored in that file at + ``/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}`` + + If not given, monitoring tables can be added via `PointingInterpolator.add_table`. + """ + + bounds_error = traits.Bool( + default_value=True, + help="If true, raises an exception when trying to extrapolate out of the given table", + ).tag(config=True) + + extrapolate = traits.Bool( + help="If bounds_error is False, this flag will specify whether values outside" + "the available values are filled with nan (False) or extrapolated (True).", + default_value=False, + ).tag(config=True) + + def __init__(self, h5file=None, **kwargs): + super().__init__(**kwargs) + + if h5file is not None and not isinstance(h5file, tables.File): + raise TypeError("h5file must be a tables.File") + self.h5file = h5file + + self.interp_options: dict[str, Any] = dict(assume_sorted=True, copy=False) + if self.bounds_error: + self.interp_options["bounds_error"] = True + elif self.extrapolate: + self.interp_options["bounds_error"] = False + self.interp_options["fill_value"] = "extrapolate" + else: + self.interp_options["bounds_error"] = False + self.interp_options["fill_value"] = np.nan + + self._alt_interpolators = {} + self._az_interpolators = {} + + def add_table(self, tel_id, pointing_table): + """ + Add a table to this interpolator + + Parameters + ---------- + tel_id : int + Telescope id + pointing_table : astropy.table.Table + Table of pointing values, expected columns + are ``time`` as ``Time`` column, ``azimuth`` and ``altitude`` + as quantity columns. + """ + missing = {"time", "azimuth", "altitude"} - set(pointing_table.colnames) + if len(missing) > 0: + raise ValueError(f"Table is missing required column(s): {missing}") + + if not isinstance(pointing_table["time"], Time): + raise TypeError("'time' column of pointing table must be astropy.time.Time") + + for col in ("azimuth", "altitude"): + unit = pointing_table[col].unit + if unit is None or not u.rad.is_equivalent(unit): + raise ValueError(f"{col} must have units compatible with 'rad'") + + # sort first, so it's not done twice for each interpolator + pointing_table.sort("time") + # interpolate in mjd TAI. Float64 mjd is precise enough for pointing + # and TAI is contiguous, so no issues with leap seconds. + mjd = pointing_table["time"].tai.mjd + + az = pointing_table["azimuth"].quantity.to_value(u.rad) + # prepare azimuth for interpolation by "unwrapping": i.e. turning + # [359, 1] into [359, 361]. This assumes that if we get values like + # [359, 1] the telescope moved 2 degrees through 0, not 358 degrees + # the other way around. This should be true for all telescopes given + # the sampling speed of pointing values and their maximum movement speed. + # No telescope can turn more than 180° in 2 seconds. + az = np.unwrap(az) + alt = pointing_table["altitude"].quantity.to_value(u.rad) + + self._az_interpolators[tel_id] = interp1d(mjd, az, **self.interp_options) + self._alt_interpolators[tel_id] = interp1d(mjd, alt, **self.interp_options) + + def _read_pointing_table(self, tel_id): + pointing_table = read_table( + self.h5file, + f"/dl0/monitoring/telescope/pointing/tel_{tel_id:03d}", + ) + self.add_table(tel_id, pointing_table) + + def __call__(self, tel_id, time): + """ + Interpolate alt/az for given time and tel_id. + + Parameters + ---------- + tel_id : int + telescope id + time : astropy.time.Time + time for which to interpolate the pointing + + Returns + ------- + altitude : astropy.units.Quantity[deg] + interpolated altitude angle + azimuth : astropy.units.Quantity[deg] + interpolated azimuth angle + """ + if tel_id not in self._az_interpolators: + if self.h5file is not None: + self._read_pointing_table(tel_id) + else: + raise KeyError(f"No pointing table available for tel_id {tel_id}") + + mjd = time.tai.mjd + az = u.Quantity(self._az_interpolators[tel_id](mjd), u.rad, copy=False) + alt = u.Quantity(self._alt_interpolators[tel_id](mjd), u.rad, copy=False) + return alt, az diff --git a/src/ctapipe/io/tableloader.py b/src/ctapipe/io/tableloader.py index 74a3903064d..36a954fb7de 100644 --- a/src/ctapipe/io/tableloader.py +++ b/src/ctapipe/io/tableloader.py @@ -11,11 +11,10 @@ from astropy.table import Table, hstack, vstack from astropy.utils.decorators import lazyproperty -from ctapipe.instrument.optics import FocalLengthKind - from ..core import Component, Provenance, traits -from ..instrument import SubarrayDescription +from ..instrument import FocalLengthKind, SubarrayDescription from .astropy_helpers import join_allow_empty, read_table +from .pointing import PointingInterpolator __all__ = ["TableLoader"] @@ -30,7 +29,8 @@ SIMULATION_CONFIG_TABLE = "/configuration/simulation/run" SHOWER_DISTRIBUTION_TABLE = "/simulation/service/shower_distribution" OBSERVATION_TABLE = "/configuration/observation/observation_block" -FIXED_POINTING_GROUP = "configuration/telescope/pointing" +FIXED_POINTING_GROUP = "/configuration/telescope/pointing" +POINTING_GROUP = "/dl0/monitoring/telescope/pointing" DL2_SUBARRAY_GROUP = "/dl2/event/subarray" DL2_TELESCOPE_GROUP = "/dl2/event/telescope" @@ -204,7 +204,7 @@ class TableLoader(Component): pointing = traits.Bool( True, - help="Load pointing information and join to events", + help="Load pointing information and interpolate / join to events", ).tag(config=True) focal_length_choice = traits.UseEnum( @@ -248,6 +248,11 @@ def __init__(self, input_url=None, h5file=None, **kwargs): Provenance().add_input_file(self.input_url, role="Event data") + self._pointing_interpolator = PointingInterpolator( + h5file=self.h5file, + parent=self, + ) + def close(self): """Close the underlying hdf5 file""" if self._should_close: @@ -595,15 +600,21 @@ def _read_telescope_events_for_id( ) table = _join_telescope_events(table, impacts) - if len(table) > 0 and pointing and FIXED_POINTING_GROUP in self.h5file.root: - pointing = read_table( - self.h5file, f"{FIXED_POINTING_GROUP}/tel_{tel_id:03d}" - ) - # we know that we only have the tel_id we are looking for, remove optimize joining - del pointing["tel_id"] - table = join_allow_empty( - table, pointing, ["obs_id"], "left", keep_order=True - ) + if len(table) > 0 and pointing: + # prefer monitoring pointing + if POINTING_GROUP in self.h5file.root: + alt, az = self._pointing_interpolator(tel_id, table["time"]) + table["telescope_pointing_altitude"] = alt + table["telescope_pointing_azimuth"] = az + elif FIXED_POINTING_GROUP in self.h5file.root: + pointing_table = read_table( + self.h5file, f"{FIXED_POINTING_GROUP}/tel_{tel_id:03d}" + ) + # we know that we only have the tel_id we are looking for, remove optimize joining + del pointing_table["tel_id"] + table = join_allow_empty( + table, pointing_table, ["obs_id"], "left", keep_order=True + ) return table diff --git a/src/ctapipe/io/tests/test_hdf5eventsource.py b/src/ctapipe/io/tests/test_hdf5eventsource.py index f426e024413..bc8b1aa33d9 100644 --- a/src/ctapipe/io/tests/test_hdf5eventsource.py +++ b/src/ctapipe/io/tests/test_hdf5eventsource.py @@ -251,3 +251,14 @@ def test_dl1_camera_frame(dl1_camera_frame_file): assert sim.true_parameters.hillas.intensity is not None assert tel_id is not None, "did not test any events" + + +def test_interpolate_pointing(dl1_mon_pointing_file): + from ctapipe.io import HDF5EventSource + + with HDF5EventSource(dl1_mon_pointing_file) as source: + for e in source: + assert set(e.pointing.tel.keys()) == set(e.trigger.tels_with_trigger) + for pointing in e.pointing.tel.values(): + assert not np.isnan(pointing.azimuth) + assert not np.isnan(pointing.altitude) diff --git a/src/ctapipe/io/tests/test_pointing.py b/src/ctapipe/io/tests/test_pointing.py new file mode 100644 index 00000000000..10913053c18 --- /dev/null +++ b/src/ctapipe/io/tests/test_pointing.py @@ -0,0 +1,154 @@ +import astropy.units as u +import numpy as np +import pytest +import tables +from astropy.table import Table +from astropy.time import Time + + +def test_simple(): + """Test pointing interpolation""" + from ctapipe.io.pointing import PointingInterpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, + "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, + "altitude": np.linspace(70.0, 60.0, 6) * u.deg, + }, + ) + + interpolator = PointingInterpolator() + interpolator.add_table(1, table) + + alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) + assert u.isclose(alt, 69 * u.deg) + assert u.isclose(az, 1 * u.deg) + + with pytest.raises(KeyError): + interpolator(tel_id=2, time=t0 + 1 * u.s) + + +def test_azimuth_switchover(): + """Test pointing interpolation""" + from ctapipe.io.pointing import PointingInterpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + [0, 1, 2] * u.s, + "azimuth": [359, 1, 3] * u.deg, + "altitude": [60, 61, 62] * u.deg, + }, + ) + + interpolator = PointingInterpolator() + interpolator.add_table(1, table) + + alt, az = interpolator(tel_id=1, time=t0 + 0.5 * u.s) + assert u.isclose(az, 360 * u.deg) + assert u.isclose(alt, 60.5 * u.deg) + + +def test_invalid_input(): + """Test invalid pointing tables raise nice errors""" + from ctapipe.io.pointing import PointingInterpolator + + wrong_time = Table( + { + "time": [1, 2, 3] * u.s, + "azimuth": [1, 2, 3] * u.deg, + "altitude": [1, 2, 3] * u.deg, + } + ) + + interpolator = PointingInterpolator() + with pytest.raises(TypeError, match="astropy.time.Time"): + interpolator.add_table(1, wrong_time) + + wrong_unit = Table( + { + "time": Time(1.7e9 + np.arange(3), format="unix"), + "azimuth": [1, 2, 3] * u.m, + "altitude": [1, 2, 3] * u.deg, + } + ) + with pytest.raises(ValueError, match="compatible with 'rad'"): + interpolator.add_table(1, wrong_unit) + + wrong_unit = Table( + { + "time": Time(1.7e9 + np.arange(3), format="unix"), + "azimuth": [1, 2, 3] * u.deg, + "altitude": [1, 2, 3], + } + ) + with pytest.raises(ValueError, match="compatible with 'rad'"): + interpolator.add_table(1, wrong_unit) + + +def test_hdf5(tmp_path): + from ctapipe.io import write_table + from ctapipe.io.pointing import PointingInterpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, + "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, + "altitude": np.linspace(70.0, 60.0, 6) * u.deg, + }, + ) + + path = tmp_path / "pointing.h5" + write_table(table, path, "/dl0/monitoring/telescope/pointing/tel_001") + with tables.open_file(path) as h5file: + interpolator = PointingInterpolator(h5file) + alt, az = interpolator(tel_id=1, time=t0 + 1 * u.s) + assert u.isclose(alt, 69 * u.deg) + assert u.isclose(az, 1 * u.deg) + + +def test_bounds(): + """Test invalid pointing tables raise nice errors""" + from ctapipe.io.pointing import PointingInterpolator + + t0 = Time("2022-01-01T00:00:00") + + table = Table( + { + "time": t0 + np.arange(0.0, 10.1, 2.0) * u.s, + "azimuth": np.linspace(0.0, 10.0, 6) * u.deg, + "altitude": np.linspace(70.0, 60.0, 6) * u.deg, + }, + ) + + interpolator = PointingInterpolator() + interpolator.add_table(1, table) + + with pytest.raises(ValueError, match="below the interpolation range"): + interpolator(tel_id=1, time=t0 - 0.1 * u.s) + + with pytest.raises(ValueError, match="above the interpolation range"): + interpolator(tel_id=1, time=t0 + 10.2 * u.s) + + interpolator = PointingInterpolator(bounds_error=False) + interpolator.add_table(1, table) + for dt in (-0.1, 10.1) * u.s: + alt, az = interpolator(tel_id=1, time=t0 + dt) + assert np.isnan(alt.value) + assert np.isnan(az.value) + + interpolator = PointingInterpolator(bounds_error=False, extrapolate=True) + interpolator.add_table(1, table) + alt, az = interpolator(tel_id=1, time=t0 - 1 * u.s) + assert u.isclose(alt, 71 * u.deg) + assert u.isclose(az, -1 * u.deg) + + alt, az = interpolator(tel_id=1, time=t0 + 11 * u.s) + assert u.isclose(alt, 59 * u.deg) + assert u.isclose(az, 11 * u.deg) diff --git a/src/ctapipe/io/tests/test_table_loader.py b/src/ctapipe/io/tests/test_table_loader.py index adfd59b3e6f..398b58d4277 100644 --- a/src/ctapipe/io/tests/test_table_loader.py +++ b/src/ctapipe/io/tests/test_table_loader.py @@ -416,3 +416,13 @@ def test_order_merged(): for tel, table in tables.items(): mask = np.isin(tel_trigger["tel_id"], loader.subarray.get_tel_ids(tel)) check_equal_array_event_order(table, tel_trigger[mask]) + + +def test_interpolate_pointing(dl1_mon_pointing_file, tmp_path): + from ctapipe.io import TableLoader + + with TableLoader(dl1_mon_pointing_file, pointing=True) as loader: + events = loader.read_telescope_events([4]) + assert len(events) > 0 + assert "telescope_pointing_azimuth" in events.colnames + assert "telescope_pointing_altitude" in events.colnames diff --git a/src/ctapipe/reco/sklearn.py b/src/ctapipe/reco/sklearn.py index 24c565c4f03..b9005a462fe 100644 --- a/src/ctapipe/reco/sklearn.py +++ b/src/ctapipe/reco/sklearn.py @@ -823,13 +823,20 @@ def predict_table(self, key, table: Table) -> dict[ReconstructionProperty, Table fov_lon = table["hillas_fov_lon"].quantity + disp * np.cos(psi) fov_lat = table["hillas_fov_lat"].quantity + disp * np.sin(psi) - # FIXME: Assume constant and parallel pointing for each run - self.log.warning("Assuming constant and parallel pointing for each run") + # 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"] + alt, az = telescope_to_horizontal( lon=fov_lon, lat=fov_lat, - pointing_alt=table["subarray_pointing_lat"], - pointing_az=table["subarray_pointing_lon"], + pointing_alt=pointing_alt, + pointing_az=pointing_az, ) altaz_result = Table(