From 3c4309a552448aa7ff4c2546cea479a0f9ea2bcf Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 13 Jun 2023 12:14:25 +0200 Subject: [PATCH 1/5] Add support for hillas parameters in TelescopeFrame to CameraDisplay --- ctapipe/visualization/mpl_camera.py | 30 ++++++++++++++++++++++------- 1 file changed, 23 insertions(+), 7 deletions(-) diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index 51a6b98fc81..46fc078a1f2 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/ctapipe/visualization/mpl_camera.py @@ -12,6 +12,7 @@ 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 @@ -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,16 +481,29 @@ 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 - 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) + 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} " + ) el = self.add_ellipse( centroid=(cen_x, cen_y), - length=length * 2, - width=width * 2, + length=n_sigma * length * 2, + width=n_sigma * width * 2, angle=hillas_parameters.psi.to_value(u.rad), **kwargs, ) From 68b8d955be302300df9d0027664e93102b98c1e1 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 13 Jun 2023 17:01:00 +0200 Subject: [PATCH 2/5] Add test for hillas in telescope frame --- ctapipe/visualization/tests/test_mpl.py | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) 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""" From b299e01c751ef7d7bac3a8fcb8a8ae5c1a7db639 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 13 Jun 2023 17:01:33 +0200 Subject: [PATCH 3/5] Make sure overlay label is not clashing with ellipse --- ctapipe/visualization/mpl_camera.py | 47 ++++++++++++++++++++++++----- 1 file changed, 40 insertions(+), 7 deletions(-) diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index 46fc078a1f2..68f4a7cf64c 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/ctapipe/visualization/mpl_camera.py @@ -7,6 +7,7 @@ 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 @@ -500,27 +501,59 @@ def overlay_moments( f", got: {hillas_parameters} " ) + 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=hillas_parameters.psi.to_value(u.rad), + angle=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( - cen_x, - cen_y, + label_x, + label_y, "({:.02f},{:.02f})\n[w={:.02f},l={:.02f}]".format( - hillas_parameters.x, - hillas_parameters.y, - hillas_parameters.width, - hillas_parameters.length, + cen_x * self.unit, + cen_y * self.unit, + hillas_parameters.width.to(self.unit), + hillas_parameters.length.to(self.unit), ), color=el.get_edgecolor(), + va="bottom", + ha="center", + rotation=rotation, + rotation_mode="anchor", ) self._axes_overlays.append(text) From 6c38cb2ad51c2790f6441ef2a915e39289de5e34 Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Tue, 13 Jun 2023 17:09:25 +0200 Subject: [PATCH 4/5] add changelog --- docs/changes/2347.feature.rst | 3 +++ 1 file changed, 3 insertions(+) create mode 100644 docs/changes/2347.feature.rst 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. From 0bc6b46206d34906dcf57e414cf51193244e8ffe Mon Sep 17 00:00:00 2001 From: Maximilian Linhoff Date: Wed, 14 Jun 2023 14:32:56 +0200 Subject: [PATCH 5/5] Add support for TelescopeFrame hillas parameters to Bokeh camera display --- ctapipe/visualization/bokeh.py | 87 +++++++++++++++-------------- ctapipe/visualization/mpl_camera.py | 76 +++++-------------------- ctapipe/visualization/utils.py | 84 ++++++++++++++++++++++++++++ 3 files changed, 144 insertions(+), 103 deletions(-) create mode 100644 ctapipe/visualization/utils.py 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 68f4a7cf64c..739ebb8703d 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/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", ) 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