Coverage for rta_reconstruction/utils/coordinates.py: 68%
149 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-11-02 09:59 +0000
1import warnings
2from collections import OrderedDict
4import astropy.units as u
5import numpy as np
6from astropy.coordinates import (
7 AltAz,
8 Angle,
9 BaseCoordinateFrame,
10 BaseRepresentation,
11 CartesianRepresentation,
12 CoordinateAttribute,
13 DynamicMatrixTransform,
14 EarthLocation,
15 EarthLocationAttribute,
16 FunctionTransform,
17 QuantityAttribute,
18 RepresentationMapping,
19 SkyCoord,
20 TimeAttribute,
21 UnitSphericalRepresentation,
22 frame_transform_graph,
23)
24from astropy.coordinates.matrix_utilities import matrix_transpose, rotation_matrix
25from astropy.time import Time
26from numpy import broadcast_arrays
29def cart2pol(x, y):
30 rho = np.sqrt(x**2 + y**2)
31 phi = np.arctan2(y, x)
32 return (rho, phi)
35def pol2cart(rho, phi):
36 x = rho * np.cos(phi)
37 y = rho * np.sin(phi)
38 return (x, y)
41def camera_to_shower_coordinates(x, y, cog_x, cog_y, psi):
42 """
43 TODO: stack coordinates and project in // as in pyhiperta
44 Return longitudinal and transverse coordinates for x and y
45 for a given set of hillas parameters
47 Parameters
48 ----------
49 x: u.Quantity[length]
50 x coordinate in camera coordinates
51 y: u.Quantity[length]
52 y coordinate in camera coordinates
53 cog_x: u.Quantity[length]
54 x coordinate of center of gravity
55 cog_y: u.Quantity[length]
56 y coordinate of center of gravity
57 psi: Angle
58 orientation angle
60 Returns
61 -------
62 longitudinal: astropy.units.Quantity
63 longitudinal coordinates (along the shower axis)
64 transverse: astropy.units.Quantity
65 transverse coordinates (perpendicular to the shower axis)
66 """
67 cos_psi = np.cos(psi)
68 sin_psi = np.sin(psi)
70 delta_x = x - cog_x
71 delta_y = y - cog_y
73 longi = delta_x * cos_psi + delta_y * sin_psi
74 trans = delta_x * -sin_psi + delta_y * cos_psi
76 return longi, trans
79# TODO: if this is a copy, can't we use the original ?
80class PlanarRepresentation(BaseRepresentation):
81 """
82 Representation of a point in a 2D plane.
83 This is essentially a copy of the Cartesian representation used
84 in astropy.
86 Parameters
87 ----------
89 x, y : `~astropy.units.Quantity`
90 The x and y coordinates of the point(s). If ``x`` and ``y``have
91 different shapes, they should be broadcastable.
93 copy : bool, optional
94 If True arrays will be copied rather than referenced.
96 """
98 attr_classes = OrderedDict([("x", u.Quantity), ("y", u.Quantity)])
100 def __init__(self, x, y, copy=True, **kwargs):
101 if x is None or y is None:
102 raise ValueError("x and y are required to instantiate CartesianRepresentation")
104 if not isinstance(x, self.attr_classes["x"]):
105 raise TypeError("x should be a {}".format(self.attr_classes["x"].__name__))
107 if not isinstance(y, self.attr_classes["y"]):
108 raise TypeError("y should be a {}".format(self.attr_classes["y"].__name__))
110 x = self.attr_classes["x"](x, copy=copy)
111 y = self.attr_classes["y"](y, copy=copy)
113 if not (x.unit.physical_type == y.unit.physical_type):
114 raise u.UnitsError("x and y should have matching physical types")
116 try:
117 x, y = broadcast_arrays(x, y, subok=True)
118 except ValueError:
119 raise ValueError("Input parameters x and y cannot be broadcast")
121 self._x = x
122 self._y = y
123 self._differentials = {}
125 @property
126 def x(self):
127 """
128 The x component of the point(s).
129 """
130 return self._x
132 @property
133 def y(self):
134 """
135 The y component of the point(s).
136 """
137 return self._y
139 @property
140 def xy(self):
141 return u.Quantity((self._x, self._y))
143 @classmethod
144 def from_cartesian(cls, cartesian):
145 return cls(x=cartesian.x, y=cartesian.y)
147 def to_cartesian(self):
148 return CartesianRepresentation(x=self._x, y=self._y, z=0 * self._x.unit)
150 @property
151 def components(self):
152 return "x", "y"
155class TelescopeFrame(BaseCoordinateFrame):
156 """
157 Telescope coordinate frame.
159 A Frame using a UnitSphericalRepresentation.
160 This is basically the same as a HorizonCoordinate, but the
161 origin is at the telescope's pointing direction.
163 This is used to specify coordinates in the field of view of a
164 telescope that is independent of the optical properties of the telescope.
166 ``fov_lon`` is aligned with azimuth and ``fov_lat`` is aligned with altitude
167 of the horizontal coordinate frame as implemented in ``astropy.coordinates.AltAz``.
169 This is what astropy calls a SkyOffsetCoordinate.
171 Attributes
172 ----------
174 telescope_pointing: SkyCoord[AltAz]
175 Coordinate of the telescope pointing in AltAz
176 obstime: Tiem
177 Observation time
178 location: EarthLocation
179 Location of the telescope
180 """
182 frame_specific_representation_info = {
183 UnitSphericalRepresentation: [
184 RepresentationMapping("lon", "fov_lon"),
185 RepresentationMapping("lat", "fov_lat"),
186 ]
187 }
188 default_representation = UnitSphericalRepresentation
190 telescope_pointing = CoordinateAttribute(default=None, frame=AltAz)
192 obstime = TimeAttribute(default=None)
193 location = EarthLocationAttribute(default=None)
195 def __init__(self, *args, **kwargs):
196 super().__init__(*args, **kwargs)
198 # make sure telescope coordinate is in range [-180°, 180°]
199 if isinstance(self._data, UnitSphericalRepresentation):
200 self._data.lon.wrap_angle = Angle(180, unit=u.deg)
203@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, TelescopeFrame)
204def telescope_to_telescope(from_telescope_coord, to_telescope_frame):
205 """Transform between two skyoffset frames."""
207 intermediate_from = from_telescope_coord.transform_to(from_telescope_coord.telescope_pointing)
208 intermediate_to = intermediate_from.transform_to(to_telescope_frame.telescope_pointing)
209 return intermediate_to.transform_to(to_telescope_frame)
212@frame_transform_graph.transform(DynamicMatrixTransform, AltAz, TelescopeFrame)
213def altaz_to_telescope(altaz_coord, telescope_frame):
214 """Convert a reference coordinate to an sky offset frame."""
216 # Define rotation matrices along the position angle vector, and
217 # relative to the telescope_pointing.
218 telescope_pointing = telescope_frame.telescope_pointing.represent_as(UnitSphericalRepresentation)
219 mat1 = rotation_matrix(-telescope_pointing.lat, "y")
220 mat2 = rotation_matrix(telescope_pointing.lon, "z")
221 return mat1 @ mat2
224@frame_transform_graph.transform(DynamicMatrixTransform, TelescopeFrame, AltAz)
225def telescope_to_altaz(telescope_coord, altaz_frame):
226 """Convert an sky offset frame coordinate to the reference frame"""
228 # use the forward transform, but just invert it
229 mat = altaz_to_telescope(altaz_frame, telescope_coord)
230 # transpose is the inverse because mat is a rotation matrix
231 return matrix_transpose(mat)
234class CameraFrame(BaseCoordinateFrame):
235 """
236 Camera coordinate frame.
238 The camera frame is a 2d cartesian frame,
239 describing position of objects in the focal plane of the telescope.
241 The frame is defined as in H.E.S.S., starting at the horizon,
242 the telescope is pointed to magnetic north in azimuth and then up to zenith.
244 Now, x points north and y points west, so in this orientation, the
245 camera coordinates line up with the CORSIKA ground coordinate system.
247 MAGIC and FACT use a different camera coordinate system:
248 Standing at the dish, looking at the camera, x points right, y points up.
249 To transform MAGIC/FACT to ctapipe, do x' = -y, y' = -x.
251 Attributes
252 ----------
254 focal_length : u.Quantity[length]
255 Focal length of the telescope as a unit quantity (usually meters)
256 rotation : u.Quantity[angle]
257 Rotation angle of the camera (0 deg in most cases)
258 telescope_pointing : SkyCoord[AltAz]
259 Pointing direction of the telescope as SkyCoord in AltAz
260 obstime : Time
261 Observation time
262 location : EarthLocation
263 location of the telescope
264 """
266 default_representation = PlanarRepresentation
268 focal_length = QuantityAttribute(default=0, unit=u.m)
269 rotation = QuantityAttribute(default=0 * u.deg, unit=u.rad)
270 telescope_pointing = CoordinateAttribute(frame=AltAz, default=None)
272 obstime = TimeAttribute(default=None)
273 location = EarthLocationAttribute(default=None)
276@frame_transform_graph.transform(FunctionTransform, CameraFrame, TelescopeFrame)
277def camera_to_telescope(camera_coord, telescope_frame):
278 """
279 Transformation between CameraFrame and TelescopeFrame.
280 Is called when a SkyCoord is transformed from CameraFrame into TelescopeFrame
281 """
282 x_pos = camera_coord.cartesian.x
283 y_pos = camera_coord.cartesian.y
285 rot = camera_coord.rotation
286 if rot == 0: # if no rotation applied save a few cycles
287 x_rotated = x_pos
288 y_rotated = y_pos
289 else:
290 cosrot = np.cos(rot)
291 sinrot = np.sin(rot)
292 x_rotated = x_pos * cosrot - y_pos * sinrot
293 y_rotated = x_pos * sinrot + y_pos * cosrot
295 focal_length = camera_coord.focal_length
296 if focal_length.shape == () and focal_length.value == 0:
297 raise ValueError("Coordinate in CameraFrame is missing focal_length information")
299 # this assumes an equidistant mapping function of the telescope optics
300 # or a small angle approximation of most other mapping functions
301 # this could be replaced by actually defining the mapping function
302 # as an Attribute of `CameraFrame` that maps f(r, focal_length) -> theta,
303 # where theta is the angle to the optical axis and r is the distance
304 # to the camera center in the focal plane
305 fov_lat = u.Quantity((x_rotated / focal_length).to_value(u.dimensionless_unscaled), u.rad, copy=False)
306 fov_lon = u.Quantity((y_rotated / focal_length).to_value(u.dimensionless_unscaled), u.rad, copy=False)
308 representation = UnitSphericalRepresentation(lat=fov_lat, lon=fov_lon)
310 return telescope_frame.realize_frame(representation)
313@frame_transform_graph.transform(FunctionTransform, TelescopeFrame, CameraFrame)
314def telescope_to_camera(telescope_coord, camera_frame):
315 """
316 Transformation between TelescopeFrame and CameraFrame
318 Is called when a SkyCoord is transformed from TelescopeFrame into CameraFrame
319 """
320 x_pos = telescope_coord.fov_lat
321 y_pos = telescope_coord.fov_lon
322 # reverse the rotation applied to get to this system
323 rot = -camera_frame.rotation
325 if rot.value == 0.0: # if no rotation applied save a few cycles
326 x_rotated = x_pos
327 y_rotated = y_pos
328 else: # or else rotate all positions around the camera centre
329 cosrot = np.cos(rot)
330 sinrot = np.sin(rot)
331 x_rotated = x_pos * cosrot - y_pos * sinrot
332 y_rotated = x_pos * sinrot + y_pos * cosrot
334 focal_length = camera_frame.focal_length
335 if focal_length.shape == () and focal_length.value == 0:
336 raise ValueError("CameraFrame is missing focal_length information")
338 # this assumes an equidistant mapping function of the telescope optics
339 # or a small angle approximation of most other mapping functions
340 # this could be replaced by actually defining the mapping function
341 # as an Attribute of `CameraFrame` that maps f(theta, focal_length) -> r,
342 # where theta is the angle to the optical axis and r is the distance
343 # to the camera center in the focal plane
344 x_rotated = x_rotated.to_value(u.rad) * focal_length
345 y_rotated = y_rotated.to_value(u.rad) * focal_length
347 representation = CartesianRepresentation(x_rotated, y_rotated, 0 * u.m)
349 return camera_frame.realize_frame(representation)
352def sky_to_camera(alt, az, focal, pointing_alt, pointing_az):
353 """
354 Coordinate transform from aky position (alt, az) (in angles)
355 to camera coordinates (x, y) in distance.
357 Parameters
358 ----------
359 alt: astropy Quantity
360 az: astropy Quantity
361 focal: astropy Quantity
362 pointing_alt: pointing altitude in angle unit
363 pointing_az: pointing altitude in angle unit
365 Returns
366 -------
367 camera frame: `astropy.coordinates.sky_coordinate.SkyCoord`
368 """
369 # TODO: what are all those hard-coded values ?!?
370 LST1_LOCATION = EarthLocation(
371 lon=-17.89149701 * u.deg,
372 lat=28.76152611 * u.deg,
373 # height of central pin + distance from pin to elevation axis
374 height=2184 * u.m + 15.883 * u.m,
375 )
376 horizon_frame = AltAz(location=LST1_LOCATION, obstime=Time("2018-11-01T02:00"))
378 pointing_direction = SkyCoord(
379 alt=np.clip(pointing_alt, -90.0 * u.deg, 90.0 * u.deg), az=pointing_az, frame=horizon_frame
380 )
382 camera_frame = CameraFrame(focal_length=focal, telescope_pointing=pointing_direction)
384 event_direction = SkyCoord(alt=np.clip(alt, -90.0 * u.deg, 90.0 * u.deg), az=az, frame=horizon_frame)
386 camera_pos = event_direction.transform_to(camera_frame)
388 return camera_pos
391def camera_to_altaz(pos_x, pos_y, focal, pointing_alt, pointing_az, obstime=None):
392 """
393 Compute camera to Horizontal frame (Altitude-Azimuth system).
394 For MC assume the default ObsTime.
396 Parameters
397 ----------
398 pos_x: `~astropy.units.Quantity`
399 X coordinate in camera (distance)
400 pos_y: `~astropy.units.Quantity`
401 Y coordinate in camera (distance)
402 focal: `~astropy.units.Quantity`
403 telescope focal (distance)
404 pointing_alt: `~astropy.units.Quantity`
405 pointing altitude in angle unit
406 pointing_az: `~astropy.units.Quantity`
407 pointing altitude in angle unit
408 obstime: `~astropy.time.Time`
410 Returns
411 -------
412 sky frame: `astropy.coordinates.SkyCoord`
413 in AltAz frame
415 Examples
416 --------
417 >>> import astropy.units as u
418 >>> import numpy as np
419 >>> pos_x = np.array([0, 0]) * u.m
420 >>> pos_y = np.array([0, 0]) * u.m
421 >>> focal = 28 * u.m
422 >>> pointing_alt = np.array([1.0, 1.0]) * u.rad
423 >>> pointing_az = np.array([0.2, 0.5]) * u.rad
424 >>> sky_coords = utils.camera_to_altaz(pos_x, pos_y, focal, pointing_alt, pointing_az)
425 """
426 # TODO: better logging/warning !
427 if not obstime:
428 warnings.warn("No time given. To be use only for MC data.")
430 # TODO: what are all those hard-coded values ?!?
431 LST1_LOCATION = EarthLocation(
432 lon=-17.89149701 * u.deg,
433 lat=28.76152611 * u.deg,
434 # height of central pin + distance from pin to elevation axis
435 height=2184 * u.m + 15.883 * u.m,
436 )
437 horizon_frame = AltAz(location=LST1_LOCATION, obstime=obstime)
439 pointing_direction = SkyCoord(
440 alt=np.clip(pointing_alt, -90.0 * u.deg, 90.0 * u.deg), az=pointing_az, frame=horizon_frame
441 )
443 camera_frame = CameraFrame(focal_length=focal, telescope_pointing=pointing_direction)
445 camera_coord = SkyCoord(pos_x, pos_y, frame=camera_frame)
447 horizon = camera_coord.transform_to(horizon_frame)
449 return horizon