Merge pull request #1911 from cta-observatory/feature/better_array_di…

…splay

improve array display

ctapipe/coordinates/__init__.py

@@ -10,7 +10,12 @@
 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 @@
     "NominalFrame",
     "GroundFrame",
     "TiltedGroundFrame",
+    "EastingNorthingFrame",
     "MissingFrameAttributeWarning",
     "project_to_ground",
 ]

ctapipe/coordinates/ground_frames.py

@@ -14,21 +14,25 @@
 - 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

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()

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 @@ def __repr__(self):
 
     @property
     def tel(self):
-        """ for backward compatibility"""
+        """for backward compatibility"""
         return self.tels
 
     @property
@@ -132,7 +130,7 @@ def info(self, printer=print):
 
     @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 @@ def tel_coords(self):
 
     @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 @@ def to_table(self, kind="subarray"):
         meta = {
             "ORIGIN": "ctapipe.instrument.SubarrayDescription",
             "SUBARRAY": self.name,
-            "SOFT_VER": ctapipe.__version__,
+            "SOFT_VER": CTAPIPE_VERSION,
             "TAB_TYPE": kind,
         }
 
@@ -277,8 +275,8 @@ def to_table(self, kind="subarray"):
             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 @@ def peek(self):
         """
         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):