Add support for TelescopeFrame hillas parameters to Bokeh camera display

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

@@ -7,15 +7,14 @@
 
 import numpy as np
 from astropy import units as u
-from astropy.coordinates import Angle
 from matplotlib import pyplot as plt
 from matplotlib.collections import PatchCollection
 from matplotlib.colors import LogNorm, Normalize, SymLogNorm
 from matplotlib.patches import Circle, Ellipse, RegularPolygon
 
-from ..containers import CameraHillasParametersContainer, HillasParametersContainer
 from ..coordinates import get_representation_component_names
 from ..instrument import PixelShape
+from .utils import build_hillas_overlay
 
 __all__ = ["CameraDisplay"]
 
@@ -482,77 +481,32 @@ def overlay_moments(
         if not keep_old:
             self.clear_overlays()
 
-        try:
-            length = hillas_parameters.length.to_value(self.unit)
-            width = hillas_parameters.width.to_value(self.unit)
-        except u.UnitsError:
-            raise ValueError("hillas_parameters must be in same frame as geometry")
-
-        # strip off any units
-        if isinstance(hillas_parameters, HillasParametersContainer):
-            cen_x = hillas_parameters.fov_lon.to_value(self.unit)
-            cen_y = hillas_parameters.fov_lat.to_value(self.unit)
-        elif isinstance(hillas_parameters, CameraHillasParametersContainer):
-            cen_x = hillas_parameters.x.to_value(self.unit)
-            cen_y = hillas_parameters.y.to_value(self.unit)
-        else:
-            raise TypeError(
-                "hillas_parameters must be a (Camera)HillasParametersContainer"
-                f", got: {hillas_parameters} "
-            )
+        params = build_hillas_overlay(
+            hillas_parameters,
+            self.unit,
+            n_sigma=n_sigma,
+            with_label=with_label,
+        )
 
-        psi_rad = hillas_parameters.psi.to_value(u.rad)
         el = self.add_ellipse(
-            centroid=(cen_x, cen_y),
-            length=n_sigma * length * 2,
-            width=n_sigma * width * 2,
-            angle=psi_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:
-            # the following code dealing with x, y, angle
-            # results in the minimal rotation of the text and puts the
-            # label just outside the ellipse
-            psi_deg = Angle(hillas_parameters.psi).wrap_at(180 * u.deg).to_value(u.deg)
-            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 = cen_x + r * np.cos(psi_rad + 0.5 * np.pi)
-                label_y = cen_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 = cen_x + r * np.cos(psi_rad)
-                label_y = cen_y + r * np.sin(psi_rad)
-                rotation = psi_deg - 90
-            else:
-                r = 1.2 * n_sigma * length
-                label_x = cen_x - r * np.cos(psi_rad)
-                label_y = cen_y - r * np.sin(psi_rad)
-                rotation = psi_deg + 90
-
             text = self.axes.text(
-                label_x,
-                label_y,
-                "({:.02f},{:.02f})\n[w={:.02f},l={:.02f}]".format(
-                    cen_x * self.unit,
-                    cen_y * self.unit,
-                    hillas_parameters.width.to(self.unit),
-                    hillas_parameters.length.to(self.unit),
-                ),
+                params["label_x"],
+                params["label_y"],
+                params["text"],
                 color=el.get_edgecolor(),
                 va="bottom",
                 ha="center",
-                rotation=rotation,
+                rotation=params["rotation"],
                 rotation_mode="anchor",
             )
 

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