Skip to content

Commit

Permalink
notebooks processing in TelescopeFrame
Browse files Browse the repository at this point in the history
  • Loading branch information
StFroese authored and maxnoe committed Sep 22, 2023
1 parent 00627eb commit c1be72f
Show file tree
Hide file tree
Showing 5 changed files with 55 additions and 52 deletions.
11 changes: 6 additions & 5 deletions examples/algorithms/dilate_image.py
Original file line number Diff line number Diff line change
Expand Up @@ -13,20 +13,21 @@
# %matplotlib inline
from matplotlib import pyplot as plt

from ctapipe.coordinates import TelescopeFrame
from ctapipe.image import dilate, tailcuts_clean, toymodel
from ctapipe.instrument import SubarrayDescription
from ctapipe.visualization import CameraDisplay

# Load a camera from an example file
subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
geom = subarray.tel[100].camera.geometry
geom = subarray.tel[100].camera.geometry.transform_to(TelescopeFrame())

# Create a fake camera image to display:
model = toymodel.Gaussian(
x=0.2 * u.m,
y=0.0 * u.m,
width=0.05 * u.m,
length=0.15 * u.m,
x=0.2 * u.deg,
y=0.0 * u.deg,
width=0.05 * u.deg,
length=0.15 * u.deg,
psi="35d",
)

Expand Down
13 changes: 7 additions & 6 deletions examples/tutorials/calibrated_data_exploration.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@

import ctapipe
from ctapipe.calib import CameraCalibrator
from ctapipe.coordinates import TelescopeFrame
from ctapipe.image import hillas_parameters, tailcuts_clean
from ctapipe.io import EventSource
from ctapipe.utils.datasets import get_dataset_path
Expand Down Expand Up @@ -98,7 +99,7 @@

tel_id = sorted(event.r1.tel.keys())[1]
sub = source.subarray
geometry = sub.tel[tel_id].camera.geometry
geometry = sub.tel[tel_id].camera.geometry.transform_to(TelescopeFrame())
image = event.dl1.tel[tel_id].image

######################################################################
Expand Down Expand Up @@ -130,9 +131,8 @@
disp.overlay_moments(params, color="red", lw=3)
disp.highlight_pixels(mask, color="white", alpha=0.3, linewidth=2)

plt.xlim(params.x.to_value(u.m) - 0.5, params.x.to_value(u.m) + 0.5)
plt.ylim(params.y.to_value(u.m) - 0.5, params.y.to_value(u.m) + 0.5)

plt.xlim(params.fov_lon.to_value(u.deg) - 2, params.fov_lon.to_value(u.deg) + 2)
plt.ylim(params.fov_lat.to_value(u.deg) - 2, params.fov_lat.to_value(u.deg) + 2)
######################################################################
source.metadata

Expand Down Expand Up @@ -166,6 +166,7 @@
cams_in_event = tels_in_event.intersection(cam_ids)
first_tel_id = list(cams_in_event)[0]
tel = sub.tel[first_tel_id]
tel_geom = tel.camera.geometry.transform_to(TelescopeFrame())
print("{}s in event: {}".format(tel, cams_in_event))


Expand All @@ -174,7 +175,7 @@
#

image_sum = np.zeros_like(
tel.camera.geometry.pix_x.value
tel_geom.pix_x.value
) # just make an array of 0's in the same shape as the camera

for tel_id in cams_in_event:
Expand All @@ -187,7 +188,7 @@

plt.figure(figsize=(8, 8))

disp = CameraDisplay(tel.camera.geometry, image=image_sum)
disp = CameraDisplay(tel_geom, image=image_sum)
disp.overlay_moments(params, with_label=False)
plt.title("Sum of {}x {}".format(len(cams_in_event), tel))

Expand Down
42 changes: 23 additions & 19 deletions examples/tutorials/ctapipe_handson.py
Original file line number Diff line number Diff line change
Expand Up @@ -22,6 +22,7 @@

from ctapipe import utils
from ctapipe.calib import CameraCalibrator
from ctapipe.coordinates import TelescopeFrame
from ctapipe.image import hillas_parameters, tailcuts_clean
from ctapipe.io import EventSource, HDF5TableWriter
from ctapipe.visualization import CameraDisplay
Expand Down Expand Up @@ -129,20 +130,21 @@ def view_waveform(chan=0, pix_id=brightest_pixel):
tel.optics

