Note
Go to the end to download the full example code.
Coordinates usage in ctapipe#
7 import copy
8
9 import astropy.units as u
10 import matplotlib.pyplot as plt
11 from astropy.coordinates import AltAz, EarthLocation, SkyCoord
12 from astropy.time import Time
13
14 from ctapipe.coordinates import CameraFrame, NominalFrame, TelescopeFrame
15 from ctapipe.instrument import SubarrayDescription
16 from ctapipe.io import EventSource
17 from ctapipe.utils import get_dataset_path
18 from ctapipe.visualization import CameraDisplay
19
20 # %matplotlib inline
21
22
23 # make plots and fonts larger
24 plt.rcParams["figure.figsize"] = (12, 8)
25 plt.rcParams["font.size"] = 16
Open test dataset#
33 filename = get_dataset_path("gamma_prod5.simtel.zst")
34 source = EventSource(filename)
35
36 events = [copy.deepcopy(event) for event in source]
37 event = events[4]
38
39 layout = set(source.subarray.tel_ids)
Choose event with LST#
This ensures that the telescope is not “parked” (as it would be in an event where it is not triggered) but is actually pointing to a source.
53 print(f"Telescope with data: {event.r1.tel.keys()}")
54 tel_id = 3
Telescope with data: dict_keys([3, 4, 8, 9, 14, 15, 23, 24, 25, 29, 32, 36, 42, 46, 52, 56, 103, 104, 109, 110, 118, 119, 120, 126, 127, 129, 130, 137, 139, 143, 161])
AltAz#
Pointing direction of telescopes or the origin of a simulated shower are
described in the AltAz
frame. This is a local, angular coordinate
frame, with angles altitude
and azimuth
. Altitude is the
measured from the Horizon (0°) to the Zenith (90°). For the azimuth,
there are different conventions. In Astropy und thus ctapipe, Azimuth is
oriented East of North (i.e., N=0°, E=90°).
73 obstime = Time("2013-11-01T03:00")
74 location = EarthLocation.of_site("Roque de los Muchachos")
75
76 altaz = AltAz(location=location, obstime=obstime)
77
78 array_pointing = SkyCoord(
79 alt=event.pointing.array_azimuth,
80 az=event.pointing.array_altitude,
81 frame=altaz,
82 )
83
84 print(array_pointing)
<SkyCoord (AltAz: obstime=2013-11-01T03:00:00.000, location=(5327448.9957829, -1718665.73869569, 3051566.90295403) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
(69.99999967, 0.)>
CameraFrame#
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.
Typical usage: Position of pixels in the focal plane.
109 geometry = source.subarray.tel[tel_id].camera.geometry
110 pix_x = geometry.pix_x
111 pix_y = geometry.pix_y
112 focal_length = source.subarray.tel[tel_id].optics.equivalent_focal_length
115 telescope_pointing = SkyCoord(
116 alt=event.pointing.tel[tel_id].altitude,
117 az=event.pointing.tel[tel_id].azimuth,
118 frame=altaz,
119 )
120
121 camera_frame = CameraFrame(
122 focal_length=focal_length,
123 rotation=0 * u.deg,
124 telescope_pointing=telescope_pointing,
125 )
126
127 cam_coords = SkyCoord(x=pix_x, y=pix_y, frame=camera_frame)
128
129 print(cam_coords)
<SkyCoord (CameraFrame: focal_length=28.0 m, rotation=0.0 rad, telescope_pointing=<AltAz Coordinate (obstime=2013-11-01T03:00:00.000, location=(5327448.9957829, -1718665.73869569, 3051566.90295403) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
(0., 69.99999967)>, obstime=None, location=None): (x, y) in m
[( 0. , 0. ), (-0.0377967 , -0.03273243),
(-0.04724547, 0.01636666), ..., ( 0.67088033, -0.96561883),
(-0.45356484, -1.08017546), (-0.50081024, -1.06380875)]>
132 plt.scatter(cam_coords.x, cam_coords.y)
133 plt.title(f"Camera type: {geometry.name}")
134 plt.xlabel(f"x / {cam_coords.x.unit}")
135 plt.ylabel(f"y / {cam_coords.y.unit}")
136 plt.axis("square")
(-1.2888595582149416, 1.2888595582149416, -1.1881999735124904, 1.3895191429173928)
The implementation of the coordinate system with astropy makes it easier to use time of the observation and location of the observing site, to understand, for example which stars are visible during a certain night and how they might be visible in the camera.
147 location = EarthLocation.of_site("Roque de los Muchachos")
148 obstime = Time("2018-11-01T04:00")
149
150 crab = SkyCoord.from_name("crab nebula")
151
152 altaz = AltAz(location=location, obstime=obstime)
153
154 pointing = crab.transform_to(altaz)
155
156 camera_frame = CameraFrame(
157 telescope_pointing=pointing,
158 focal_length=focal_length,
159 obstime=obstime,
160 location=location,
161 )
162
163
164 subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
165 cam = subarray.tel[1].camera.geometry
166 fig, ax = plt.subplots()
167 display = CameraDisplay(cam, ax=ax)
168
169 ax.set_title(
170 f"La Palma, {obstime}, az={pointing.az.deg:.1f}°, zenith={pointing.zen.deg:.1f}°, camera={geometry.name}"
171 )
172
173 for i, name in enumerate(["crab nebula", "o tau", "zet tau"]):
174 star = SkyCoord.from_name(name)
175 star_cam = star.transform_to(camera_frame)
176
177 x = star_cam.x.to_value(u.m)
178 y = star_cam.y.to_value(u.m)
179
180 ax.plot(x, y, marker="*", color=f"C{i}")
181 ax.annotate(
182 name,
183 xy=(x, y),
184 xytext=(5, 5),
185 textcoords="offset points",
186 color=f"C{i}",
187 )
188
189 plt.show()
TelescopeFrame#
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 what astropy calls a
SkyOffsetFrame
.
The axis of the telescope frame, fov_lon
and fov_lat
, are
aligned with the horizontal system’s azimuth and altitude respectively.
Pointing corrections should applied to the transformation between this frame and the camera frame.
210 telescope_frame = TelescopeFrame(
211 telescope_pointing=pointing,
212 obstime=pointing.obstime,
213 location=pointing.location,
214 )
215 telescope_coords = cam_coords.transform_to(telescope_frame)
218 wrap_angle = telescope_pointing.az + 180 * u.deg
219
220 plt.axis("equal")
221 plt.scatter(
222 telescope_coords.fov_lon.deg, telescope_coords.fov_lat.deg, alpha=0.2, color="gray"
223 )
224
225
226 for i, name in enumerate(["crab nebula", "o tau", "zet tau"]):
227 star = SkyCoord.from_name(name)
228 star_tel = star.transform_to(telescope_frame)
229
230 plt.plot(star_tel.fov_lon.deg, star_tel.fov_lat.deg, "*", ms=10)
231 plt.annotate(
232 name,
233 xy=(star_tel.fov_lon.deg, star_tel.fov_lat.deg),
234 xytext=(5, 5),
235 textcoords="offset points",
236 color=f"C{i}",
237 )
238
239 plt.xlabel("fov_lon / {}".format(telescope_coords.altaz.az.unit))
240 plt.ylabel("fov_lat / {}".format(telescope_coords.altaz.alt.unit))
Text(101.97222222222221, 0.5, 'fov_lat / deg')
NominalFrame#
Nominal coordinate frame. A Frame using a
UnitSphericalRepresentation
. This is basically the same as a
HorizonCoordinate
, but the origin is at an arbitrary position in the
sky. This is what astropy calls a SkyOffsetFrame
If the telescopes
are in divergent pointing, this Frame
can be used to transform to a
common system. - 2D reconstruction (HillasIntersector
) is performed
in this frame - 3D reconstruction (HillasReconstructor
) doesn’t need
this frame
Let’s play a bit with 3 LSTs with divergent pointing
262 location = EarthLocation.of_site("Roque de los Muchachos")
263 obstime = Time("2018-11-01T02:00")
264 altaz = AltAz(location=location, obstime=obstime)
265
266 crab = SkyCoord.from_name("crab nebula")
267
268 # let's observe crab
269 array_pointing = crab.transform_to(altaz)
270
271
272 # let the telescopes point to different positions
273 alt_offsets = u.Quantity([1, -1, -1], u.deg)
274 az_offsets = u.Quantity([0, -2, +2], u.deg)
275
276
277 tel_pointings = SkyCoord(
278 alt=array_pointing.alt + alt_offsets,
279 az=array_pointing.az + az_offsets,
280 frame=altaz,
281 )
282
283 camera_frames = CameraFrame(
284 telescope_pointing=tel_pointings, # multiple pointings, so we get multiple frames
285 focal_length=focal_length,
286 obstime=obstime,
287 location=location,
288 )
289
290 nom_frame = NominalFrame(origin=array_pointing, obstime=obstime, location=location)
293 fig, ax = plt.subplots(figsize=(15, 10))
294 ax.set_aspect(1)
295
296 for i in range(3):
297 cam_coord = SkyCoord(x=pix_x, y=pix_y, frame=camera_frames[i])
298 nom_coord = cam_coord.transform_to(nom_frame)
299
300 ax.scatter(
301 x=nom_coord.fov_lon.deg,
302 y=nom_coord.fov_lat.deg,
303 label=f"Telescope {i + 1}",
304 s=30,
305 alpha=0.15,
306 )
307
308
309 for i, name in enumerate(["Crab", "o tau", "zet tau"]):
310 s = SkyCoord.from_name(name)
311 s_nom = s.transform_to(nom_frame)
312 ax.plot(
313 s_nom.fov_lon.deg,
314 s_nom.fov_lat.deg,
315 "*",
316 ms=10,
317 )
318 ax.annotate(
319 name,
320 xy=(s_nom.fov_lon.deg, s_nom.fov_lat.deg),
321 xytext=(5, 5),
322 textcoords="offset points",
323 color=f"C{i}",
324 )
325
326
327 ax.set_xlabel("fov_lon / deg")
328 ax.set_ylabel("fov_lat / deg")
329
330 ax.legend()
331 plt.show()
GroundFrame#
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
Typical usage: positions of telescopes on the ground (x, y, z)
347 source.subarray.peek()
<ctapipe.visualization.mpl_array.ArrayDisplay object at 0x7f746fc03850>
In case a layout is selected, the following line will produce a different output from the picture above.
<ctapipe.visualization.mpl_array.ArrayDisplay object at 0x7f746c147880>
In this image all the telescope from the gamma_test.simtel.gz
file
are plotted as spheres in the GroundFrame.
TiltedGroundFrame#
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. The plane is rotated along the z_axis by the
azimuth of the pointing_direction
and then it is inclined with an
angle equal to the zenith angle of the pointing_direction
.
This frame is used for the reconstruction of the shower core position.
This image picture both the telescopes in the GroundFrame (red) and in
the TiltedGroundFrame (green) are displayed: in this case since the
azimuth of the pointing_direction
is 0 degrees, then the plane is
just tilted according to the zenith angle.
For playing with these and with more 3D models of the telescopes themselves, have a look at the CREED_VTK library.
Total running time of the script: (0 minutes 18.158 seconds)