Merge pull request #72 from fact-project/coordinates

Coordinates

docs/fact.coordinates.rst

@@ -0,0 +1,26 @@
+fact.coordinates package
+========================
+
+.. automodule:: fact.coordinates
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+Submodules
+----------
+
+fact.coordinates.camera_frame module
+------------------------------------
+
+.. automodule:: fact.coordinates.camera_frame
+    :members:
+    :undoc-members:
+    :show-inheritance:
+
+fact.coordinates.utils module
+-----------------------------
+
+.. automodule:: fact.coordinates.utils
+    :members:
+    :undoc-members:
+    :show-inheritance:

docs/fact.rst

@@ -13,6 +13,7 @@ Subpackages
 
     fact.analysis
     fact.auxservices
+    fact.coordinates
     fact.credentials
     fact.factdb
     fact.instrument
@@ -22,22 +23,23 @@ Subpackages
 Submodules
 ----------
 
-fact\.encrypt\_credentials module
----------------------------------
+fact\.io module
+---------------
 
-.. automodule:: fact.encrypt_credentials
+.. automodule:: fact.io
     :members:
     :undoc-members:
     :show-inheritance:
 
-fact\.io module
----------------
+fact\.encrypt\_credentials module
+---------------------------------
 
-.. automodule:: fact.io
+.. automodule:: fact.encrypt_credentials
     :members:
     :undoc-members:
     :show-inheritance:
 
+
 fact\.qla module
 ----------------
 

fact/VERSION

@@ -1 +1 @@
-0.12.1
+0.13.0

fact/coordinates/__init__.py

@@ -0,0 +1,12 @@
+from .camera_frame import CameraCoordinate
+from .utils import equatorial_to_camera, camera_to_equatorial
+from .utils import horizontal_to_camera, camera_to_horizontal
+
+
+__all__ = [
+    'CameraCoordinate',
+    'equatorial_to_camera',
+    'camera_to_equatorial'
+    'horizontal_to_camera',
+    'camera_to_horizontal'
+]

fact/coordinates/camera_frame.py

@@ -0,0 +1,75 @@
+from astropy.coordinates import (
+    BaseCoordinateFrame,
+    AltAz,
+    frame_transform_graph,
+    FunctionTransform
+)
+
+try:
+    from astropy.coordinates import Attribute
+except ImportError:
+    # for astropy <= 2.0.0
+    from astropy.coordinates import FrameAttribute as Attribute
+
+from astropy.coordinates.matrix_utilities import rotation_matrix
+from astropy.coordinates.representation import CartesianRepresentation
+import astropy.units as u
+
+from .representation import PlanarRepresentation
+from ..instrument.constants import FOCAL_LENGTH_MM
+import numpy as np
+
+focal_length = FOCAL_LENGTH_MM * u.mm
+
+
+class CameraCoordinate(BaseCoordinateFrame):
+    '''Astropy CoordinateFrame representing coordinates in the CameraPlane'''
+    default_representation = PlanarRepresentation
+    pointing_direction = Attribute(default=None)
+
+
+@frame_transform_graph.transform(FunctionTransform, CameraCoordinate, AltAz)
+def camera_to_altaz(camera_frame, altaz):
+    if camera_frame.pointing_direction is None:
+        raise AttributeError('Pointing Direction must be set')
+
+    x = camera_frame.x.copy()
+    y = camera_frame.y.copy()
+
+    z = 1 / np.sqrt(1 + (x / focal_length)**2 + (y / focal_length)**2)
+    x *= z / focal_length
+    y *= z / focal_length
+
+    cartesian = CartesianRepresentation(x, y, z, copy=False)
+
+    rot_z_az = rotation_matrix(-camera_frame.pointing_direction.az, 'z')
+    rot_y_zd = rotation_matrix(-camera_frame.pointing_direction.zen, 'y')
+
+    cartesian = cartesian.transform(rot_y_zd)
+    cartesian = cartesian.transform(rot_z_az)
+
+    altitude = 90 * u.deg - np.arccos(cartesian.z)
+    azimuth = np.arctan2(cartesian.y, cartesian.x)
+
+    return AltAz(
+        alt=altitude, az=azimuth,
+        location=altaz.location,
+        obstime=altaz.obstime,
+    )
+
+
+@frame_transform_graph.transform(FunctionTransform, AltAz, CameraCoordinate)
+def altaz_to_camera(altaz, camera_frame):
+    cartesian = altaz.cartesian
+
+    rot_z_az = rotation_matrix(camera_frame.pointing_direction.az, 'z')
+    rot_y_zd = rotation_matrix(camera_frame.pointing_direction.zen, 'y')
+
+    cartesian = cartesian.transform(rot_z_az)
+    cartesian = cartesian.transform(rot_y_zd)
+
+    return CameraCoordinate(
+        x=cartesian.x * focal_length / cartesian.z,
+        y=cartesian.y * focal_length / cartesian.z,
+        pointing_direction=camera_frame.pointing_direction,
+    )

fact/coordinates/representation.py

@@ -0,0 +1,80 @@
+from astropy.coordinates import BaseRepresentation, CartesianRepresentation
+import astropy.units as u
+from collections import OrderedDict
+from astropy.utils.compat.numpy import broadcast_arrays
+
+
+class PlanarRepresentation(BaseRepresentation):
+    '''
+    Representation of a point in a 2D plane. This is needed for coordinate
+    frames to store their coordinates internally.
+    This is essentially a copy of the Cartesian representation used in astropy.
+    Copied from ctapipe.
+
+    Parameters
+    ----------
+    x, y : `~astropy.units.Quantity`
+        The x and y coordinates of the point(s). If ``x`` and ``y``have
+        different shapes, they should be broadcastable.
+    copy : bool, optional
+        If True arrays will be copied rather than referenced.
+    '''
+    attr_classes = OrderedDict([('x', u.Quantity),
+                                ('y', u.Quantity)])
+
+    def __init__(self, x, y, copy=True, **kwargs):
+
+        if x is None or y is None:
+            raise ValueError(
+                'x and y are required to instantiate CartesianRepresentation'
+            )
+
+        if not isinstance(x, self.attr_classes['x']):
+            raise TypeError('x should be a {0}'.format(self.attr_classes['x'].__name__))
+
+        if not isinstance(y, self.attr_classes['y']):
+            raise TypeError('y should be a {0}'.format(self.attr_classes['y'].__name__))
+
+        x = self.attr_classes['x'](x, copy=copy)
+        y = self.attr_classes['y'](y, copy=copy)
+
+        if not (x.unit.physical_type == y.unit.physical_type):
+            raise u.UnitsError("x and y should have matching physical types")
+
+        try:
+            x, y = broadcast_arrays(x, y, subok=True)
+        except ValueError:
+            raise ValueError("Input parameters x and y cannot be broadcast")
+
+        self._x = x
+        self._y = y
+        self._differentials = {}
+
+    @property
+    def x(self):
+        '''
+        The x component of the point(s).
+        '''
+        return self._x
+
+    @property
+    def y(self):
+        '''
+        The y component of the point(s).
+        '''
+        return self._y
+
+    @property
+    def xy(self):
+        return u.Quantity((self._x, self._y))
+
+    @property
+    def components(self):
+        return 'x', 'y'
+
+    @classmethod
+    def from_cartesian(cls, cartesian):
+        return cls(x=cartesian.x, y=cartesian.y)
+
+    def to_cartesian(self):
+        return CartesianRepresentation(x=self._x, y=self._y, z=0*self._x.unit)