From b55b696b05c799c2c16fa87ef59188b2e7eddc52 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 6 Jun 2019 11:07:41 +0200 Subject: [PATCH 01/11] Add EngineeringCameraFrame --- ctapipe/coordinates/__init__.py | 3 +- ctapipe/coordinates/camera_frame.py | 48 +++++++++++++++++++ .../tests/test_engineering_frame.py | 16 +++++++ 3 files changed, 66 insertions(+), 1 deletion(-) create mode 100644 ctapipe/coordinates/tests/test_engineering_frame.py diff --git a/ctapipe/coordinates/__init__.py b/ctapipe/coordinates/__init__.py index 00f93920121..8d7bf9f744a 100644 --- a/ctapipe/coordinates/__init__.py +++ b/ctapipe/coordinates/__init__.py @@ -11,12 +11,13 @@ from .telescope_frame import TelescopeFrame from .nominal_frame import NominalFrame from .ground_frames import GroundFrame, TiltedGroundFrame, project_to_ground -from .camera_frame import CameraFrame +from .camera_frame import CameraFrame, EngineeringCameraFrame __all__ = [ 'TelescopeFrame', 'CameraFrame', + 'EngineeringCameraFrame', 'NominalFrame', 'GroundFrame', 'TiltedGroundFrame', diff --git a/ctapipe/coordinates/camera_frame.py b/ctapipe/coordinates/camera_frame.py index cff4b3ee63f..49c7a901814 100644 --- a/ctapipe/coordinates/camera_frame.py +++ b/ctapipe/coordinates/camera_frame.py @@ -11,12 +11,22 @@ CartesianRepresentation, UnitSphericalRepresentation, AltAz, + AffineTransform, ) from .telescope_frame import TelescopeFrame from .representation import PlanarRepresentation +# Go from SimTel / HESS to MAGIC/FACT/Engineering frame and back +CAMERA_FRAME_TRANSFORM_MATRIX = np.array([ + [0, -1, 0], + [-1, 0, 0], + [0, 0, 1] +]) +CAMERA_FRAME_TRANSFORM_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m) + + class CameraFrame(BaseCoordinateFrame): ''' Camera coordinate frame. @@ -58,6 +68,34 @@ class CameraFrame(BaseCoordinateFrame): location = EarthLocationAttribute(default=None) +class EngineeringCameraFrame(CameraFrame): + ''' + Engineering camera coordinate frame. + + The camera frame is a 2d cartesian frame, + describing position of objects in the focal plane of the telescope. + + The frame is defined as in MAGIC and FACT. + Standing in the dish, looking onto the camera, x points right and y points up. + + HESS, ctapipe and sim_telarray use a different camera coordinate system: + To transform H.E.S.S./ctapipe -> FACT/MAGIC, do x' = -y, y' = -x. + + Attributes + ---------- + focal_length : u.Quantity[length] + Focal length of the telescope as a unit quantity (usually meters) + rotation : u.Quantity[angle] + Rotation angle of the camera (0 deg in most cases) + telescope_pointing : SkyCoord[AltAz] + Pointing direction of the telescope as SkyCoord in AltAz + obstime : Time + Observation time + location : EarthLocation + location of the telescope + ''' + + @frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame) def camera_to_telescope(camera_coord, telescope_frame): ''' @@ -132,3 +170,13 @@ def telescope_to_camera(telescope_coord, camera_frame): ) return camera_frame.realize_frame(representation) + + +@frame_transform_graph.transform(AffineTransform, CameraFrame, EngineeringCameraFrame) +def camera_to_engineering(from_coord, to_frame): + return CAMERA_FRAME_TRANSFORM_MATRIX, CAMERA_FRAME_TRANSFORM_OFFSET + + +@frame_transform_graph.transform(AffineTransform, EngineeringCameraFrame, CameraFrame) +def engineering_to_camera(from_coord, to_frame): + return CAMERA_FRAME_TRANSFORM_MATRIX, CAMERA_FRAME_TRANSFORM_OFFSET diff --git a/ctapipe/coordinates/tests/test_engineering_frame.py b/ctapipe/coordinates/tests/test_engineering_frame.py new file mode 100644 index 00000000000..5db40a0485a --- /dev/null +++ b/ctapipe/coordinates/tests/test_engineering_frame.py @@ -0,0 +1,16 @@ +from astropy.coordinates import SkyCoord +import astropy.units as u + + +def test_conversion(): + from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame + + coord = SkyCoord(x=1 * u.m, y=2 * u.m, frame=CameraFrame()) + eng_coord = coord.transform_to(EngineeringCameraFrame()) + + assert eng_coord.x == -coord.y + assert eng_coord.y == -coord.x + + back = eng_coord.transform_to(CameraFrame()) + assert back.x == coord.x + assert back.y == coord.y From 902e15442c14333899812c2e11aacb8c1d77ee3d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 6 Jun 2019 12:33:30 +0200 Subject: [PATCH 02/11] Allow plotting in engineering frame --- ctapipe/visualization/mpl_camera.py | 41 ++++++++++++++++--------- ctapipe/visualization/tests/test_mpl.py | 15 +++++++++ 2 files changed, 42 insertions(+), 14 deletions(-) diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index e07f27a6a5e..18f042c769b 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/ctapipe/visualization/mpl_camera.py @@ -7,12 +7,15 @@ import numpy as np from astropy import units as u +from astropy.coordinates import SkyCoord from matplotlib import pyplot as plt from matplotlib.collections import PatchCollection from matplotlib.colors import Normalize, LogNorm, SymLogNorm from matplotlib.patches import Ellipse, RegularPolygon, Rectangle from numpy import sqrt +from ..coordinates import CameraFrame, EngineeringCameraFrame + __all__ = ['CameraDisplay'] logger = logging.getLogger(__name__) @@ -56,6 +59,8 @@ class CameraDisplay: rescale the vmin/vmax values when the image changes. This is set to False if `set_limits_*` is called to explicity set data limits. + engineering_frame : bool (default False) + Transform camera coordinates into the engineering frame Notes ----- @@ -94,7 +99,8 @@ def __init__( cmap=None, allow_pick=False, autoupdate=True, - autoscale=True + autoscale=True, + engineering_frame=False, ): self.axes = ax if ax is not None else plt.gca() self.geom = geometry @@ -117,25 +123,32 @@ def __init__( if not hasattr(self.geom, "mask"): self.geom.mask = np.ones_like(self.geom.pix_x.value, dtype=bool) - for xx, yy, aa in zip( - u.Quantity(self.geom.pix_x[self.geom.mask]).value, - u.Quantity(self.geom.pix_y[self.geom.mask]).value, - u.Quantity(np.array(self.geom.pix_area)[self.geom.mask]).value): - + pix_x = self.geom.pix_x[self.geom.mask] + pix_y = self.geom.pix_y[self.geom.mask] + pix_area = self.geom.pix_area[self.geom.mask] + pix_rotation = self.geom.pix_rotation + if engineering_frame: + from_frame = CameraFrame() + coord = SkyCoord(x=pix_x, y=pix_y, frame=from_frame) + trans = coord.transform_to(EngineeringCameraFrame()) + pix_x, pix_y = trans.x, trans.y + pix_rotation = 90 * u.deg - pix_rotation + + for x, y, area in zip(pix_x.value, pix_y.value, pix_area.value): if self.geom.pix_type.startswith("hex"): - rr = sqrt(aa * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON + r = sqrt(area * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON poly = RegularPolygon( - (xx, yy), 6, radius=rr, - orientation=self.geom.pix_rotation.rad, + (x, y), 6, radius=r, + orientation=pix_rotation.to_value(u.rad), fill=True, ) else: - rr = sqrt(aa) + PIXEL_EPSILON + r = sqrt(area) + PIXEL_EPSILON poly = Rectangle( - (xx - rr / 2., yy - rr / 2.), - width=rr, - height=rr, - angle=self.geom.pix_rotation.deg, + (x - r / 2, y - r / 2), + width=r, + height=r, + angle=pix_rotation.to_value(u.deg), fill=True, ) diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 1267e4faeda..734090f98c7 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -53,6 +53,21 @@ def test_camera_display_multiple(): d2.image = image +def test_camera_engineering_frame(): + """ create a figure with 2 subplots, each with a CameraDisplay """ + from ..mpl_camera import CameraDisplay + + geom = CameraGeometry.from_name("LSTCam") + fig, ax = plt.subplots(1, 1) + + CameraDisplay(geom, ax=ax, engineering_frame=True) + + geom = CameraGeometry.from_name("CHEC") + fig, ax = plt.subplots(1, 1) + + CameraDisplay(geom, ax=ax, engineering_frame=True) + + def test_array_display(): from ctapipe.visualization.mpl_array import ArrayDisplay from ctapipe.image.timing_parameters import timing_parameters From dd597683af94af88c101d6bf123171f030f8954f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 6 Jun 2019 14:38:22 +0200 Subject: [PATCH 03/11] Add transform_to to CameraGeometry and use it in plotting --- ctapipe/instrument/camera.py | 45 +++++++++++++++++++++++-- ctapipe/visualization/mpl_camera.py | 35 +++++++++---------- ctapipe/visualization/tests/test_mpl.py | 5 +-- 3 files changed, 61 insertions(+), 24 deletions(-) diff --git a/ctapipe/instrument/camera.py b/ctapipe/instrument/camera.py index c845b936ce0..50c6c9ee4fd 100644 --- a/ctapipe/instrument/camera.py +++ b/ctapipe/instrument/camera.py @@ -6,7 +6,7 @@ import numpy as np from astropy import units as u -from astropy.coordinates import Angle +from astropy.coordinates import Angle, SkyCoord from astropy.table import Table from astropy.utils import lazyproperty from scipy.spatial import cKDTree as KDTree @@ -15,6 +15,7 @@ from ctapipe.utils import get_table_dataset, find_all_matching_datasets from ctapipe.utils.linalg import rotation_matrix_2d +from ctapipe.coordinates import CameraFrame __all__ = ['CameraGeometry'] @@ -68,7 +69,7 @@ class CameraGeometry: def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type, pix_rotation="0d", cam_rotation="0d", - neighbors=None, apply_derotation=True): + neighbors=None, apply_derotation=True, frame=None): if pix_x.ndim != 1 or pix_y.ndim != 1: raise ValueError(f'Pixel coordinates must be 1 dimensional, got {pix_x.ndim}') @@ -84,6 +85,7 @@ def __init__(self, cam_id, pix_id, pix_x, pix_y, pix_area, pix_type, self.pix_rotation = Angle(pix_rotation) self.cam_rotation = Angle(cam_rotation) self._neighbors = neighbors + self.frame = frame if neighbors is not None: if isinstance(neighbors, list): @@ -125,6 +127,45 @@ def __eq__(self, other): (self.pix_y == other.pix_y).all(), ]) + def transform_to(self, frame): + ''' + Transform the pixel coordinates stored in this geometry + and the pixel and camera rotations to another camera coordinate frame. + + Parameters + ---------- + frame: ctapipe.coordinates.CameraFrame + The coordinate frame to transform to. + ''' + if self.frame is None: + self.frame = CameraFrame() + + coord = SkyCoord(x=self.pix_x, y=self.pix_y, frame=self.frame) + trans = coord.transform_to(frame) + + # also transform the unit vectors, to get rotation / mirroring + uv = SkyCoord(x=[1, 0], y=[0, 1], unit=u.m, frame=self.frame) + uv_trans = uv.transform_to(frame) + rot = np.arctan2(uv_trans[0].y, uv_trans[1].y) + det = np.linalg.det([uv_trans.x.value, uv_trans.y.value]) + + cam_rotation = rot + det * self.cam_rotation + pix_rotation = rot + det * self.pix_rotation + + return CameraGeometry( + cam_id=self.cam_id, + pix_id=self.pix_id, + pix_x=trans.x, + pix_y=trans.y, + pix_area=self.pix_area, + pix_type=self.pix_type, + pix_rotation=pix_rotation, + cam_rotation=cam_rotation, + neighbors=None, + apply_derotation=False, + frame=frame, + ) + def __hash__(self): return hash(( self.cam_id, diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index 18f042c769b..4b623e4a535 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/ctapipe/visualization/mpl_camera.py @@ -14,8 +14,6 @@ from matplotlib.patches import Ellipse, RegularPolygon, Rectangle from numpy import sqrt -from ..coordinates import CameraFrame, EngineeringCameraFrame - __all__ = ['CameraDisplay'] logger = logging.getLogger(__name__) @@ -59,8 +57,9 @@ class CameraDisplay: rescale the vmin/vmax values when the image changes. This is set to False if `set_limits_*` is called to explicity set data limits. - engineering_frame : bool (default False) - Transform camera coordinates into the engineering frame + frame : A CameraFrame instance (default None) + If given, transform pixel coordinates into this frame + before plotting the camera. Notes ----- @@ -100,10 +99,9 @@ def __init__( allow_pick=False, autoupdate=True, autoscale=True, - engineering_frame=False, + frame=None, ): self.axes = ax if ax is not None else plt.gca() - self.geom = geometry self.pixels = None self.colorbar = None self.autoupdate = autoupdate @@ -112,6 +110,10 @@ def __init__( self._active_pixel_label = None self._axes_overlays = [] + self.geom = geometry + if frame is not None: + self.geom = self.geom.transform_to(frame) + if title is None: title = geometry.cam_id @@ -123,23 +125,16 @@ def __init__( if not hasattr(self.geom, "mask"): self.geom.mask = np.ones_like(self.geom.pix_x.value, dtype=bool) - pix_x = self.geom.pix_x[self.geom.mask] - pix_y = self.geom.pix_y[self.geom.mask] - pix_area = self.geom.pix_area[self.geom.mask] - pix_rotation = self.geom.pix_rotation - if engineering_frame: - from_frame = CameraFrame() - coord = SkyCoord(x=pix_x, y=pix_y, frame=from_frame) - trans = coord.transform_to(EngineeringCameraFrame()) - pix_x, pix_y = trans.x, trans.y - pix_rotation = 90 * u.deg - pix_rotation - - for x, y, area in zip(pix_x.value, pix_y.value, pix_area.value): + pix_x = self.geom.pix_x.value[self.geom.mask] + pix_y = self.geom.pix_y.value[self.geom.mask] + pix_area = self.geom.pix_area.value[self.geom.mask] + + for x, y, area in zip(pix_x, pix_y, pix_area): if self.geom.pix_type.startswith("hex"): r = sqrt(area * 2 / 3 / sqrt(3)) + 2 * PIXEL_EPSILON poly = RegularPolygon( (x, y), 6, radius=r, - orientation=pix_rotation.to_value(u.rad), + orientation=self.geom.pix_rotation.to_value(u.rad), fill=True, ) else: @@ -148,7 +143,7 @@ def __init__( (x - r / 2, y - r / 2), width=r, height=r, - angle=pix_rotation.to_value(u.deg), + angle=self.geom.pix_rotation.to_value(u.deg), fill=True, ) diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 734090f98c7..d5018ac5912 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -6,6 +6,7 @@ from ctapipe.instrument import (CameraGeometry, SubarrayDescription, TelescopeDescription) from ctapipe.io.containers import HillasParametersContainer +from ctapipe.coordinates import EngineeringCameraFrame from numpy import ones from astropy import units as u @@ -60,12 +61,12 @@ def test_camera_engineering_frame(): geom = CameraGeometry.from_name("LSTCam") fig, ax = plt.subplots(1, 1) - CameraDisplay(geom, ax=ax, engineering_frame=True) + CameraDisplay(geom, ax=ax, frame=EngineeringCameraFrame()) geom = CameraGeometry.from_name("CHEC") fig, ax = plt.subplots(1, 1) - CameraDisplay(geom, ax=ax, engineering_frame=True) + CameraDisplay(geom, ax=ax, frame=EngineeringCameraFrame()) def test_array_display(): From 02330ebb6b3bb0eb8a3ff073b6a5d8ad9bc2bb0e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 6 Jun 2019 14:42:10 +0200 Subject: [PATCH 04/11] Rename constants --- ctapipe/coordinates/camera_frame.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/ctapipe/coordinates/camera_frame.py b/ctapipe/coordinates/camera_frame.py index 49c7a901814..c4eab94d523 100644 --- a/ctapipe/coordinates/camera_frame.py +++ b/ctapipe/coordinates/camera_frame.py @@ -19,12 +19,12 @@ # Go from SimTel / HESS to MAGIC/FACT/Engineering frame and back -CAMERA_FRAME_TRANSFORM_MATRIX = np.array([ +ENGINEERING_FRAME_TRANSFORM_MATRIX = np.array([ [0, -1, 0], [-1, 0, 0], [0, 0, 1] ]) -CAMERA_FRAME_TRANSFORM_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m) +ENGINEERING_FRAME_TRANSFORM_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m) class CameraFrame(BaseCoordinateFrame): @@ -174,9 +174,9 @@ def telescope_to_camera(telescope_coord, camera_frame): @frame_transform_graph.transform(AffineTransform, CameraFrame, EngineeringCameraFrame) def camera_to_engineering(from_coord, to_frame): - return CAMERA_FRAME_TRANSFORM_MATRIX, CAMERA_FRAME_TRANSFORM_OFFSET + return ENGINEERING_FRAME_TRANSFORM_MATRIX, ENGINEERING_FRAME_TRANSFORM_OFFSET @frame_transform_graph.transform(AffineTransform, EngineeringCameraFrame, CameraFrame) def engineering_to_camera(from_coord, to_frame): - return CAMERA_FRAME_TRANSFORM_MATRIX, CAMERA_FRAME_TRANSFORM_OFFSET + return ENGINEERING_FRAME_TRANSFORM_MATRIX, ENGINEERING_FRAME_TRANSFORM_OFFSET From 2b407575fd42af1f3adb1d91682afa6fd3b9be4d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Thu, 6 Jun 2019 16:53:41 +0200 Subject: [PATCH 05/11] Move transformation from plotting to Geometry --- ctapipe/instrument/tests/test_camera.py | 14 +++++++++++++- ctapipe/visualization/mpl_camera.py | 6 ------ ctapipe/visualization/tests/test_mpl.py | 15 --------------- 3 files changed, 13 insertions(+), 22 deletions(-) diff --git a/ctapipe/instrument/tests/test_camera.py b/ctapipe/instrument/tests/test_camera.py index 29313ffeb11..e830c88e918 100644 --- a/ctapipe/instrument/tests/test_camera.py +++ b/ctapipe/instrument/tests/test_camera.py @@ -291,4 +291,16 @@ def test_hashing(): def test_camera_from_name(camera_name): """ check we can construct all cameras from name""" camera = CameraGeometry.from_name(camera_name) - assert str(camera) == camera_name \ No newline at end of file + assert str(camera) == camera_name + + +@pytest.mark.parametrize("camera_name", CameraGeometry.get_known_camera_names()) +def test_camera_engineering_frame(camera_name): + from ctapipe.coordinates import EngineeringCameraFrame + + geom = CameraGeometry.from_name(camera_name) + trans_geom = geom.transform_to(EngineeringCameraFrame()) + + unit = geom.pix_x.unit + assert np.allclose(geom.pix_x.to_value(unit), -trans_geom.pix_y.to_value(unit)) + assert np.allclose(geom.pix_y.to_value(unit), -trans_geom.pix_x.to_value(unit)) diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index 4b623e4a535..5087630b768 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/ctapipe/visualization/mpl_camera.py @@ -57,9 +57,6 @@ class CameraDisplay: rescale the vmin/vmax values when the image changes. This is set to False if `set_limits_*` is called to explicity set data limits. - frame : A CameraFrame instance (default None) - If given, transform pixel coordinates into this frame - before plotting the camera. Notes ----- @@ -99,7 +96,6 @@ def __init__( allow_pick=False, autoupdate=True, autoscale=True, - frame=None, ): self.axes = ax if ax is not None else plt.gca() self.pixels = None @@ -111,8 +107,6 @@ def __init__( self._axes_overlays = [] self.geom = geometry - if frame is not None: - self.geom = self.geom.transform_to(frame) if title is None: title = geometry.cam_id diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index d5018ac5912..88fa430d114 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -54,21 +54,6 @@ def test_camera_display_multiple(): d2.image = image -def test_camera_engineering_frame(): - """ create a figure with 2 subplots, each with a CameraDisplay """ - from ..mpl_camera import CameraDisplay - - geom = CameraGeometry.from_name("LSTCam") - fig, ax = plt.subplots(1, 1) - - CameraDisplay(geom, ax=ax, frame=EngineeringCameraFrame()) - - geom = CameraGeometry.from_name("CHEC") - fig, ax = plt.subplots(1, 1) - - CameraDisplay(geom, ax=ax, frame=EngineeringCameraFrame()) - - def test_array_display(): from ctapipe.visualization.mpl_array import ArrayDisplay from ctapipe.image.timing_parameters import timing_parameters From 6d4deb7c8d5847b00138887f9107e9f5695cc578 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Fri, 7 Jun 2019 01:10:26 +0200 Subject: [PATCH 06/11] Add example to the docs --- docs/ctapipe_api/coordinates/index.rst | 8 ++++++++ examples/plot_camera_frames.py | 25 +++++++++++++++++++++++++ 2 files changed, 33 insertions(+) create mode 100644 examples/plot_camera_frames.py diff --git a/docs/ctapipe_api/coordinates/index.rst b/docs/ctapipe_api/coordinates/index.rst index 159a584176b..6bd31b4a368 100644 --- a/docs/ctapipe_api/coordinates/index.rst +++ b/docs/ctapipe_api/coordinates/index.rst @@ -27,6 +27,7 @@ Coordinates should be described using an `astropy.coordinates.SkyCoord` in the a The following special frames are defined for CTA: * `CameraFrame` +* `EngineeringCameraFrame` * `TelescopeFrame` * `NominalFrame` * `GroundFrame` @@ -35,6 +36,13 @@ The following special frames are defined for CTA: they can be transformed to and from any other `astropy.coordinates` frame, like `astropy.coordinates.AltAz` or `astropy.coordinates.ICRS` (RA/Dec) +The two different coordinate frames are shown here: + +.. plot:: ../examples/plot_camera_frames.py + :include-source: + +The `CameraFrame` is used internally and is used by H.E.S.S. and `sim_telarray`, +the `EngineeringCameraFrame` is used by `MAGIC` and `FACT`. Reference/API diff --git a/examples/plot_camera_frames.py b/examples/plot_camera_frames.py new file mode 100644 index 00000000000..d6717d1c98a --- /dev/null +++ b/examples/plot_camera_frames.py @@ -0,0 +1,25 @@ +from ctapipe.instrument import CameraGeometry +from ctapipe.visualization import CameraDisplay +import matplotlib.pyplot as plt +from ctapipe.coordinates import EngineeringCameraFrame +from ctapipe.image.toymodel import Gaussian +import astropy.units as u + + +fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(6, 3)) + +model = Gaussian(0 * u.m, 0.1 * u.m, 0.3 * u.m, 0.05 * u.m, 25 * u.deg) +cam = CameraGeometry.from_name('FlashCam') +image, *_ = model.generate_image(cam, 2500) + +CameraDisplay(cam, ax=axs[0], image=image) +CameraDisplay( + cam.transform_to(EngineeringCameraFrame()), + ax=axs[1], + image=image, +) + +axs[0].set_title('CameraFrame') +axs[1].set_title('EngineeringCameraFrame') + +plt.show() From 1b8ca1d04f6d3000576cc8447e1b2e61c73b6cc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Fri, 7 Jun 2019 01:27:03 +0200 Subject: [PATCH 07/11] Remove unused imports --- ctapipe/visualization/mpl_camera.py | 1 - ctapipe/visualization/tests/test_mpl.py | 1 - 2 files changed, 2 deletions(-) diff --git a/ctapipe/visualization/mpl_camera.py b/ctapipe/visualization/mpl_camera.py index 5087630b768..9af7ae2ad2a 100644 --- a/ctapipe/visualization/mpl_camera.py +++ b/ctapipe/visualization/mpl_camera.py @@ -7,7 +7,6 @@ import numpy as np from astropy import units as u -from astropy.coordinates import SkyCoord from matplotlib import pyplot as plt from matplotlib.collections import PatchCollection from matplotlib.colors import Normalize, LogNorm, SymLogNorm diff --git a/ctapipe/visualization/tests/test_mpl.py b/ctapipe/visualization/tests/test_mpl.py index 88fa430d114..1267e4faeda 100644 --- a/ctapipe/visualization/tests/test_mpl.py +++ b/ctapipe/visualization/tests/test_mpl.py @@ -6,7 +6,6 @@ from ctapipe.instrument import (CameraGeometry, SubarrayDescription, TelescopeDescription) from ctapipe.io.containers import HillasParametersContainer -from ctapipe.coordinates import EngineeringCameraFrame from numpy import ones from astropy import units as u From e98b2ef0ea0f242c686681e1274b8d4d496e2597 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Mon, 10 Jun 2019 17:49:38 +0200 Subject: [PATCH 08/11] Add dual-mirror case --- ctapipe/coordinates/camera_frame.py | 30 ++++++++++++++++--- .../tests/test_engineering_frame.py | 25 +++++++++++----- 2 files changed, 44 insertions(+), 11 deletions(-) diff --git a/ctapipe/coordinates/camera_frame.py b/ctapipe/coordinates/camera_frame.py index c4eab94d523..3399a5ff6a4 100644 --- a/ctapipe/coordinates/camera_frame.py +++ b/ctapipe/coordinates/camera_frame.py @@ -4,6 +4,7 @@ BaseCoordinateFrame, CoordinateAttribute, QuantityAttribute, + Attribute, TimeAttribute, EarthLocationAttribute, FunctionTransform, @@ -18,13 +19,28 @@ from .representation import PlanarRepresentation +class MirrorAttribute(Attribute): + def convert_input(self, value): + if value in (1, 2): + return value, False + + raise ValueError('Only 1 or 2 mirrors supported') + + # Go from SimTel / HESS to MAGIC/FACT/Engineering frame and back -ENGINEERING_FRAME_TRANSFORM_MATRIX = np.array([ +CAMERA_TO_ENGINEERING_1M_MATRIX = np.array([ [0, -1, 0], [-1, 0, 0], [0, 0, 1] ]) -ENGINEERING_FRAME_TRANSFORM_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m) +ENGINEERING_1M_TO_CAMERA_MATRIX = CAMERA_TO_ENGINEERING_1M_MATRIX +CAMERA_TO_ENGINEERING_2M_MATRIX = np.array([ + [0, 1, 0], + [-1, 0, 0], + [0, 0, 1] +]) +ENGINEERING_2M_TO_CAMERA_MATRIX = CAMERA_TO_ENGINEERING_2M_MATRIX.T +ZERO_OFFSET = CartesianRepresentation(0, 0, 0, unit=u.m) class CameraFrame(BaseCoordinateFrame): @@ -94,6 +110,7 @@ class EngineeringCameraFrame(CameraFrame): location : EarthLocation location of the telescope ''' + n_mirrors = MirrorAttribute(default=1) @frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame) @@ -174,9 +191,14 @@ def telescope_to_camera(telescope_coord, camera_frame): @frame_transform_graph.transform(AffineTransform, CameraFrame, EngineeringCameraFrame) def camera_to_engineering(from_coord, to_frame): - return ENGINEERING_FRAME_TRANSFORM_MATRIX, ENGINEERING_FRAME_TRANSFORM_OFFSET + if to_frame.n_mirrors == 1: + return CAMERA_TO_ENGINEERING_1M_MATRIX, ZERO_OFFSET + + return CAMERA_TO_ENGINEERING_2M_MATRIX, ZERO_OFFSET @frame_transform_graph.transform(AffineTransform, EngineeringCameraFrame, CameraFrame) def engineering_to_camera(from_coord, to_frame): - return ENGINEERING_FRAME_TRANSFORM_MATRIX, ENGINEERING_FRAME_TRANSFORM_OFFSET + if from_coord.n_mirrors == 1: + return ENGINEERING_1M_TO_CAMERA_MATRIX, ZERO_OFFSET + return ENGINEERING_2M_TO_CAMERA_MATRIX, ZERO_OFFSET diff --git a/ctapipe/coordinates/tests/test_engineering_frame.py b/ctapipe/coordinates/tests/test_engineering_frame.py index 5db40a0485a..23b9be01c0b 100644 --- a/ctapipe/coordinates/tests/test_engineering_frame.py +++ b/ctapipe/coordinates/tests/test_engineering_frame.py @@ -5,12 +5,23 @@ def test_conversion(): from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame - coord = SkyCoord(x=1 * u.m, y=2 * u.m, frame=CameraFrame()) - eng_coord = coord.transform_to(EngineeringCameraFrame()) + coords = SkyCoord(x=[3, 1] * u.m, y=[2, 4] * u.m, frame=CameraFrame()) - assert eng_coord.x == -coord.y - assert eng_coord.y == -coord.x + for coord in coords: + eng_coord = coord.transform_to(EngineeringCameraFrame()) - back = eng_coord.transform_to(CameraFrame()) - assert back.x == coord.x - assert back.y == coord.y + assert eng_coord.x == -coord.y + assert eng_coord.y == -coord.x + + back = eng_coord.transform_to(CameraFrame()) + assert back.x == coord.x + assert back.y == coord.y + + eng_coord = coord.transform_to(EngineeringCameraFrame(n_mirrors=2)) + + assert eng_coord.x == coord.y + assert eng_coord.y == -coord.x + + back = eng_coord.transform_to(CameraFrame()) + assert back.x == coord.x + assert back.y == coord.y From e461683fae4021beab3a5ce66c338bfb0ddbb8db Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Wed, 12 Jun 2019 11:37:54 +0200 Subject: [PATCH 09/11] Style fixes --- .../tests/test_engineering_frame.py | 6 ++++ ctapipe/instrument/tests/test_camera.py | 4 ++- examples/plot_camera_frames.py | 35 ++++++++++++------- 3 files changed, 31 insertions(+), 14 deletions(-) diff --git a/ctapipe/coordinates/tests/test_engineering_frame.py b/ctapipe/coordinates/tests/test_engineering_frame.py index 23b9be01c0b..4d454d74822 100644 --- a/ctapipe/coordinates/tests/test_engineering_frame.py +++ b/ctapipe/coordinates/tests/test_engineering_frame.py @@ -1,8 +1,14 @@ +''' +Tests for the conversion between camera coordinate frames +''' from astropy.coordinates import SkyCoord import astropy.units as u def test_conversion(): + ''' + Test conversion between CameraFrame and EngineeringCameraFrame + ''' from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame coords = SkyCoord(x=[3, 1] * u.m, y=[2, 4] * u.m, frame=CameraFrame()) diff --git a/ctapipe/instrument/tests/test_camera.py b/ctapipe/instrument/tests/test_camera.py index e830c88e918..81220613467 100644 --- a/ctapipe/instrument/tests/test_camera.py +++ b/ctapipe/instrument/tests/test_camera.py @@ -287,6 +287,7 @@ def test_hashing(): assert len(set([cam1, cam2, cam3])) == 2 + @pytest.mark.parametrize("camera_name", CameraGeometry.get_known_camera_names()) def test_camera_from_name(camera_name): """ check we can construct all cameras from name""" @@ -295,7 +296,8 @@ def test_camera_from_name(camera_name): @pytest.mark.parametrize("camera_name", CameraGeometry.get_known_camera_names()) -def test_camera_engineering_frame(camera_name): +def test_camera_coordinate_transform(camera_name): + '''test conversion of the coordinates stored in a camera frame''' from ctapipe.coordinates import EngineeringCameraFrame geom = CameraGeometry.from_name(camera_name) diff --git a/examples/plot_camera_frames.py b/examples/plot_camera_frames.py index d6717d1c98a..a36ef581051 100644 --- a/examples/plot_camera_frames.py +++ b/examples/plot_camera_frames.py @@ -1,3 +1,7 @@ +''' +Plot the same event in two camera displays showing the +different coordinate frames for camera coordinates. +''' from ctapipe.instrument import CameraGeometry from ctapipe.visualization import CameraDisplay import matplotlib.pyplot as plt @@ -6,20 +10,25 @@ import astropy.units as u -fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(6, 3)) +def main(): + fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(6, 3)) -model = Gaussian(0 * u.m, 0.1 * u.m, 0.3 * u.m, 0.05 * u.m, 25 * u.deg) -cam = CameraGeometry.from_name('FlashCam') -image, *_ = model.generate_image(cam, 2500) + model = Gaussian(0 * u.m, 0.1 * u.m, 0.3 * u.m, 0.05 * u.m, 25 * u.deg) + cam = CameraGeometry.from_name('FlashCam') + image, *_ = model.generate_image(cam, 2500) -CameraDisplay(cam, ax=axs[0], image=image) -CameraDisplay( - cam.transform_to(EngineeringCameraFrame()), - ax=axs[1], - image=image, -) + CameraDisplay(cam, ax=axs[0], image=image) + CameraDisplay( + cam.transform_to(EngineeringCameraFrame()), + ax=axs[1], + image=image, + ) -axs[0].set_title('CameraFrame') -axs[1].set_title('EngineeringCameraFrame') + axs[0].set_title('CameraFrame') + axs[1].set_title('EngineeringCameraFrame') -plt.show() + plt.show() + + +if __name__ == '__main__': + main() From edfe26f169cfafaa38bafa0cd83614f3977307fe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Wed, 12 Jun 2019 11:40:07 +0200 Subject: [PATCH 10/11] Update docs --- docs/ctapipe_api/coordinates/index.rst | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/docs/ctapipe_api/coordinates/index.rst b/docs/ctapipe_api/coordinates/index.rst index 6bd31b4a368..3d4f14220c4 100644 --- a/docs/ctapipe_api/coordinates/index.rst +++ b/docs/ctapipe_api/coordinates/index.rst @@ -41,8 +41,10 @@ The two different coordinate frames are shown here: .. plot:: ../examples/plot_camera_frames.py :include-source: -The `CameraFrame` is used internally and is used by H.E.S.S. and `sim_telarray`, -the `EngineeringCameraFrame` is used by `MAGIC` and `FACT`. +The `CameraFrame` is used internally in `ctapipe` and comes from `sim_telarray`. +It abstracts the transformation differences between 1 and 2 mirror telescopes away. +The `EngineeringCameraFrame` is used by `MAGIC`, `FACT` and the `H.E.S.S.` analysis +software. Reference/API From 0343142ceaa576da0185e7260a58ec33ca2c83d4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Maximilian=20N=C3=B6the?= Date: Wed, 12 Jun 2019 11:41:37 +0200 Subject: [PATCH 11/11] Add missing docstring --- ctapipe/coordinates/camera_frame.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/ctapipe/coordinates/camera_frame.py b/ctapipe/coordinates/camera_frame.py index 3399a5ff6a4..be1b52b645f 100644 --- a/ctapipe/coordinates/camera_frame.py +++ b/ctapipe/coordinates/camera_frame.py @@ -20,7 +20,10 @@ class MirrorAttribute(Attribute): + '''A frame Attribute that can only store the integers 1 and 2''' + def convert_input(self, value): + '''make sure input is 1 or 2''' if value in (1, 2): return value, False