from collections import OrderedDict
import astropy.units as u
import numpy as np
from astropy.coordinates import (
AltAz,
BaseCoordinateFrame,
BaseRepresentation,
CartesianRepresentation,
CoordinateAttribute,
EarthLocation,
EarthLocationAttribute,
QuantityAttribute,
SkyCoord,
TimeAttribute,
)
from astropy.time import Time
from numpy import broadcast_arrays
# copied from https://github.com/cta-observatory/ctapipe/blob/v0.19.2/ctapipe/coordinates/representation.py#L13
# to avoid dependency on entire ctapipe
[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"
# Copied from https://github.com/cta-observatory/ctapipe/blob/v0.19.2/ctapipe/coordinates/camera_frame.py#L44
# to avoid dependency on entire ctapipe
[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)
# Copied from https://github.com/cta-observatory/ctapipe_io_lst/blob/v0.22.6/src/ctapipe_io_lst/constants.py#L30
# to avoid dependency on entire lstchain
#: location of lst-1 as `~astropy.coordinates.EarthLocation`
#: Taken from Abelardo's Coordinates of LST-1 & MAGIC presentation
#: https://redmine.cta-observatory.org/attachments/65827
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,
)
# copied from https://github.com/cta-observatory/cta-lstchain/blob/v0.10.7/lstchain/reco/utils.py#L608
# to avoid dependency on entire lstchain
[docs]
def clip_alt(alt):
"""
Make sure altitude is not larger than 90 deg (it happens in some MC files for zenith=0),
to keep astropy happy
"""
return np.clip(alt, -90.0 * u.deg, 90.0 * u.deg)
# copied from https://github.com/cta-observatory/cta-lstchain/blob/v0.10.7/lstchain/reco/utils.py#L316
# to avoid dependency on entire lstchain
[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`
"""
pointing_direction = SkyCoord(
alt=clip_alt(pointing_alt),
az=pointing_az,
frame=AltAz(location=LST1_LOCATION, obstime=Time("2018-11-01T02:00")),
)
camera_frame = CameraFrame(focal_length=focal, telescope_pointing=pointing_direction)
event_direction = SkyCoord(
alt=clip_alt(alt), az=az, frame=AltAz(location=LST1_LOCATION, obstime=Time("2018-11-01T02:00"))
)
camera_pos = event_direction.transform_to(camera_frame)
return camera_pos