Merge pull request #2347 from cta-observatory/overlay_telescope_frame

Add support for hillas parameters in TelescopeFrame to CameraDisplay

ctapipe/visualization/bokeh.py

@@ -1,31 +1,30 @@
 import sys
 import tempfile
 from abc import ABCMeta
-import matplotlib.pyplot as plt
-from matplotlib.colors import to_hex
 
+import astropy.units as u
+import matplotlib.pyplot as plt
 import numpy as np
-
-from bokeh.io import output_notebook, push_notebook, show, output_file
-from bokeh.plotting import figure
+from bokeh.io import output_file, output_notebook, push_notebook, show
 from bokeh.models import (
-    ColumnDataSource,
-    TapTool,
+    BoxZoomTool,
+    CategoricalColorMapper,
     ColorBar,
-    LinearColorMapper,
-    LogColorMapper,
+    ColumnDataSource,
     ContinuousColorMapper,
-    CategoricalColorMapper,
-    HoverTool,
-    BoxZoomTool,
     Ellipse,
+    HoverTool,
     Label,
+    LinearColorMapper,
+    LogColorMapper,
+    TapTool,
 )
-from bokeh.palettes import Viridis256, Magma256, Inferno256, Greys256, d3
-import astropy.units as u
+from bokeh.palettes import Greys256, Inferno256, Magma256, Viridis256, d3
+from bokeh.plotting import figure
+from matplotlib.colors import to_hex
 
 from ..instrument import CameraGeometry, PixelShape
-
+from .utils import build_hillas_overlay
 
 PLOTARGS = dict(tools="", toolbar_location=None, outline_line_color="#595959")
 
@@ -66,10 +65,11 @@ def generate_hex_vertices(geom):
     phi += geom.pix_rotation.rad + np.deg2rad(30)
 
     # we need the circumcircle radius, pixel_width is incircle diameter
-    r = 2 / np.sqrt(3) * geom.pixel_width.value / 2
+    unit = geom.pix_x.unit
+    r = 2 / np.sqrt(3) * geom.pixel_width.to_value(unit) / 2
 
-    x = geom.pix_x.value
-    y = geom.pix_y.value
+    x = geom.pix_x.to_value(unit)
+    y = geom.pix_y.to_value(unit)
 
     return (
         x[:, np.newaxis] + r[:, np.newaxis] * np.cos(phi)[np.newaxis],
@@ -79,9 +79,10 @@ def generate_hex_vertices(geom):
 
 def generate_square_vertices(geom):
     """Generate vertices of pixels for a square grid camera geometry"""
-    width = geom.pixel_width.value / 2
-    x = geom.pix_x.value
-    y = geom.pix_y.value
+    unit = geom.pix_x.unit
+    width = geom.pixel_width.to_value(unit) / 2
+    x = geom.pix_x.to_value(unit)
+    y = geom.pix_y.to_value(unit)
 
     x_offset = width[:, np.newaxis] * np.array([-1, -1, 1, 1])
     y_offset = width[:, np.newaxis] * np.array([1, -1, -1, 1])
@@ -327,15 +328,18 @@ def _init_datasource(self, image=None):
             line_alpha=np.zeros(self._geometry.n_pixels),
         )
 
+        self._unit = self._geometry.pix_x.unit
+
         if self._geometry.pix_type == PixelShape.HEXAGON:
             x, y = generate_hex_vertices(self._geometry)
 
         elif self._geometry.pix_type == PixelShape.SQUARE:
             x, y = generate_square_vertices(self._geometry)
 
         elif self._geometry.pix_type == PixelShape.CIRCLE:
