Source code for ctapipe.coordinates.ground_frames

"""This module defines the important coordinate systems to be used in
reconstruction with the CTA pipeline and the transformations between
this different systems. Frames and transformations are defined using
the astropy.coordinates framework. This module defines transformations
for ground based cartesian and planar systems.

For examples on usage see examples/

This code is based on the coordinate transformations performed in the
read_hess code


- Tests Tests Tests!
import astropy.units as u
import numpy as np
from astropy.coordinates import (
from astropy.units.quantity import Quantity

__all__ = [

def _altaz_to_earthlocation(altaz):
    local_itrs = altaz.transform_to(ITRS(location=altaz.location))
    itrs = ITRS(local_itrs.cartesian + altaz.location.get_itrs().cartesian)
    return itrs.earth_location

def _earthlocation_to_altaz(location, reference_location):
    # See
    # for why this is necessary and we cannot just do
    # `get_itrs().transform_to(AltAz())`
    itrs_cart = location.get_itrs().cartesian
    itrs_ref_cart = reference_location.get_itrs().cartesian
    local_itrs = ITRS(itrs_cart - itrs_ref_cart, location=reference_location)
    return local_itrs.transform_to(AltAz(location=reference_location))

[docs]class GroundFrame(BaseCoordinateFrame): """Ground coordinate frame. The ground coordinate frame is a simple cartesian frame describing the 3 dimensional position of objects compared to the array ground level in relation to the nomial centre of the array. Typically this frame will be used for describing the position on telescopes and equipment. In this frame x points north, y points west and z is meters above array center. Frame attributes: None """ default_representation = CartesianRepresentation reference_location = EarthLocationAttribute()
[docs] def to_earth_location(self): """Convert this GroundFrame coordinate into an `astropy.coordinates.EarthLocation` This requires that the ``reference_location`` is set. """ # in astropy, x points north, y points east, so we need a minus for y. cart = CartesianRepresentation(self.x, -self.y, self.z) altaz = AltAz(cart, location=self.reference_location) return _altaz_to_earthlocation(altaz)
[docs] @classmethod def from_earth_location(cls, location, reference_location): """ Convert `~astropy.coordinates.EarthLocation` into Groundframe. Parameters ---------- location : astropy.coordinates.EarthLocation The location to convert reference_location : astropy.coordinates.EarthLocation Reference location for the GroundFrame Returns ------- ground_frame : GroundFrame EarthLocation converted to GroundFrame """ altaz = _earthlocation_to_altaz(location, reference_location) x, y, z = # in astropy, x points north, y points east, so we need a minus for y. return GroundFrame(x=x, y=-y, z=z, reference_location=reference_location)
@property def observation_level(self): """Height of the reference location""" return self.reference_location.height
[docs]class EastingNorthingFrame(BaseCoordinateFrame): """GroundFrame but in standard Easting/Northing coordinates instead of SimTel/Corsika conventions Frame attributes: None """ default_representation = CartesianRepresentation reference_location = EarthLocationAttribute() frame_specific_representation_info = { CartesianRepresentation: [ RepresentationMapping("x", "easting"), RepresentationMapping("y", "northing"), RepresentationMapping("z", "height"), ] }
[docs] def to_earth_location(self): # in astropy, x points north, y points east cart = CartesianRepresentation(self.northing, self.easting, self.height) altaz = AltAz(cart, location=self.reference_location) return _altaz_to_earthlocation(altaz)
[docs] @classmethod def from_earth_location(cls, location, reference_location): altaz = _earthlocation_to_altaz(location, reference_location) x, y, z = # in astropy, x points north, y points east return GroundFrame( northing=x, easting=y, height=z, reference_location=reference_location, )
@property def observation_level(self): return self.reference_location.height
[docs]class TiltedGroundFrame(BaseCoordinateFrame): """Tilted ground coordinate frame. The tilted ground coordinate frame is a cartesian system describing the 2 dimensional projected positions of objects in a tilted plane described by pointing_direction Typically this frame will be used for the reconstruction of the shower core position Frame attributes: * ``pointing_direction`` Alt,Az direction of the tilted reference plane """ default_representation = CartesianRepresentation # Pointing direction of the tilted system (alt,az), # could be the telescope pointing direction or the reconstructed shower # direction pointing_direction = CoordinateAttribute(default=None, frame=AltAz)
def _get_shower_trans_matrix(azimuth, altitude, inverse=False): """Get Transformation matrix for conversion from the ground system to the Tilted system and back again (This function is directly lifted from read_hess, probably could be streamlined using python functionality) Parameters ---------- azimuth: float or ndarray Azimuth angle in radians of the tilted system used altitude: float or ndarray Altitude angle in radiuan of the tilted system used Returns ------- trans: 3x3 ndarray transformation matrix """ cos_z = np.sin(altitude) # this is the same as np.cos(zenith) but faster sin_z = np.cos(altitude) cos_az = np.cos(azimuth) sin_az = np.sin(azimuth) trans = np.array( [ [cos_z * cos_az, -cos_z * sin_az, -sin_z], [sin_az, cos_az, np.zeros_like(sin_z)], [sin_z * cos_az, -sin_z * sin_az, cos_z], ], dtype=np.float64, ) if inverse: return np.swapaxes(trans, 0, 1) return trans def _get_xyz(coord): """ Essentially the same as, but much faster by avoiding some astropy bottlenecks. """ # this is a speed optimization. Much faster to use data if already a # Cartesian object if isinstance(, CartesianRepresentation): cart = else: cart = coord.cartesian # this is ~5x faster then return u.Quantity([cart.x, cart.y, cart.z]) @frame_transform_graph.transform(FunctionTransform, GroundFrame, TiltedGroundFrame) def ground_to_tilted(ground_coord, tilted_frame): """ Transformation from ground system to tilted ground system Parameters ---------- ground_coord: `astropy.coordinates.SkyCoord` Coordinate in GroundFrame tilted_frame: `ctapipe.coordinates.TiltedFrame` Frame to transform to Returns ------- SkyCoordinate transformed to `tilted_frame` coordinates """ xyz_grd = _get_xyz(ground_coord) # convert to rad first and substract. Faster than .zen altitude = tilted_frame.pointing_direction.alt.to_value(u.rad) azimuth = rotation_matrix = _get_shower_trans_matrix(azimuth, altitude) vec = np.einsum("ij...,j...->i...", rotation_matrix, xyz_grd) representation = CartesianRepresentation(*vec) return tilted_frame.realize_frame(representation) @frame_transform_graph.transform(FunctionTransform, TiltedGroundFrame, GroundFrame) def tilted_to_ground(tilted_coord, ground_frame): """ Transformation from tilted ground system to ground system Parameters ---------- tilted_coord: `astropy.coordinates.SkyCoord` TiltedGroundFrame system ground_frame: `astropy.coordinates.SkyCoord` GroundFrame system Returns ------- GroundFrame coordinates """ xyz_tilt = _get_xyz(tilted_coord) altitude = tilted_coord.pointing_direction.alt.to_value(u.rad) azimuth = rotation_matrix = _get_shower_trans_matrix(azimuth, altitude, inverse=True) vec = np.einsum("ij...,j...->i...", rotation_matrix, xyz_tilt) representation = CartesianRepresentation(*vec) grd = ground_frame.realize_frame(representation) return grd
[docs]def project_to_ground(tilt_system): """Project position in the tilted system onto the ground. This is needed as the standard transformation will return the 3d position of the tilted frame. This projection may ultimately be the standard use case so may be implemented in the tilted to ground transformation Parameters ---------- tilt_system: `astropy.coordinates.SkyCoord` coorinate in the the tilted ground system Returns ------- Projection of tilted system onto the ground (GroundSystem) """ ground_system = tilt_system.transform_to(GroundFrame) unit = ground_system.x.unit x_initial = ground_system.x.value y_initial = ground_system.y.value z_initial = ground_system.z.value trans = _get_shower_trans_matrix(, tilt_system.pointing_direction.alt.to_value(u.rad), ) x_projected = x_initial - trans[2][0] * z_initial / trans[2][2] y_projected = y_initial - trans[2][1] * z_initial / trans[2][2] return GroundFrame( x=u.Quantity(x_projected, unit), y=u.Quantity(y_projected, unit), z=u.Quantity(0, unit), )
@frame_transform_graph.transform(FunctionTransform, GroundFrame, GroundFrame) def ground_to_ground(ground_coords, ground_frame): """Null transform for going from ground to ground, since there are no attributes of the GroundSystem""" return ground_coords # Matrices for transforming between GroundFrame and EastingNorthingFrame NO_OFFSET = CartesianRepresentation(Quantity([0, 0, 0], u.m)) GROUND_TO_EASTNORTH = np.asarray([[0, -1, 0], [1, 0, 0], [0, 0, 1]]) @frame_transform_graph.transform(AffineTransform, GroundFrame, EastingNorthingFrame) def ground_to_easting_northing(ground_coords, eastnorth_frame): """ convert GroundFrame points into eastings/northings for plotting purposes """ return GROUND_TO_EASTNORTH, NO_OFFSET @frame_transform_graph.transform(AffineTransform, EastingNorthingFrame, GroundFrame) def easting_northing_to_ground(eastnorth_coords, ground_frame): """ convert eastings/northings back to GroundFrame """ return GROUND_TO_EASTNORTH.T, NO_OFFSET