Repository: cta-observatory/ctapipe (python)

Issue:
Error when `SubarrayDescription.select_subarray` is given unknown telescope ids could be better
When using more than one allowed telescope in your config file and the pedestal dataset only contains data for a specific telescope id you get the following, imprecise error message:
```
ERROR [ctapipe.ctapipe-process] (tool.run): Caught unexpected exception: 2
Traceback (most recent call last):
File "/ctasoft/ctapipe/ctapipe/core/tool.py", line 336, in run
    self.setup()
  File "/ctasoft/ctapipe/ctapipe/tools/process.py", line 150, in setup
    self.event_source = EventSource(parent=self)
  File "/ctasoft/ctapipe/ctapipe/io/simteleventsource.py", line 243, in __init__
    self._subarray_info = self.prepare_subarray_info(
  File "/ctasoft/ctapipe/ctapipe/io/simteleventsource.py", line 364, in prepare_subarray_info
    subarray = subarray.select_subarray(self.allowed_tels)
  File "/ctasoft/ctapipe/ctapipe/instrument/subarray.py", line 319, in select_subarray
    tel_positions = {tid: self.positions[tid] for tid in tel_ids}
  File "/ctasoft/ctapipe/ctapipe/instrument/subarray.py", line 319, in <dictcomp>
    tel_positions = {tid: self.positions[tid] for tid in tel_ids}
KeyError: 2
```

An error message telling the user what telescope id's are allowed would be helpful for future use.


Pull request:
Better error message in case unknown telescopes are selected, fixes #1924

Patch:
diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py
index d32680ed..5a210b55 100644
--- a/ctapipe/coordinates/__init__.py
+++ b/ctapipe/coordinates/__init__.py
@@ -10,12 +10,7 @@ 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,
-    EastingNorthingFrame,
-)
+from .ground_frames import GroundFrame, TiltedGroundFrame, project_to_ground
 from .camera_frame import CameraFrame, EngineeringCameraFrame
 
 
@@ -26,7 +21,6 @@ __all__ = [
     "NominalFrame",
     "GroundFrame",
     "TiltedGroundFrame",
-    "EastingNorthingFrame",
     "MissingFrameAttributeWarning",
     "project_to_ground",
 ]
diff --git a/ctapipe/coordinates/ground_frames.py b/ctapipe/coordinates/ground_frames.py
index 8b19bb4d..7285607b 100644
--- a/ctapipe/coordinates/ground_frames.py
+++ b/ctapipe/coordinates/ground_frames.py
@@ -14,25 +14,21 @@ 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,
-    CoordinateAttribute,
     FunctionTransform,
-    RepresentationMapping,
-    frame_transform_graph,
-    AffineTransform,
+    CoordinateAttribute,
+    AltAz,
 )
+from astropy.coordinates import frame_transform_graph
 from numpy import cos, sin
 
 __all__ = [
     "GroundFrame",
     "TiltedGroundFrame",
     "project_to_ground",
-    "EastingNorthingFrame",
 ]
 
 
@@ -51,25 +47,6 @@ 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
@@ -223,34 +200,3 @@ 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/__init__.py b/ctapipe/instrument/__init__.py
index 79442154..561e992e 100644
--- a/ctapipe/instrument/__init__.py
+++ b/ctapipe/instrument/__init__.py
@@ -2,7 +2,7 @@ from .camera import CameraDescription, CameraGeometry, CameraReadout, PixelShape
 from .atmosphere import get_atmosphere_profile_functions
 from .telescope import TelescopeDescription
 from .optics import OpticsDescription
-from .subarray import SubarrayDescription
+from .subarray import SubarrayDescription, UnknownTelescopeID
 from .guess import guess_telescope
 
 
@@ -11,9 +11,10 @@ __all__ = [
     "CameraGeometry",
     "CameraReadout",
     "get_atmosphere_profile_functions",
-    "TelescopeDescription",
-    "OpticsDescription",
-    "SubarrayDescription",
     "guess_telescope",
+    "OpticsDescription",
     "PixelShape",
+    "SubarrayDescription",
+    "TelescopeDescription",
+    "UnknownTelescopeID",
 ]
diff --git a/ctapipe/instrument/subarray.py b/ctapipe/instrument/subarray.py
index cc7c5b6d..1ae60d23 100644
--- a/ctapipe/instrument/subarray.py
+++ b/ctapipe/instrument/subarray.py
@@ -1,34 +1,41 @@
 """
 Description of Arrays or Subarrays of telescopes
 """
-import warnings
-from contextlib import ExitStack
-from copy import copy
-from itertools import groupby
 from typing import Dict, List, Union
+from contextlib import ExitStack
+import warnings
 
 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
 
-from .. import __version__ as CTAPIPE_VERSION
-from ..coordinates import CameraFrame, GroundFrame
-from .camera import CameraDescription, CameraGeometry, CameraReadout
-from .optics import OpticsDescription
+import ctapipe
+
+from ..coordinates import GroundFrame, CameraFrame
 from .telescope import TelescopeDescription
+from .camera import CameraDescription, CameraReadout, CameraGeometry
+from .optics import OpticsDescription
+
 
 __all__ = ["SubarrayDescription"]
 
 
+class UnknownTelescopeID(KeyError):
+    """Raised when an unknown telescope id is encountered"""
+
+
 def _group_consecutives(sequence):
     """
     Turn consequtive lists into ranges (used in SubarrayDescription.info())
 
     from https://codereview.stackexchange.com/questions/214820/codewars-range-extraction
     """