######################################################################
tel.camera.geometry.pix_x
geom = tel.camera.geometry.transform_to(TelescopeFrame())
geom.pix_x

######################################################################
tel.camera.geometry.to_table()
geom.to_table()

######################################################################
tel.optics.mirror_area


######################################################################
disp = CameraDisplay(tel.camera.geometry)
disp = CameraDisplay(geom)

######################################################################
disp = CameraDisplay(tel.camera.geometry)
disp = CameraDisplay(geom)
disp.image = r0tel.waveform[
0, :, 10
] # display channel 0, sample 0 (try others like 10)
Expand Down Expand Up @@ -180,10 +182,10 @@ def view_waveform(chan=0, pix_id=brightest_pixel):
dl1tel.peak_time

######################################################################
CameraDisplay(tel.camera.geometry, image=dl1tel.image)
CameraDisplay(geom, image=dl1tel.image)

######################################################################
CameraDisplay(tel.camera.geometry, image=dl1tel.peak_time)
CameraDisplay(geom, image=dl1tel.peak_time)


######################################################################
Expand All @@ -192,33 +194,33 @@ def view_waveform(chan=0, pix_id=brightest_pixel):


image = dl1tel.image
mask = tailcuts_clean(tel.camera.geometry, image, picture_thresh=10, boundary_thresh=5)
mask = tailcuts_clean(geom, image, picture_thresh=10, boundary_thresh=5)
mask

######################################################################
CameraDisplay(tel.camera.geometry, image=mask)
CameraDisplay(geom, image=mask)

######################################################################
cleaned = image.copy()
cleaned[~mask] = 0

######################################################################
disp = CameraDisplay(tel.camera.geometry, image=cleaned)
disp = CameraDisplay(geom, image=cleaned)
disp.cmap = plt.cm.coolwarm
disp.add_colorbar()
plt.xlim(0.5, 1.0)
plt.ylim(-1.0, 0.0)
plt.xlim(-2.0, -1.0)
plt.ylim(1.0, 2.0)

######################################################################
params = hillas_parameters(tel.camera.geometry, cleaned)
params = hillas_parameters(geom, cleaned)
print(params)

######################################################################
disp = CameraDisplay(tel.camera.geometry, image=cleaned)
disp = CameraDisplay(geom, image=cleaned)
disp.cmap = plt.cm.coolwarm
disp.add_colorbar()
plt.xlim(0.5, 1.0)
plt.ylim(-1.0, 0.0)
plt.xlim(-2.0, -1.0)
plt.ylim(1.0, 2.0)
disp.overlay_moments(params, color="white", lw=2)


Expand Down Expand Up @@ -258,9 +260,10 @@ def view_waveform(chan=0, pix_id=brightest_pixel):

for tel_id, tel_data in event.dl1.tel.items():
tel = source.subarray.tel[tel_id]
mask = tailcuts_clean(tel.camera.geometry, tel_data.image)
geom = tel.camera.geometry.transform_to(TelescopeFrame())
mask = tailcuts_clean(geom, tel_data.image)
if np.count_nonzero(mask) > 0:
params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])
params = hillas_parameters(geom[mask], tel_data.image[mask])


######################################################################
Expand All @@ -272,8 +275,9 @@ def view_waveform(chan=0, pix_id=brightest_pixel):

for tel_id, tel_data in event.dl1.tel.items():
tel = source.subarray.tel[tel_id]
mask = tailcuts_clean(tel.camera.geometry, tel_data.image)
params = hillas_parameters(tel.camera.geometry[mask], tel_data.image[mask])
geom = tel.camera.geometry.transform_to(TelescopeFrame())
mask = tailcuts_clean(geom, tel_data.image)
params = hillas_parameters(geom[mask], tel_data.image[mask])
writer.write("hillas", params)


Expand Down
6 changes: 4 additions & 2 deletions examples/tutorials/ctapipe_overview.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,6 +25,7 @@
from traitlets.config import Config

from ctapipe.calib import CameraCalibrator
from ctapipe.coordinates import TelescopeFrame
from ctapipe.image import (
ImageProcessor,
camera_to_shower_coordinates,
Expand Down Expand Up @@ -203,7 +204,7 @@
tel_id = 130

######################################################################
geometry = source.subarray.tel[tel_id].camera.geometry
geometry = source.subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame())
dl1 = event.dl1.tel[tel_id]

