Coverage for rta_to_lstchain_formator/utils/coordinates.py: 45%

1from collections import OrderedDict 

2 

3import astropy.units as u 

4import numpy as np 

5from astropy.coordinates import ( 

6 AltAz, 

7 BaseCoordinateFrame, 

8 BaseRepresentation, 

9 CartesianRepresentation, 

10 CoordinateAttribute, 

11 EarthLocation, 

12 EarthLocationAttribute, 

13 QuantityAttribute, 

14 SkyCoord, 

15 TimeAttribute, 

16) 

17from astropy.time import Time 

18from numpy import broadcast_arrays 

19 

20 

21# copied from https://github.com/cta-observatory/ctapipe/blob/v0.19.2/ctapipe/coordinates/representation.py#L13 

22# to avoid dependency on entire ctapipe 

23class PlanarRepresentation(BaseRepresentation): 

24 """ 

25 Representation of a point in a 2D plane. 

26 This is essentially a copy of the Cartesian representation used 

27 in astropy. 

28 

29 Parameters 

30 ---------- 

31 

32 x, y : `~astropy.units.Quantity` 

33 The x and y coordinates of the point(s). If ``x`` and ``y``have 

34 different shapes, they should be broadcastable. 

35 

36 copy : bool, optional 

37 If True arrays will be copied rather than referenced. 

38 

39 """ 

40 

41 attr_classes = OrderedDict([("x", u.Quantity), ("y", u.Quantity)]) 

42 

43 def __init__(self, x, y, copy=True, **kwargs): 

44 if x is None or y is None: 

45 raise ValueError("x and y are required to instantiate CartesianRepresentation") 

46 

47 if not isinstance(x, self.attr_classes["x"]): 

48 raise TypeError("x should be a {}".format(self.attr_classes["x"].__name__)) 

49 

50 if not isinstance(y, self.attr_classes["y"]): 

51 raise TypeError("y should be a {}".format(self.attr_classes["y"].__name__)) 

52 

53 x = self.attr_classes["x"](x, copy=copy) 

54 y = self.attr_classes["y"](y, copy=copy) 

55 

56 if not (x.unit.physical_type == y.unit.physical_type): 

57 raise u.UnitsError("x and y should have matching physical types") 

58 

59 try: 

60 x, y = broadcast_arrays(x, y, subok=True) 

61 except ValueError: 

62 raise ValueError("Input parameters x and y cannot be broadcast") 

63 

64 self._x = x 

65 self._y = y 

66 self._differentials = {} 

67 

68 @property 

69 def x(self): 

70 """ 

71 The x component of the point(s). 

72 """ 

73 return self._x 

74 

75 @property 

76 def y(self): 

77 """ 

78 The y component of the point(s). 

79 """ 

80 return self._y 

81 

82 @property 

83 def xy(self): 

84 return u.Quantity((self._x, self._y)) 

85 

86 @classmethod 

87 def from_cartesian(cls, cartesian): 

88 return cls(x=cartesian.x, y=cartesian.y) 

89 

90 def to_cartesian(self): 

91 return CartesianRepresentation(x=self._x, y=self._y, z=0 * self._x.unit) 

92 

93 @property 

94 def components(self): 

95 return "x", "y" 

96 

97 

98# Copied from https://github.com/cta-observatory/ctapipe/blob/v0.19.2/ctapipe/coordinates/camera_frame.py#L44 

99# to avoid dependency on entire ctapipe 

100class CameraFrame(BaseCoordinateFrame): 

