Add EngineeringCameraFrame

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',

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
+ENGINEERING_FRAME_TRANSFORM_MATRIX = np.array([
+    [0, -1, 0],
+    [-1, 0, 0],
+    [0, 0, 1]
+])
+ENGINEERING_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 ENGINEERING_FRAME_TRANSFORM_MATRIX, ENGINEERING_FRAME_TRANSFORM_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

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

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,

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
+    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))

ctapipe/visualization/mpl_camera.py

@@ -94,10 +94,9 @@ def __init__(
             cmap=None,
             allow_pick=False,
             autoupdate=True,
-            autoscale=True
+            autoscale=True,
     ):
         self.axes = ax if ax is not None else plt.gca()
-        self.geom = geometry
         self.pixels = None
         self.colorbar = None
         self.autoupdate = autoupdate
@@ -106,6 +105,8 @@ def __init__(
         self._active_pixel_label = None
         self._axes_overlays = []
 
+        self.geom = geometry
+
         if title is None:
             title = geometry.cam_id
 
@@ -117,25 +118,25 @@ 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.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"):
-                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=self.geom.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=self.geom.pix_rotation.to_value(u.deg),
                     fill=True,
                 )
 

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

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()