Repository: cta-observatory/ctapipe
Pull request #1911 (2022-05-19 15:08:37)

improve array display
* Adds a `EastingNorthingFrame` to properly display array on ground (since the Corsika/Simtel coordinates are confusing with north going on the X-axis).
* Improve ArrayDisplay to have correct frame attribute names on the x and y axes, and add a radial grid
* use an `EastingNorthingFrame` `ArrayDisplay` in `SubarrayDescription.peek()`

Looks like this:

<img width="654" alt="image" src="https://user-images.githubusercontent.com/11677812/169329269-9fb21961-bbec-4ad0-b72b-eaf5af18cbde.png">



Patch:
diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py
index 5a210b55..d32680ed 100644
--- a/ctapipe/coordinates/__init__.py
+++ b/ctapipe/coordinates/__init__.py
@@ -10,7 +10,12 @@ from astropy.coordinates import (
 import warnings
 from .telescope_frame import TelescopeFrame
 from .nominal_frame import NominalFrame
-from .ground_frames import GroundFrame, TiltedGroundFrame, project_to_ground
+from .ground_frames import (
+    GroundFrame,
+    TiltedGroundFrame,
+    project_to_ground,
+    EastingNorthingFrame,
+)
 from .camera_frame import CameraFrame, EngineeringCameraFrame
 
 
@@ -21,6 +26,7 @@ __all__ = [
     "NominalFrame",
     "GroundFrame",
     "TiltedGroundFrame",
+    "EastingNorthingFrame",
     "MissingFrameAttributeWarning",
     "project_to_ground",
 ]
diff --git a/ctapipe/coordinates/ground_frames.py b/ctapipe/coordinates/ground_frames.py
index 7285607b..8b19bb4d 100644
--- a/ctapipe/coordinates/ground_frames.py
+++ b/ctapipe/coordinates/ground_frames.py
@@ -14,21 +14,25 @@ TODO:
 - Tests Tests Tests!
 """
 import astropy.units as u
+from astropy.units.quantity import Quantity
 import numpy as np
 from astropy.coordinates import (
+    AltAz,
     BaseCoordinateFrame,
     CartesianRepresentation,
-    FunctionTransform,
     CoordinateAttribute,
-    AltAz,
+    FunctionTransform,
+    RepresentationMapping,
+    frame_transform_graph,
+    AffineTransform,
 )
-from astropy.coordinates import frame_transform_graph
 from numpy import cos, sin
 
 __all__ = [
     "GroundFrame",
     "TiltedGroundFrame",
     "project_to_ground",
+    "EastingNorthingFrame",
 ]
 
 
@@ -47,6 +51,25 @@ class GroundFrame(BaseCoordinateFrame):
     default_representation = CartesianRepresentation
 
 
+class EastingNorthingFrame(BaseCoordinateFrame):
+    """GroundFrame but in standard Easting/Northing coordinates instead of
+    SimTel/Corsika conventions
+
+    Frame attributes: None
+
+    """
+
+    default_representation = CartesianRepresentation
+
+    frame_specific_representation_info = {
+        CartesianRepresentation: [
+            RepresentationMapping("x", "easting"),
+            RepresentationMapping("y", "northing"),
+            RepresentationMapping("z", "height"),
+        ]
+    }
+
+
 class TiltedGroundFrame(BaseCoordinateFrame):
     """Tilted ground coordinate frame.  The tilted ground coordinate frame
     is a cartesian system describing the 2 dimensional projected
@@ -200,3 +223,34 @@ def project_to_ground(tilt_system):
     y_projected = y_initial - trans[2][1] * z_initial / trans[2][2]
 
     return GroundFrame(x=x_projected * unit, y=y_projected * unit, z=0 * unit)
+
+
+@frame_transform_graph.transform(FunctionTransform, GroundFrame, GroundFrame)
+def ground_to_ground(ground_coords, ground_frame):
+    """Null transform for going from ground to ground, since there are no
+    attributes of the GroundSystem"""
+    return ground_coords
+
+
+# Matrices for transforming between GroundFrame and EastingNorthingFrame
+NO_OFFSET = CartesianRepresentation(Quantity([0, 0, 0], u.m))
+GROUND_TO_EASTNORTH = np.asarray([[0, -1, 0], [1, 0, 0], [0, 0, 1]])
+
+
+@frame_transform_graph.transform(AffineTransform, GroundFrame, EastingNorthingFrame)
+def ground_to_easting_northing(ground_coords, eastnorth_frame):
+    """
+    convert GroundFrame points into eastings/northings for plotting purposes
+
+    """
+
+    return GROUND_TO_EASTNORTH, NO_OFFSET
+
+
+@frame_transform_graph.transform(AffineTransform, EastingNorthingFrame, GroundFrame)
+def easting_northing_to_ground(eastnorth_coords, ground_frame):
+    """
+    convert  eastings/northings back to GroundFrame
+
+    """
+    return GROUND_TO_EASTNORTH.T, NO_OFFSET
diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py
index 557040b4..cc7c5b6d 100644
--- a/ctapipe/instrument/subarray.py
+++ b/ctapipe/instrument/subarray.py
@@ -1,26 +1,24 @@
 """
 Description of Arrays or Subarrays of telescopes
 """
-from typing import Dict, List, Union
-from contextlib import ExitStack
 import warnings
+from contextlib import ExitStack
+from copy import copy
+from itertools import groupby
+from typing import Dict, List, Union
 
 import numpy as np
+import tables
 from astropy import units as u
 from astropy.coordinates import SkyCoord
 from astropy.table import QTable, Table
 from astropy.utils import lazyproperty
-import tables
-from copy import copy
-from itertools import groupby
 
-import ctapipe
-
-from ..coordinates import GroundFrame, CameraFrame
-from .telescope import TelescopeDescription
-from .camera import CameraDescription, CameraReadout, CameraGeometry
+from .. import __version__ as CTAPIPE_VERSION
+from ..coordinates import CameraFrame, GroundFrame
+from .camera import CameraDescription, CameraGeometry, CameraReadout
 from .optics import OpticsDescription
-
+from .telescope import TelescopeDescription
 
 __all__ = ["SubarrayDescription"]
 
@@ -90,7 +88,7 @@ class SubarrayDescription:
 
     @property
     def tel(self):
-        """ for backward compatibility"""
+        """for backward compatibility"""
         return self.tels
 
     @property
@@ -132,7 +130,7 @@ class SubarrayDescription:
 
     @lazyproperty
     def tel_coords(self):
-        """ returns telescope positions as astropy.coordinates.SkyCoord"""
+        """returns telescope positions as astropy.coordinates.SkyCoord"""
 
         pos_x = np.array([p[0].to("m").value for p in self.positions.values()]) * u.m
         pos_y = np.array([p[1].to("m").value for p in self.positions.values()]) * u.m
@@ -142,7 +140,7 @@ class SubarrayDescription:
 
     @lazyproperty
     def tel_ids(self):
-        """ telescope IDs as an array"""
+        """telescope IDs as an array"""
         return np.array(list(self.tel.keys()))
 
     @lazyproperty
@@ -238,7 +236,7 @@ class SubarrayDescription:
         meta = {
             "ORIGIN": "ctapipe.instrument.SubarrayDescription",
             "SUBARRAY": self.name,
-            "SOFT_VER": ctapipe.__version__,
+            "SOFT_VER": CTAPIPE_VERSION,
             "TAB_TYPE": kind,
         }
 
@@ -277,8 +275,8 @@ class SubarrayDescription:
             unique_types = self.telescope_types
 
             mirror_area = u.Quantity(
-                [t.optics.mirror_area.to_value(u.m ** 2) for t in unique_types],
-                u.m ** 2,
+                [t.optics.mirror_area.to_value(u.m**2) for t in unique_types],
+                u.m**2,
             )
             focal_length = u.Quantity(
                 [t.optics.equivalent_focal_length.to_value(u.m) for t in unique_types],
@@ -332,47 +330,29 @@ class SubarrayDescription:
         """
         Draw a quick matplotlib plot of the array
         """
+        from ctapipe.coordinates.ground_frames import EastingNorthingFrame
+        from ctapipe.visualization import ArrayDisplay
         from matplotlib import pyplot as plt
-        from astropy.visualization import quantity_support
-
-        types = set(self.tels.values())
-        tab = self.to_table()
 
         plt.figure(figsize=(8, 8))
-
-        with quantity_support():
-            for tel_type in types:
-                tels = tab[tab["tel_description"] == str(tel_type)]["tel_id"]
-                sub = self.select_subarray(tels, name=tel_type)
-                tel_coords = sub.tel_coords
-                radius = np.array(
-                    [
-                        np.sqrt(tel.optics.mirror_area / np.pi).value
-                        for tel in sub.tels.values()
-                    ]
-                )
-
-                plt.scatter(
-                    tel_coords.x, tel_coords.y, s=radius * 8, alpha=0.5, label=tel_type
-                )
-
-            plt.legend(loc="best")
-            plt.title(self.name)
-            plt.tight_layout()
+        ad = ArrayDisplay(subarray=self, frame=EastingNorthingFrame(), tel_scale=0.75)
+        ad.add_labels()
+        plt.title(self.name)
+        plt.tight_layout()
 
     @lazyproperty
     def telescope_types(self) -> List[TelescopeDescription]:
-        """ list of telescope types in the array"""
+        """list of telescope types in the array"""
         return list({tel for tel in self.tel.values()})
 
     @lazyproperty
     def camera_types(self) -> List[CameraDescription]:
-        """ list of camera types in the array """
+        """list of camera types in the array"""
         return list({tel.camera for tel in self.tel.values()})
 
     @lazyproperty
     def optics_types(self) -> List[OpticsDescription]:
-        """ list of optics types in the array """
+        """list of optics types in the array"""
         return list({tel.optics for tel in self.tel.values()})
 
     def get_tel_ids_for_type(self, tel_type):
diff --git a/ctapipe/visualization/mpl_array.py b/ctapipe/visualization/mpl_array.py
index 5d895ec8..1c3b3e2c 100644
--- a/ctapipe/visualization/mpl_array.py
+++ b/ctapipe/visualization/mpl_array.py
@@ -62,10 +62,12 @@ class ArrayDisplay:
 
         self.frame = frame
         self.subarray = subarray
+        self.axes = axes or plt.gca()
 
         # get the telescope positions. If a new frame is set, this will
         # transform to the new frame.
-        self.tel_coords = subarray.tel_coords.transform_to(frame)
+        self.tel_coords = subarray.tel_coords.transform_to(frame).cartesian
+        self.unit = self.tel_coords.x.unit
 
         # set up colors per telescope type
         tel_types = [str(tel) for tel in subarray.tels.values()]
@@ -76,6 +78,10 @@ class ArrayDisplay:
                 for tel in subarray.tel.values()
             ]
 
+            self.radii = radius
+        else:
+            self.radii = np.ones(len(tel_types)) * radius
+
         if title is None:
             title = subarray.name
 
@@ -116,15 +122,22 @@ class ArrayDisplay:
             )
         plt.legend(handles=legend_elements)
 
