diff --git a/ctapipe/visualization/bokeh.py b/ctapipe/visualization/bokeh.py index c99b83aba39..cdfeceebf7b 100644 --- a/ctapipe/visualization/bokeh.py +++ b/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,6 +328,8 @@ 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) @@ -334,8 +337,9 @@ def _init_datasource(self, image=None): 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: diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index 51a6b98fc81..739ebb8703d 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/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,6 +472,8 @@ 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) @@ -478,17 +481,18 @@ def overlay_moments( 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, ) @@ -496,15 +500,14 @@ def overlay_moments( 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) diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index fb8e69cf712..31867d5d158 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/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""" diff --git a/ctapipe/visualization/utils.py b/ctapipe/visualization/utils.py new file mode 100644 index 00000000000..94988439f64 --- /dev/null +++ b/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 diff --git a/docs/changes/2347.feature.rst b/docs/changes/2347.feature.rst new file mode 100644 index 00000000000..7838705f765 --- /dev/null +++ b/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.