import warnings
from collections import OrderedDict
import astropy.units as u
import numpy as np
from astropy.coordinates import (
AltAz,
Angle,
BaseCoordinateFrame,
BaseRepresentation,
CartesianRepresentation,
CoordinateAttribute,
DynamicMatrixTransform,
EarthLocation,
EarthLocationAttribute,
FunctionTransform,
QuantityAttribute,
RepresentationMapping,
SkyCoord,
TimeAttribute,
UnitSphericalRepresentation,
frame_transform_graph,
)
from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
from astropy.time import Time
from numpy import broadcast_arrays
[docs]
def cart2pol(x, y):
rho = np.sqrt(x**2 + y**2)
phi = np.arctan2(y, x)
return (rho, phi)
[docs]
def pol2cart(rho, phi):
x = rho * np.cos(phi)
y = rho * np.sin(phi)
return (x, y)
[docs]
def camera_to_shower_coordinates(x, y, cog_x, cog_y, psi):
"""
TODO: stack coordinates and project in // as in pyhiperta
Return longitudinal and transverse coordinates for x and y
for a given set of hillas parameters
Parameters
----------
x: u.Quantity[length]
x coordinate in camera coordinates
y: u.Quantity[length]
y coordinate in camera coordinates
cog_x: u.Quantity[length]
x coordinate of center of gravity
cog_y: u.Quantity[length]
y coordinate of center of gravity
psi: Angle
orientation angle
Returns
-------
longitudinal: astropy.units.Quantity
longitudinal coordinates (along the shower axis)
transverse: astropy.units.Quantity
transverse coordinates (perpendicular to the shower axis)
"""
cos_psi = np.cos(psi)
sin_psi = np.sin(psi)
delta_x = x - cog_x
delta_y = y - cog_y
longi = delta_x * cos_psi + delta_y * sin_psi
trans = delta_x * -sin_psi + delta_y * cos_psi
return longi, trans
# TODO: if this is a copy, can't we use the original ?
[docs]
class PlanarRepresentation(BaseRepresentation):
"""
Representation of a point in a 2D plane.
This is essentially a copy of the Cartesian representation used
in astropy.
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)])
[docs]
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 {}".format(self.attr_classes["x"].__name__))
if not isinstance(y, self.attr_classes["y"]):
raise TypeError("y should be a {}".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))
[docs]
@classmethod
def from_cartesian(cls, cartesian):
return cls(x=cartesian.x, y=cartesian.y)
[docs]
def to_cartesian(self):
return CartesianRepresentation(x=self._x, y=self._y, z=0 * self._x.unit)
@property
def components(self):
return "x", "y"
[docs]
class TelescopeFrame(BaseCoordinateFrame):
"""
Telescope coordinate frame.
A Frame using a UnitSphericalRepresentation.
This is basically the same as a HorizonCoordinate, but the
origin is at the telescope's pointing direction.
This is used to specify coordinates in the field of view of a
telescope that is independent of the optical properties of the telescope.
``fov_lon`` is aligned with azimuth and ``fov_lat`` is aligned with altitude
of the horizontal coordinate frame as implemented in ``astropy.coordinates.AltAz``.
This is what astropy calls a SkyOffsetCoordinate.
Attributes
----------
telescope_pointing: SkyCoord[AltAz]
Coordinate of the telescope pointing in AltAz
obstime: Tiem
Observation time
location: EarthLocation
Location of the telescope
"""
frame_specific_representation_info = {
UnitSphericalRepresentation: [
RepresentationMapping("lon", "fov_lon"),
RepresentationMapping("lat", "fov_lat"),
]
}
default_representation = UnitSphericalRepresentation
telescope_pointing = CoordinateAttribute(default=None, frame=AltAz)
obstime = TimeAttribute(default=None)
location = EarthLocationAttribute(default=None)
[docs]
def __init__(self, *args, **kwargs):
super().__init__(*args, **kwargs)
# make sure telescope coordinate is in range [-180°, 180°]
if isinstance(self._data, UnitSphericalRepresentation):
self._data.lon.wrap_angle = Angle(180, unit=u.deg)
[docs]
@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, TelescopeFrame)
def telescope_to_telescope(from_telescope_coord, to_telescope_frame):
"""Transform between two skyoffset frames."""
intermediate_from = from_telescope_coord.transform_to(from_telescope_coord.telescope_pointing)
intermediate_to = intermediate_from.transform_to(to_telescope_frame.telescope_pointing)
return intermediate_to.transform_to(to_telescope_frame)
[docs]
@frame_transform_graph.transform(DynamicMatrixTransform, AltAz, TelescopeFrame)
def altaz_to_telescope(altaz_coord, telescope_frame):
"""Convert a reference coordinate to an sky offset frame."""
# Define rotation matrices along the position angle vector, and
# relative to the telescope_pointing.
telescope_pointing = telescope_frame.telescope_pointing.represent_as(UnitSphericalRepresentation)
mat1 = rotation_matrix(-telescope_pointing.lat, "y")
mat2 = rotation_matrix(telescope_pointing.lon, "z")
return mat1 @ mat2
[docs]
@frame_transform_graph.transform(DynamicMatrixTransform, TelescopeFrame, AltAz)
def telescope_to_altaz(telescope_coord, altaz_frame):
"""Convert an sky offset frame coordinate to the reference frame"""
# use the forward transform, but just invert it
mat = altaz_to_telescope(altaz_frame, telescope_coord)
# transpose is the inverse because mat is a rotation matrix
return matrix_transpose(mat)
[docs]
class CameraFrame(BaseCoordinateFrame):
"""
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 H.E.S.S., starting at the horizon,
the telescope is pointed to magnetic north in azimuth and then up to zenith.
Now, x points north and y points west, so in this orientation, the
camera coordinates line up with the CORSIKA ground coordinate system.
MAGIC and FACT use a different camera coordinate system:
Standing at the dish, looking at the camera, x points right, y points up.
To transform MAGIC/FACT to ctapipe, 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
"""
default_representation = PlanarRepresentation
focal_length = QuantityAttribute(default=0, unit=u.m)
rotation = QuantityAttribute(default=0 * u.deg, unit=u.rad)
telescope_pointing = CoordinateAttribute(frame=AltAz, default=None)
obstime = TimeAttribute(default=None)
location = EarthLocationAttribute(default=None)
[docs]
@frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame)
def camera_to_telescope(camera_coord, telescope_frame):
"""
Transformation between CameraFrame and TelescopeFrame.
Is called when a SkyCoord is transformed from CameraFrame into TelescopeFrame
"""
x_pos = camera_coord.cartesian.x
y_pos = camera_coord.cartesian.y
rot = camera_coord.rotation
if rot == 0: # if no rotation applied save a few cycles
x_rotated = x_pos
y_rotated = y_pos
else:
cosrot = np.cos(rot)
sinrot = np.sin(rot)
x_rotated = x_pos * cosrot - y_pos * sinrot
y_rotated = x_pos * sinrot + y_pos * cosrot
focal_length = camera_coord.focal_length
if focal_length.shape == () and focal_length.value == 0:
raise ValueError("Coordinate in CameraFrame is missing focal_length information")
# this assumes an equidistant mapping function of the telescope optics
# or a small angle approximation of most other mapping functions
# this could be replaced by actually defining the mapping function
# as an Attribute of `CameraFrame` that maps f(r, focal_length) -> theta,
# where theta is the angle to the optical axis and r is the distance
# to the camera center in the focal plane
fov_lat = u.Quantity((x_rotated / focal_length).to_value(u.dimensionless_unscaled), u.rad, copy=False)
fov_lon = u.Quantity((y_rotated / focal_length).to_value(u.dimensionless_unscaled), u.rad, copy=False)
representation = UnitSphericalRepresentation(lat=fov_lat, lon=fov_lon)
return telescope_frame.realize_frame(representation)
[docs]
@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, CameraFrame)
def telescope_to_camera(telescope_coord, camera_frame):
"""
Transformation between TelescopeFrame and CameraFrame
Is called when a SkyCoord is transformed from TelescopeFrame into CameraFrame
"""
x_pos = telescope_coord.fov_lat
y_pos = telescope_coord.fov_lon
# reverse the rotation applied to get to this system
rot = -camera_frame.rotation
if rot.value == 0.0: # if no rotation applied save a few cycles
x_rotated = x_pos
y_rotated = y_pos
else: # or else rotate all positions around the camera centre
cosrot = np.cos(rot)
sinrot = np.sin(rot)
x_rotated = x_pos * cosrot - y_pos * sinrot
y_rotated = x_pos * sinrot + y_pos * cosrot
focal_length = camera_frame.focal_length
if focal_length.shape == () and focal_length.value == 0:
raise ValueError("CameraFrame is missing focal_length information")
# this assumes an equidistant mapping function of the telescope optics
# or a small angle approximation of most other mapping functions
# this could be replaced by actually defining the mapping function
# as an Attribute of `CameraFrame` that maps f(theta, focal_length) -> r,
# where theta is the angle to the optical axis and r is the distance
# to the camera center in the focal plane
x_rotated = x_rotated.to_value(u.rad) * focal_length
y_rotated = y_rotated.to_value(u.rad) * focal_length
representation = CartesianRepresentation(x_rotated, y_rotated, 0 * u.m)
return camera_frame.realize_frame(representation)
[docs]
def sky_to_camera(alt, az, focal, pointing_alt, pointing_az):
"""
Coordinate transform from aky position (alt, az) (in angles)
to camera coordinates (x, y) in distance.
Parameters
----------
alt: astropy Quantity
az: astropy Quantity
focal: astropy Quantity
pointing_alt: pointing altitude in angle unit
pointing_az: pointing altitude in angle unit
Returns
-------
camera frame: `astropy.coordinates.sky_coordinate.SkyCoord`
"""
# TODO: what are all those hard-coded values ?!?
LST1_LOCATION = EarthLocation(
lon=-17.89149701 * u.deg,
lat=28.76152611 * u.deg,
# height of central pin + distance from pin to elevation axis
height=2184 * u.m + 15.883 * u.m,
)
horizon_frame = AltAz(location=LST1_LOCATION, obstime=Time("2018-11-01T02:00"))
pointing_direction = SkyCoord(
alt=np.clip(pointing_alt, -90.0 * u.deg, 90.0 * u.deg), az=pointing_az, frame=horizon_frame
)
camera_frame = CameraFrame(focal_length=focal, telescope_pointing=pointing_direction)
event_direction = SkyCoord(alt=np.clip(alt, -90.0 * u.deg, 90.0 * u.deg), az=az, frame=horizon_frame)
camera_pos = event_direction.transform_to(camera_frame)
return camera_pos
[docs]
def camera_to_altaz(pos_x, pos_y, focal, pointing_alt, pointing_az, obstime=None):
"""
Compute camera to Horizontal frame (Altitude-Azimuth system).
For MC assume the default ObsTime.
Parameters
----------
pos_x: `~astropy.units.Quantity`
X coordinate in camera (distance)
pos_y: `~astropy.units.Quantity`
Y coordinate in camera (distance)
focal: `~astropy.units.Quantity`
telescope focal (distance)
pointing_alt: `~astropy.units.Quantity`
pointing altitude in angle unit
pointing_az: `~astropy.units.Quantity`
pointing altitude in angle unit
obstime: `~astropy.time.Time`
Returns
-------
sky frame: `astropy.coordinates.SkyCoord`
in AltAz frame
Examples
--------
>>> import astropy.units as u
>>> import numpy as np
>>> pos_x = np.array([0, 0]) * u.m
>>> pos_y = np.array([0, 0]) * u.m
>>> focal = 28 * u.m
>>> pointing_alt = np.array([1.0, 1.0]) * u.rad
>>> pointing_az = np.array([0.2, 0.5]) * u.rad
>>> sky_coords = utils.camera_to_altaz(pos_x, pos_y, focal, pointing_alt, pointing_az)
"""
# TODO: better logging/warning !
if not obstime:
warnings.warn("No time given. To be use only for MC data.")
# TODO: what are all those hard-coded values ?!?
LST1_LOCATION = EarthLocation(
lon=-17.89149701 * u.deg,
lat=28.76152611 * u.deg,
# height of central pin + distance from pin to elevation axis
height=2184 * u.m + 15.883 * u.m,
)
horizon_frame = AltAz(location=LST1_LOCATION, obstime=obstime)
pointing_direction = SkyCoord(
alt=np.clip(pointing_alt, -90.0 * u.deg, 90.0 * u.deg), az=pointing_az, frame=horizon_frame
)
camera_frame = CameraFrame(focal_length=focal, telescope_pointing=pointing_direction)
camera_coord = SkyCoord(pos_x, pos_y, frame=camera_frame)
horizon = camera_coord.transform_to(horizon_frame)
return horizon