101 """ 

102 Camera coordinate frame. 

103 

104 The camera frame is a 2d cartesian frame, 

105 describing position of objects in the focal plane of the telescope. 

106 

107 The frame is defined as in H.E.S.S., starting at the horizon, 

108 the telescope is pointed to magnetic north in azimuth and then up to zenith. 

109 

110 Now, x points north and y points west, so in this orientation, the 

111 camera coordinates line up with the CORSIKA ground coordinate system. 

112 

113 MAGIC and FACT use a different camera coordinate system: 

114 Standing at the dish, looking at the camera, x points right, y points up. 

115 To transform MAGIC/FACT to ctapipe, do x' = -y, y' = -x. 

116 

117 Attributes 

118 ---------- 

119 

120 focal_length : u.Quantity[length] 

121 Focal length of the telescope as a unit quantity (usually meters) 

122 rotation : u.Quantity[angle] 

123 Rotation angle of the camera (0 deg in most cases) 

124 telescope_pointing : SkyCoord[AltAz] 

125 Pointing direction of the telescope as SkyCoord in AltAz 

126 obstime : Time 

127 Observation time 

128 location : EarthLocation 

129 location of the telescope 

130 """ 

131 

132 default_representation = PlanarRepresentation 

133 

134 focal_length = QuantityAttribute(default=0, unit=u.m) 

135 rotation = QuantityAttribute(default=0 * u.deg, unit=u.rad) 

136 telescope_pointing = CoordinateAttribute(frame=AltAz, default=None) 

137 

138 obstime = TimeAttribute(default=None) 

139 location = EarthLocationAttribute(default=None) 

140 

141 

142# Copied from https://github.com/cta-observatory/ctapipe_io_lst/blob/v0.22.6/src/ctapipe_io_lst/constants.py#L30 

143# to avoid dependency on entire lstchain 

144#: location of lst-1 as `~astropy.coordinates.EarthLocation` 

145#: Taken from Abelardo's Coordinates of LST-1 & MAGIC presentation 

146#: https://redmine.cta-observatory.org/attachments/65827 

147LST1_LOCATION = EarthLocation( 

148 lon=-17.89149701 * u.deg, 

149 lat=28.76152611 * u.deg, 

150 # height of central pin + distance from pin to elevation axis 

151 height=2184 * u.m + 15.883 * u.m, 

152) 

153 

154 

155# copied from https://github.com/cta-observatory/cta-lstchain/blob/v0.10.7/lstchain/reco/utils.py#L608 

156# to avoid dependency on entire lstchain 

157def clip_alt(alt): 

158 """ 

159 Make sure altitude is not larger than 90 deg (it happens in some MC files for zenith=0), 

160 to keep astropy happy 

161 """ 

162 return np.clip(alt, -90.0 * u.deg, 90.0 * u.deg) 

163 

164 

165# copied from https://github.com/cta-observatory/cta-lstchain/blob/v0.10.7/lstchain/reco/utils.py#L316 

166# to avoid dependency on entire lstchain 

167def sky_to_camera(alt, az, focal, pointing_alt, pointing_az): 

168 """ 

169 Coordinate transform from aky position (alt, az) (in angles) 

170 to camera coordinates (x, y) in distance. 

171 

172 Parameters 

173 ---------- 

174 alt: astropy Quantity 

175 az: astropy Quantity 

176 focal: astropy Quantity 

177 pointing_alt: pointing altitude in angle unit 

178 pointing_az: pointing altitude in angle unit 

179 

180 Returns 

181 ------- 

182 camera frame: `astropy.coordinates.sky_coordinate.SkyCoord` 

183 """ 

184 pointing_direction = SkyCoord( 

185 alt=clip_alt(pointing_alt), 

186 az=pointing_az, 

187 frame=AltAz(location=LST1_LOCATION, obstime=Time("2018-11-01T02:00")), 

188 ) 

189 

190 camera_frame = CameraFrame(focal_length=focal, telescope_pointing=pointing_direction) 

191 

192 event_direction = SkyCoord( 

193 alt=clip_alt(alt), az=az, frame=AltAz(location=LST1_LOCATION, obstime=Time("2018-11-01T02:00")) 

194 ) 

195 

196 camera_pos = event_direction.transform_to(camera_frame) 

197 

198 return camera_pos