Repository: cta-observatory/ctapipe
Pull request #2537 (2024-04-22 11:46:45)

Add function to get point on shower axis in altaz
This can be used to compute the shower axis in telescope / camera frame, e.g. like this:

```python
import astropy.units as u
from astropy.coordinates import AltAz
import matplotlib.pyplot as plt

from ctapipe.calib import CameraCalibrator
from ctapipe.io import EventSource
from ctapipe.coordinates import get_point_on_shower_axis, TelescopeFrame
from ctapipe.visualization import CameraDisplay


plt.style.use("dark_background")
plt.rcParams["savefig.facecolor"] = "0.15"

path = "dataset://lst_prod3_calibration_and_mcphotons.simtel.zst"

if __name__ == "__main__":
    with EventSource(path, focal_length_choice="EQUIVALENT") as source:
        it = iter(source)
        for _ in range(2):
            e = next(it)
        subarray = source.subarray

    calib = CameraCalibrator(subarray)
    calib(e)

    shower = e.simulation.shower
    distance = 5 * u.km

    fig, axs = plt.subplots(2, 2, layout="constrained", figsize=(8, 8))

    telescope_positions = source.subarray.tel_coords
    axis_point = get_point_on_shower_axis(
        core_x=shower.core_x,
        core_y=shower.core_y,
        alt=shower.alt,
        az=shower.az,
        telescope_position=telescope_positions,
        distance=distance,
    )

    for ax, (tel_id, dl1) in zip(axs.flat, e.dl1.tel.items()):
        tel_index = source.subarray.tel_ids_to_indices(tel_id)[0]

        pointing = AltAz(
            alt=e.pointing.tel[tel_id].altitude, az=e.pointing.tel[tel_id].azimuth
        )
        tel_frame = TelescopeFrame(telescope_pointing=pointing)

        geometry = source.subarray.tel[tel_id].camera.geometry
        disp = CameraDisplay(geometry.transform_to(tel_frame), image=dl1.image, ax=ax)

        ax.set_title(f"LST-{tel_id}")
        ax.set_axis_off()
        ax.set_ylabel("")
        ax.set_xlabel("")

        source_tel_frame = AltAz(alt=shower.alt, az=shower.az).transform_to(tel_frame)
        p = axis_point[tel_index].transform_to(tel_frame)

        disp.axes.plot(
            [source_tel_frame.fov_lon.deg, p.fov_lon.deg],
            [source_tel_frame.fov_lat.deg, p.fov_lat.deg],
            color="xkcd:baby blue",
            marker=".",
        )
        disp.overlay_coordinate(source_tel_frame, color="xkcd:light yellow")

    fig.suptitle(
        f"$E={{}}${shower.energy.to(u.GeV):.0f}, Impact = ({shower.core_x:.0f}, {shower.core_y:.0f})"
    )

    plt.savefig("shower_axis.png", dpi=200)
```

![shower_axis](https://github.com/cta-observatory/ctapipe/assets/5488440/37e5adf9-ce83-4fe1-bb6a-972aed4c0cd0)



ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 0 Issues
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 0 Code Smells

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/100.svg?sanitize=true) 100.00% Coverage (92.80% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.70% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2537)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 0 Issues
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 0 Code Smells

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/100.svg?sanitize=true) 100.00% Coverage (92.80% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.70% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2537)


ctao-dpps-sonarqube[bot]: ![Passed](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/QualityGateBadge/passed.svg?sanitize=true)

# Analysis Details
## 0 Issues
- ![Bug](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/bug.svg?sanitize=true) 0 Bugs
- ![Vulnerability](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/vulnerability.svg?sanitize=true) 0 Vulnerabilities
- ![Code Smell](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/common/code_smell.svg?sanitize=true) 0 Code Smells