+        self.add_radial_grid()
+
+        # create the plot
         self.tel_colors = tel_color
         self.autoupdate = autoupdate
         self.telescopes = PatchCollection(patches, match_original=True)
         self.telescopes.set_linewidth(2.0)
 
-        self.axes = axes or plt.gca()
         self.axes.add_collection(self.telescopes)
         self.axes.set_aspect(1.0)
         self.axes.set_title(title)
+        xunit = self.tel_coords.x.unit.to_string("latex")
+        yunit = self.tel_coords.y.unit.to_string("latex")
+        xname, yname, _ = frame.get_representation_component_names().keys()
+        self.axes.set_xlabel(f"{xname} [{xunit}] $\\rightarrow$")
+        self.axes.set_ylabel(f"{yname} [{yunit}] $\\rightarrow$")
         self._labels = []
         self._quiver = None
         self.axes.autoscale_view()
@@ -136,12 +149,50 @@ class ArrayDisplay:
 
     @values.setter
     def values(self, values):
-        """ set the telescope colors to display  """
+        """set the telescope colors to display"""
         self.telescopes.set_array(np.ma.masked_invalid(values))
         self._update()
 
+    def add_radial_grid(self, spacing=100 * u.m):
+        """add some dotted rings for distance estimation. The number of rings
+        is estimated automatically from the spacing and the array footprint.
+
+        Parameters
+        ----------
+        spacing: Quantity
+            spacing between rings
+
+        """
+
+        n_circles = np.round(
+            (np.sqrt(self.subarray.footprint / np.pi) / spacing).to_value(""),
+            0,
+        )
+        circle_radii = np.arange(1, n_circles + 2, 1) * spacing.to_value(self.unit)
+        circle_patches = PatchCollection(
+            [
+                Circle(
+                    xy=(0, 0),
+                    radius=r,
+                    fill=False,
+                    fc="none",
+                    linestyle="dotted",
+                    color="gray",
+                    alpha=0.1,
+                    lw=1,
+                )
+                for r in circle_radii
+            ],
+            color="#eeeeee",
+            ls="dotted",
+            fc="none",
+            lw=3,
+        )
+
+        self.axes.add_collection(circle_patches)
+
     def set_vector_uv(self, uu, vv, c=None, **kwargs):