geometry, dl1
Expand Down Expand Up @@ -294,7 +295,7 @@

######################################################################
long, trans = camera_to_shower_coordinates(
geometry.pix_x, geometry.pix_y, hillas.x, hillas.y, hillas.psi
geometry.pix_x, geometry.pix_y, hillas.fov_lon, hillas.fov_lat, hillas.psi
)

plt.plot(long[clean], dl1.peak_time[clean], "o")
Expand Down Expand Up @@ -540,6 +541,7 @@
)
image += new_image

geometry = geometry.transform_to(TelescopeFrame())
######################################################################
clean = tailcuts_clean(
geometry,
Expand Down
35 changes: 15 additions & 20 deletions examples/visualization/camera_display.py
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
from matplotlib.animation import FuncAnimation
from matplotlib.colors import PowerNorm

from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame, TelescopeFrame
from ctapipe.coordinates import EngineeringCameraFrame, TelescopeFrame
from ctapipe.image import hillas_parameters, tailcuts_clean, toymodel
from ctapipe.instrument import SubarrayDescription
from ctapipe.io import EventSource
Expand Down Expand Up @@ -100,7 +100,7 @@
#

geom_camframe = geom
geom = geom_camframe.transform_to(EngineeringCameraFrame())
geom = geom_camframe.transform_to(TelescopeFrame())


######################################################################
Expand Down Expand Up @@ -197,9 +197,8 @@
)
disp.highlight_pixels(mask, alpha=1, linewidth=3, color="green")

ax[1].set_ylim(-0.5, 0.5)
ax[1].set_xlim(-0.5, 0.5)

ax[1].set_ylim(-2, 2)
ax[1].set_xlim(-2, 2)

######################################################################
# Drawing a Hillas-parameter ellipse
Expand Down Expand Up @@ -227,7 +226,7 @@
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#
# This depends on the coordinate frame of the ``CameraGeometry``. Here we
# will sepcify the coordinate the ``EngineerngCameraFrame``, but if you
# will sepcify the coordinate in the ``CameraFrame``, but if you
# have enough information to do the coordinate transform, you could use
# ``ICRS`` coordinates and overlay star positions. ``CameraDisplay`` will
# convert the coordinate you pass in to the ``Frame`` of the display
Expand All @@ -242,12 +241,8 @@
plt.figure(figsize=(6, 6))
disp = CameraDisplay(geom, image=image, cmap="gray_r")

coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom.frame)
coord_in_another_frame = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=CameraFrame())
coord = c.SkyCoord(x=0.5 * u.m, y=0.7 * u.m, frame=geom_camframe.frame)
disp.overlay_coordinate(coord, markersize=20, marker="*")
disp.overlay_coordinate(
coord_in_another_frame, markersize=20, marker="*", keep_old=True
)


######################################################################
Expand All @@ -260,11 +255,11 @@


subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
geom = subarray.tel[1].camera.geometry
geom = subarray.tel[1].camera.geometry.transform_to(TelescopeFrame())

fov = 1.0
maxwid = 0.05
maxlen = 0.1
fov = 2.0
maxwid = 0.1
maxlen = 0.5

fig, ax = plt.subplots(1, 1, figsize=(8, 6))
disp = CameraDisplay(geom, ax=ax) # we only need one display (it can be re-used)
Expand All @@ -281,10 +276,10 @@ def update(frame):
intens = width * length * (5e4 + 1e5 * np.random.exponential(2))

model = toymodel.Gaussian(
x=x * u.m,
y=y * u.m,
width=width * u.m,
length=length * u.m,
x=x * u.deg,
y=y * u.deg,
width=width * u.deg,
length=length * u.deg,
psi=angle * u.deg,
)
image, _, _ = model.generate_image(
Expand Down Expand Up @@ -343,7 +338,7 @@ def update(frame):
event = next(iter(source))

tel_id = list(event.r0.tel.keys())[0]
geom = source.subarray.tel[tel_id].camera.geometry
geom = source.subarray.tel[tel_id].camera.geometry.transform_to(TelescopeFrame())
waveform = event.r0.tel[tel_id].waveform
n_chan, n_pix, n_samp = waveform.shape

Expand Down

0 comments on commit c1be72f

Please sign in to comment.