## Coverage and Duplications
- ![Coverage](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/CoverageChart/100.svg?sanitize=true) 100.00% Coverage (92.80% Estimated after merge)
- ![Duplications](https://sonar-cta-dpps.zeuthen.desy.de/static/communityBranchPlugin/checks/Duplications/3.svg?sanitize=true) 0.00% Duplicated Code (0.70% Estimated after merge)

**Project ID:** cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs

[View in SonarQube](https://sonar-cta-dpps.zeuthen.desy.de/dashboard?id=cta-observatory_ctapipe_AY52EYhuvuGcMFidNyUs&pullRequest=2537)


maxnoe: > Looks good to me. Can you show an example with larger images?

Here is a 5 TeV shower:

![Uploading impact_5tev.png…]()


Patch:
diff --git a/docs/changes/2537.feature.rst b/docs/changes/2537.feature.rst
new file mode 100644
index 00000000..2804cc27
--- /dev/null
+++ b/docs/changes/2537.feature.rst
@@ -0,0 +1,3 @@
+Add function ``ctapipe.coordinates.get_point_on_shower_axis``
+that computes a point on the shower axis in alt/az as seen
+from a telescope.
diff --git a/src/ctapipe/coordinates/__init__.py b/src/ctapipe/coordinates/__init__.py
index 388dd3e4..f975b41e 100644
--- a/src/ctapipe/coordinates/__init__.py
+++ b/src/ctapipe/coordinates/__init__.py
@@ -20,7 +20,7 @@ from .ground_frames import (
 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
+from .utils import altaz_to_righthanded_cartesian, get_point_on_shower_axis
 
 __all__ = [
     "TelescopeFrame",
@@ -35,6 +35,7 @@ __all__ = [
     "altaz_to_righthanded_cartesian",
     "impact_distance",
     "shower_impact_distance",
+    "get_point_on_shower_axis",
 ]
 
 
diff --git a/src/ctapipe/coordinates/ground_frames.py b/src/ctapipe/coordinates/ground_frames.py
index e0b8b8c6..77a6e8a5 100644
--- a/src/ctapipe/coordinates/ground_frames.py
+++ b/src/ctapipe/coordinates/ground_frames.py
@@ -46,11 +46,15 @@ def _earthlocation_to_altaz(location, reference_location):
 
 
 class GroundFrame(BaseCoordinateFrame):
-    """Ground coordinate frame.  The ground coordinate frame is a simple
-    cartesian frame describing the 3 dimensional position of objects
-    compared to the array ground level in relation to the nomial
-    centre of the array.  Typically this frame will be used for
-    describing the position on telescopes and equipment.
+    """Ground coordinate frame.
+
+    The ground coordinate frame is a simple cartesian frame describing the
+    3 dimensional position of objects compared to the array ground level
+    in relation to the nominal center of the array.
+
+    Typically this frame will be used for describing the position of telescopes,
+    equipment and shower impact coordinates.
+
     In this frame x points north, y points west and z is meters above array center.
 
     Frame attributes: None
diff --git a/src/ctapipe/coordinates/utils.py b/src/ctapipe/coordinates/utils.py
index 3d46c64c..33c766ae 100644
--- a/src/ctapipe/coordinates/utils.py
+++ b/src/ctapipe/coordinates/utils.py
@@ -1,12 +1,21 @@
 import astropy.units as u
+import numpy as np
+from astropy.coordinates import AltAz
+from erfa.ufunc import p2s as cartesian_to_spherical
 from erfa.ufunc import s2p as spherical_to_cartesian
 
+from .ground_frames import _get_xyz
+
 __all__ = [
     "altaz_to_righthanded_cartesian",
+    "get_point_on_shower_axis",
 ]
 
 
-def altaz_to_righthanded_cartesian(alt, az):
+_zero_m = u.Quantity(0, u.m)
+
+
+def altaz_to_righthanded_cartesian(alt, az, distance=1.0):
     """Turns an alt/az coordinate into a 3d direction in a right-handed coordinate
     system.  This is because azimuths are in a left-handed system.
 
@@ -19,8 +28,55 @@ def altaz_to_righthanded_cartesian(alt, az):
     az: u.Quantity
         azimuth
     """
+    # this is an optimization, shaves of ca. 50 % compared to handing over units
+    # to the erfa function
     if hasattr(az, "unit"):
         az = az.to_value(u.rad)
         alt = alt.to_value(u.rad)
 
-    return spherical_to_cartesian(-az, alt, 1.0)
+    unit = None
+    if hasattr(distance, "unit"):
+        unit = distance.unit
+        distance = distance.value
+
+    result = spherical_to_cartesian(-az, alt, distance)
+    if unit is not None:
+        return u.Quantity(result, unit=unit, copy=False)
+    return result
+
+
+@u.quantity_input(core_x=u.m, core_y=u.m, alt=u.rad, az=u.rad, distance=u.km)
+def get_point_on_shower_axis(core_x, core_y, alt, az, telescope_position, distance):
+    """
+    Get a point on the shower axis in AltAz frame as seen by a telescope at the given position.
+
+    Parameters
+    ----------
+    core_x : u.Quantity[length]
+        Impact x-coordinate
+    core_y : u.Quantity[length]
+        Impact y-coordinate
+    alt : u.Quantity[angle]
+        Altitude of primary
+    az : u.Quantity[angle]
+        Azimuth of primary
+    telescope_position : GroundFrame
+        Telescope position
+    distance : u.Quantity[length]
+        Distance along the shower axis from the ground of the point returned.
+
+    Returns
+    -------
+    coord : AltAz
+        The AltAz coordinate of a point on the shower axis at the given distance
+        from the impact point.
+    """
+    impact = u.Quantity([core_x, core_y, _zero_m], unit=u.m)
+    # move up the shower axis by slant_distance
+    point = impact + altaz_to_righthanded_cartesian(alt=alt, az=az, distance=distance)
+
+    # offset by telescope positions and convert to sperical
+    # to get local AltAz for each telescope
+    cartesian = point[np.newaxis, :] - _get_xyz(telescope_position).T
+    lon, lat, _ = cartesian_to_spherical(cartesian)
+    return AltAz(alt=lat, az=-lon, copy=False)
diff --git a/src/ctapipe/io/simteleventsource.py b/src/ctapipe/io/simteleventsource.py
index 4cda9ae6..133aeeae 100644
--- a/src/ctapipe/io/simteleventsource.py
+++ b/src/ctapipe/io/simteleventsource.py
@@ -630,7 +630,7 @@ class SimTelEventSource(EventSource):
 
                 # TODO: switch to warning or even an exception once
                 # we start relying on this.
-                self.log.info(
+                self.log.debug(
                     "Could not determine telescope from sim_telarray metadata,"
                     " guessing using builtin lookup-table: %d: %s",
                     tel_id,


Test patch:
diff --git a/src/ctapipe/coordinates/tests/test_utils.py b/src/ctapipe/coordinates/tests/test_utils.py
new file mode 100644
index 00000000..1ebd419c
--- /dev/null
+++ b/src/ctapipe/coordinates/tests/test_utils.py
@@ -0,0 +1,55 @@
+import astropy.units as u
+import numpy as np
+import pytest
+from astropy.coordinates import AltAz
+
+
+def test_point_on_shower_axis_far(subarray_prod5_paranal):
+    """Test for get_point_on_shower_axis"""
+    from ctapipe.coordinates import get_point_on_shower_axis
+
+    core_x = 50 * u.m
+    core_y = 100 * u.m
+    alt = 68 * u.deg
+    az = 30 * u.deg
+    # for a very large distance, should be identical to the shower direction
+    distance = 1000 * u.km
+
+    point = get_point_on_shower_axis(
+        core_x=core_x,
+        core_y=core_y,
+        alt=alt,
+        az=az,
+        telescope_position=subarray_prod5_paranal.tel_coords,
+        distance=distance,
+    )
+
+    np.testing.assert_allclose(point.alt, alt, rtol=1e-3)
+    np.testing.assert_allclose(point.az, az, rtol=1e-2)
+
+
+def test_single_telescope(subarray_prod5_paranal):
+    from ctapipe.coordinates import (
+        MissingFrameAttributeWarning,
+        get_point_on_shower_axis,
+    )
+
+    core_x = 50 * u.m
+    core_y = 100 * u.m
+    alt = 68 * u.deg
+    az = 30 * u.deg
+    distance = 10 * u.km
+
+    point = get_point_on_shower_axis(
+        core_x=core_x,
+        core_y=core_y,
+        alt=alt,
+        az=az,
+        telescope_position=subarray_prod5_paranal.tel_coords[0],
+        distance=distance,
+    )
+
+    source = AltAz(alt=alt, az=az)
+    # 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)