-        """ sets the vector field U,V and color for all telescopes
+        """sets the vector field U,V and color for all telescopes
 
         Parameters
         ----------
@@ -289,9 +340,17 @@ class ArrayDisplay:
     def add_labels(self):
         px = self.tel_coords.x.to_value("m")
         py = self.tel_coords.y.to_value("m")
-        for tel, x, y in zip(self.subarray.tels, px, py):
+        for tel, x, y, r in zip(self.subarray.tels, px, py, self.radii):
             name = str(tel)
-            lab = self.axes.text(x, y, name, fontsize=8, clip_on=True)
+            lab = self.axes.text(
+                x,
+                y - r * 1.8,
+                name,
+                fontsize=8,
+                clip_on=True,
+                horizontalalignment="center",
+                verticalalignment="top",
+            )
             self._labels.append(lab)
 
     def remove_labels(self):
@@ -300,7 +359,7 @@ class ArrayDisplay:
         self._labels = []
 
     def _update(self):
-        """ signal a redraw if necessary """
+        """signal a redraw if necessary"""
         if self.autoupdate:
             plt.draw()
 


Test patch:
diff --git a/ctapipe/coordinates/tests/test_coordinates.py b/ctapipe/coordinates/tests/test_coordinates.py
index 23206d67..ef75a58d 100644
--- a/ctapipe/coordinates/tests/test_coordinates.py
+++ b/ctapipe/coordinates/tests/test_coordinates.py
@@ -257,6 +257,7 @@ def test_camera_focal_length_array():
 
 
 def test_ground_frame_roundtrip():