+    sequence = sorted(sequence)
     for _, g in groupby(enumerate(sequence), lambda i_x: i_x[0] - i_x[1]):
         r = [x for _, x in g]
         if len(r) > 2:
@@ -236,7 +243,7 @@ class SubarrayDescription:
         meta = {
             "ORIGIN": "ctapipe.instrument.SubarrayDescription",
             "SUBARRAY": self.name,
-            "SOFT_VER": CTAPIPE_VERSION,
+            "SOFT_VER": ctapipe.__version__,
             "TAB_TYPE": kind,
         }
 
@@ -314,11 +321,15 @@ class SubarrayDescription:
         SubarrayDescription
         """
 
+        unknown_tel_ids = set(tel_ids).difference(self.tel.keys())
+        if len(unknown_tel_ids) > 0:
+            known = _range_extraction(self.tel.keys())
+            raise UnknownTelescopeID(f"{unknown_tel_ids}, known telescopes: {known}")
+
         tel_positions = {tid: self.positions[tid] for tid in tel_ids}
         tel_descriptions = {tid: self.tel[tid] for tid in tel_ids}
 
         if not name:
-            tel_ids = sorted(tel_ids)
             name = self.name + "_" + _range_extraction(tel_ids)
 
         newsub = SubarrayDescription(
@@ -330,15 +341,33 @@ 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))
-        ad = ArrayDisplay(subarray=self, frame=EastingNorthingFrame(), tel_scale=0.75)
-        ad.add_labels()
-        plt.title(self.name)
-        plt.tight_layout()
+
+        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()
 
     @lazyproperty
     def telescope_types(self) -> List[TelescopeDescription]:
diff --git a/ctapipe/visualization/mpl_array.py b/ctapipe/visualization/mpl_array.py
index 1c3b3e2c..5d895ec8 100644
--- a/ctapipe/visualization/mpl_array.py
+++ b/ctapipe/visualization/mpl_array.py
@@ -62,12 +62,10 @@ 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).cartesian
-        self.unit = self.tel_coords.x.unit
+        self.tel_coords = subarray.tel_coords.transform_to(frame)
 
         # set up colors per telescope type
         tel_types = [str(tel) for tel in subarray.tels.values()]
@@ -78,10 +76,6 @@ 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
 
@@ -122,22 +116,15 @@ 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()
@@ -149,50 +136,12 @@ 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
         ----------
@@ -340,17 +289,9 @@ 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, r in zip(self.subarray.tels, px, py, self.radii):
+        for tel, x, y in zip(self.subarray.tels, px, py):
             name = str(tel)
-            lab = self.axes.text(
-                x,
-                y - r * 1.8,
-                name,
-                fontsize=8,
-                clip_on=True,
-                horizontalalignment="center",
-                verticalalignment="top",
-            )
+            lab = self.axes.text(x, y, name, fontsize=8, clip_on=True)
             self._labels.append(lab)
 
     def remove_labels(self):
@@ -359,7 +300,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 ef75a58d..23206d67 100644
--- a/ctapipe/coordinates/tests/test_coordinates.py
+++ b/ctapipe/coordinates/tests/test_coordinates.py
@@ -257,7 +257,6 @@ 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())
@@ -269,22 +268,3 @@ 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
diff --git a/ctapipe/instrument/tests/test_subarray.py b/ctapipe/instrument/tests/test_subarray.py
index 333957cf..c594730c 100644
--- a/ctapipe/instrument/tests/test_subarray.py
+++ b/ctapipe/instrument/tests/test_subarray.py
@@ -16,7 +16,7 @@ from ctapipe.instrument import (
 
 
 def example_subarray(n_tels=10):
-    """ generate a simple subarray for testing purposes """
+    """generate a simple subarray for testing purposes"""
     rng = np.random.default_rng(0)
 
     pos = {}
@@ -32,7 +32,7 @@ def example_subarray(n_tels=10):
 
 
 def test_subarray_description():
-    """ Test SubarrayDescription functionality """
+    """Test SubarrayDescription functionality"""
     n_tels = 10
     sub = example_subarray(n_tels)
     sub.peek()
@@ -67,14 +67,14 @@ def test_subarray_description():
 
 
 def test_to_table(example_subarray):
-    """ Check that we can generate astropy Tables from the SubarrayDescription """
+    """Check that we can generate astropy Tables from the SubarrayDescription"""
     sub = example_subarray
     assert len(sub.to_table(kind="subarray")) == sub.num_tels
     assert len(sub.to_table(kind="optics")) == len(sub.optics_types)
 
 
 def test_tel_indexing(example_subarray):
-    """ Check that we can convert between telescope_id and telescope_index """
+    """Check that we can convert between telescope_id and telescope_index"""
     sub = example_subarray
 
     assert sub.tel_indices[1] == 0  # first tel_id is in slot 0
@@ -218,3 +218,10 @@ def test_get_tel_ids(example_subarray):
     # test invalid telescope type
     with pytest.raises(Exception):
         tel_ids = subarray.get_tel_ids(["It's a-me, Mario!"])
+
+
+def test_unknown_telescopes(example_subarray):
+    from ctapipe.instrument import UnknownTelescopeID
+
+    with pytest.raises(UnknownTelescopeID):
+        example_subarray.select_subarray([300, 201])