-            x, y = self._geometry.pix_x.value, self._geometry.pix_y.value
-            data["radius"] = self._geometry.pixel_width / 2
+            x = self._geometry.pix_x.to_value(self._unit)
+            y = self._geometry.pix_y.to_value(self._unit)
+            data["radius"] = self._geometry.pixel_width.to_value(self._unit) / 2
         else:
             raise NotImplementedError(
                 f"Unsupported pixel shape {self._geometry.pix_type}"
@@ -464,7 +468,7 @@ def add_ellipse(self, centroid, length, width, angle, asymmetry=0.0, **kwargs):
         return ellipse
 
     def overlay_moments(
-        self, hillas_parameters, with_label=True, keep_old=False, **kwargs
+        self, hillas_parameters, with_label=True, keep_old=False, n_sigma=1, **kwargs
     ):
         """helper to overlay ellipse from a `HillasParametersContainer` structure
 
@@ -483,30 +487,29 @@ def overlay_moments(
         if not keep_old:
             self.clear_overlays()
 
-        # strip off any units
-        cen_x = u.Quantity(hillas_parameters.x).value
-        cen_y = u.Quantity(hillas_parameters.y).value
-        length = u.Quantity(hillas_parameters.length).value
-        width = u.Quantity(hillas_parameters.width).value
+        params = build_hillas_overlay(
+            hillas_parameters,
+            self._unit,
+            n_sigma=n_sigma,
+            with_label=with_label,
+        )
 
         el = self.add_ellipse(
-            centroid=(cen_x, cen_y),
-            length=length * 2,
-            width=width * 2,
-            angle=hillas_parameters.psi.to_value(u.rad),
+            centroid=(params["cog_x"], params["cog_y"]),
+            length=2 * n_sigma * params["length"],
+            width=2 * n_sigma * params["width"],
+            angle=params["psi_rad"],
             **kwargs,
         )
 
         if with_label:
             label = Label(
-                x=cen_x,
-                y=cen_y,
-                text="({:.02f},{:.02f})\n[w={:.02f},l={:.02f}]".format(
-                    hillas_parameters.x,
-                    hillas_parameters.y,
-                    hillas_parameters.width,
-                    hillas_parameters.length,
-                ),
+                x=params["label_x"],
+                y=params["label_y"],
+                text=params["text"],
+                angle=params["rotation"],
+                angle_units="deg",
+                text_align="center",
                 text_color=el.line_color,
             )
             self.figure.add_layout(label, "center")
@@ -644,7 +647,7 @@ def _init_datasource(self, subarray, values, *, radius, frame, scale, alpha):
         for i, telescope_id in enumerate(telescope_ids):
             telescope = subarray.tel[telescope_id]
             tel_types.append(str(telescope))
-            mirror_area = telescope.optics.mirror_area.to_value(u.m ** 2)
+            mirror_area = telescope.optics.mirror_area.to_value(u.m**2)
             mirror_radii[i] = np.sqrt(mirror_area) / np.pi
 
         if values is None:

ctapipe/visualization/mpl_camera.py

@@ -14,6 +14,7 @@
 
 from ..coordinates import get_representation_component_names
 from ..instrument import PixelShape
+from .utils import build_hillas_overlay
 
 __all__ = ["CameraDisplay"]
 
@@ -459,7 +460,7 @@ def overlay_coordinate(self, coord, keep_old=False, **kwargs):
         self._axes_overlays.append(plot)
 
     def overlay_moments(
-        self, hillas_parameters, with_label=True, keep_old=False, **kwargs
+        self, hillas_parameters, with_label=True, keep_old=False, n_sigma=1, **kwargs
     ):
         """helper to overlay ellipse from a `~ctapipe.containers.HillasParametersContainer` structure
 
@@ -471,40 +472,42 @@ def overlay_moments(
             If True, show coordinates of centroid and width and length
         keep_old: bool
             If True, to not remove old overlays
+        n_sigma: float
+            How many sigmas to use for the ellipse
         kwargs: key=value
             any style keywords to pass to matplotlib (e.g. color='red'
             or linewidth=6)
         """
         if not keep_old:
             self.clear_overlays()
 
-        # strip off any units
-        cen_x = u.Quantity(hillas_parameters.x).to_value(self.unit)
-        cen_y = u.Quantity(hillas_parameters.y).to_value(self.unit)
-        length = u.Quantity(hillas_parameters.length).to_value(self.unit)
-        width = u.Quantity(hillas_parameters.width).to_value(self.unit)
+        params = build_hillas_overlay(
+            hillas_parameters,
+            self.unit,
+            n_sigma=n_sigma,
+            with_label=with_label,
+        )
 
         el = self.add_ellipse(
-            centroid=(cen_x, cen_y),
-            length=length * 2,
-            width=width * 2,
-            angle=hillas_parameters.psi.to_value(u.rad),
+            centroid=(params["cog_x"], params["cog_y"]),
+            length=n_sigma * params["length"] * 2,
+            width=n_sigma * params["width"] * 2,
+            angle=params["psi_rad"],
             **kwargs,
         )
 
         self._axes_overlays.append(el)
 
         if with_label:
             text = self.axes.text(
-                cen_x,
-                cen_y,
-                "({:.02f},{:.02f})\n[w={:.02f},l={:.02f}]".format(
-                    hillas_parameters.x,
-                    hillas_parameters.y,
-                    hillas_parameters.width,
-                    hillas_parameters.length,
-                ),
+                params["label_x"],
+                params["label_y"],
+                params["text"],
                 color=el.get_edgecolor(),
+                va="bottom",
+                ha="center",
+                rotation=params["rotation"],
+                rotation_mode="anchor",
             )
 
             self._axes_overlays.append(text)

ctapipe/visualization/tests/test_mpl.py

@@ -90,7 +90,7 @@ def test_camera_display_single(prod5_lst_cam, tmp_path):
     fig.savefig(tmp_path / "result.png")
 
 
-def test_hillas_overlay(prod5_lst_cam, tmp_path):
+def test_hillas_overlay_camera_frame(prod5_lst_cam, tmp_path):
     from ctapipe.visualization import CameraDisplay
 
     fig, ax = plt.subplots()
@@ -103,6 +103,23 @@ def test_hillas_overlay(prod5_lst_cam, tmp_path):
     fig.savefig(tmp_path / "result.png")
 
 
+def test_hillas_overlay(prod5_lst_cam, tmp_path):
+    from ctapipe.visualization import CameraDisplay
+
+    fig, ax = plt.subplots()
+    disp = CameraDisplay(prod5_lst_cam.transform_to(TelescopeFrame()), ax=ax)
+    hillas = HillasParametersContainer(
+        fov_lon=0.1 * u.deg,
+        fov_lat=-0.1 * u.deg,
+        length=0.5 * u.deg,
+        width=0.2 * u.deg,
+        psi=120 * u.deg,
+    )
+
+    disp.overlay_moments(hillas, color="w")
+    fig.savefig(tmp_path / "result.png")
+
+
 @pytest.mark.parametrize("pix_type", PixelShape.__members__.values())
 def test_pixel_shapes(pix_type, prod5_lst_cam, tmp_path):
     """test CameraDisplay functionality"""

ctapipe/visualization/utils.py

@@ -0,0 +1,84 @@
+import astropy.units as u
+import numpy as np
+from astropy.coordinates import Angle
+
+from ..containers import CameraHillasParametersContainer, HillasParametersContainer
+
+
+def build_hillas_overlay(hillas, unit, with_label=True, n_sigma=1):
+    """
+    Get position, rotation and text for the hillas parameters label
+    """
+    try:
+        length = hillas.length.to_value(unit)
+        width = hillas.width.to_value(unit)
+    except u.UnitsError:
+        raise ValueError("hillas must be in same frame as geometry")
+
+    # strip off any units
+    if isinstance(hillas, HillasParametersContainer):
+        cog_x = hillas.fov_lon.to_value(unit)
+        cog_y = hillas.fov_lat.to_value(unit)
+    elif isinstance(hillas, CameraHillasParametersContainer):
+        cog_x = hillas.x.to_value(unit)
+        cog_y = hillas.y.to_value(unit)
+    else:
+        raise TypeError(
+            "hillas must be a (Camera)HillasParametersContainer" f", got: {hillas} "
+        )
+
+    psi_rad = hillas.psi.to_value(u.rad)
+    psi_deg = Angle(hillas.psi).wrap_at(180 * u.deg).to_value(u.deg)
+
+    ret = dict(
+        cog_x=cog_x,
+        cog_y=cog_y,
+        width=width,
+        length=length,
+        psi_rad=psi_rad,
+    )
+
+    if not with_label:
+        return ret
+
+    # the following code dealing with x, y, angle
+    # results in the minimal rotation of the text and puts the
+    # label just outside the ellipse
+    if psi_deg < -135:
+        psi_deg += 180
+        psi_rad += np.pi
+    elif psi_deg > 135:
+        psi_deg -= 180
+        psi_rad -= np.pi
+
+    if -45 < psi_deg <= 45:
+        r = 1.2 * n_sigma * width
+        label_x = cog_x + r * np.cos(psi_rad + 0.5 * np.pi)
+        label_y = cog_y + r * np.sin(psi_rad + 0.5 * np.pi)
+        rotation = psi_deg
+    elif 45 < psi_deg <= 135:
+        r = 1.2 * n_sigma * length
+        label_x = cog_x + r * np.cos(psi_rad)
+        label_y = cog_y + r * np.sin(psi_rad)
+        rotation = psi_deg - 90
+    else:
+        r = 1.2 * n_sigma * length
+        label_x = cog_x - r * np.cos(psi_rad)
+        label_y = cog_y - r * np.sin(psi_rad)
+        rotation = psi_deg + 90
+
+    ret["rotation"] = rotation
+    ret["label_x"] = label_x
+    ret["label_y"] = label_y
+
+    if unit == u.deg:
+        sep = ""
+    else:
+        sep = " "
+
+    ret["text"] = (
+        f"({cog_x:.2f}{sep}{unit:unicode}, {cog_y:.2f}{sep}{unit:unicode})\n"
+        f"[w={width:.2f}{sep}{unit:unicode},l={length:.2f}{sep}{unit:unicode}]"
+    )
+
+    return ret

docs/changes/2347.feature.rst

@@ -0,0 +1,3 @@
+Add support for Hillas parameters in ``TelescopeFrame`` to
+``CameraDisplay.overlay_moments`` and make sure that the
+label text does not overlap with the ellipse.