+    """test transform from sky to ground roundtrip"""
     from ctapipe.coordinates import GroundFrame, TiltedGroundFrame
 
     normal = SkyCoord(alt=70 * u.deg, az=0 * u.deg, frame=AltAz())
@@ -268,3 +269,22 @@ def test_ground_frame_roundtrip():
     assert u.isclose(coord.x, back.x, atol=1e-12 * u.m)
     assert u.isclose(coord.y, back.y, atol=1e-12 * u.m)
     assert u.isclose(coord.z, back.z, atol=1e-12 * u.m)
+
+
+def test_ground_to_eastnorth_roundtrip():
+    """Check Ground to EastingNorthing and the round-trip"""
+    from ctapipe.coordinates import GroundFrame, EastingNorthingFrame
+
+    ground = SkyCoord(
+        x=[1, 2, 3] * u.m, y=[-2, 5, 2] * u.m, z=[1, -1, 2] * u.m, frame=GroundFrame()
+    )
+    eastnorth = ground.transform_to(EastingNorthingFrame())
+    ground2 = eastnorth.transform_to(GroundFrame())
+
+    assert u.isclose(eastnorth.easting, [2, -5, -2] * u.m).all()
+    assert u.isclose(eastnorth.northing, [1, 2, 3] * u.m).all()
+    assert u.isclose(eastnorth.height, [1, -1, 2] * u.m).all()
+
+    assert u.isclose(ground.x, ground2.x).all()
+    assert u.isclose(ground.y, ground2.y).all()
+    assert u.isclose(ground.z, ground2.z).all()
\ No newline at end of file