Displaying Camera Images¶
[1]:
import astropy.coordinates as c
import astropy.units as u
import matplotlib.pylab as plt
import numpy as np
from ctapipe.coordinates import CameraFrame, EngineeringCameraFrame, TelescopeFrame
from ctapipe.image import hillas_parameters, tailcuts_clean, toymodel
from ctapipe.instrument import SubarrayDescription
from ctapipe.visualization import CameraDisplay
First, let’s create a fake Cherenkov image from a given CameraGeometry
and fill it with some data that we can draw later.
[2]:
# load an example camera geometry from a simulation file
subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
geom = subarray.tel[100].camera.geometry
# 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,
psi="35d",
)
image, sig, bg = model.generate_image(geom, intensity=1500, nsb_level_pe=10)
mask = tailcuts_clean(geom, image, picture_thresh=15, boundary_thresh=5)
[3]:
geom
[3]:
CameraGeometry(name='NectarCam', pix_type=PixelShape.HEXAGON, npix=1855, cam_rot=0.000 deg, pix_rot=40.893 deg, frame=<CameraFrame Frame (focal_length=16.445049285888672 m, rotation=0.0 rad, telescope_pointing=None, obstime=None, location=None)>)
Displaying Images¶
The simplest plot is just to generate a CameraDisplay with an image in its constructor. A figure and axis will be created automatically
[4]:
CameraDisplay(geom)
[4]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7ff6bc25daf0>
You can also specify the initial image
, cmap
and norm
(colomap and normalization, see below), title
to use. You can specify ax
if you want to draw the camera on an existing matplotlib Axes
object (otherwise one is created).
To change other options, or to change options dynamically, you can call the relevant functions of the CameraDisplay
object that is returned. For example to add a color bar, call add_colorbar()
, or to change the color scale, modify the cmap
or norm
properties directly.
Choosing a coordinate frame¶
The CameraGeometry
object contains a ctapipe.coordinates.Frame
used by CameraDisplay
to draw the camera in the correct orientation and distance units. The default frame is the CameraFrame
, which will display the camera in units of meters and with an orientation that the top of the camera (when parked) is aligned to the X-axis. To show the camera in another orientation, it’s useful to apply a coordinate transform to the CameraGeometry
before passing it to the
CameraDisplay
. The following Frames
are supported: * EngineeringCameraFrame
: similar to CameraFrame, but with the top of the camera aligned to the Y axis * TelescopeFrame
: In degrees (on the sky) coordinates relative to the telescope Alt/Az pointing position, with the Alt axis pointing upward.
[5]:
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
CameraDisplay(geom, image=image, ax=ax[0])
CameraDisplay(geom.transform_to(EngineeringCameraFrame()), image=image, ax=ax[1])
CameraDisplay(geom.transform_to(TelescopeFrame()), image=image, ax=ax[2])
[5]:
<ctapipe.visualization.mpl_camera.CameraDisplay at 0x7ff679a1c610>
Note the the name of the Frame appears in the lower-right corner
For the rest of this demo, let’s use the TelescopeFrame
[6]:
geom_camframe = geom
geom = geom_camframe.transform_to(EngineeringCameraFrame())
Changing the color map and scale¶
CameraDisplay supports any matplotlib color map It is highly recommended to use a perceptually uniform map, unless you have a good reason not to.
[7]:
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
for ii, cmap in enumerate(["PuOr_r", "rainbow", "twilight"]):
disp = CameraDisplay(geom, image=image, ax=ax[ii], title=cmap)
disp.add_colorbar()
disp.cmap = cmap
By default the minimum and maximum of the color bar are set automatically by the data in the image. To choose fixed limits, use:`
[8]:
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
for ii, minmax in enumerate([(10, 50), (-10, 10), (1, 100)]):
disp = CameraDisplay(geom, image=image, ax=ax[ii], title=minmax)
disp.add_colorbar()
disp.set_limits_minmax(minmax[0], minmax[1])
Or you can set the maximum limit by percentile of the charge distribution:
[9]:
fig, ax = plt.subplots(1, 3, figsize=(15, 4))
for ii, pct in enumerate([30, 50, 90]):
disp = CameraDisplay(geom, image=image, ax=ax[ii], title=f"{pct} %")
disp.add_colorbar()
disp.set_limits_percent(pct)
Using different normalizations¶
You can choose from several preset normalizations (lin, log, symlog) and also provide a custom normalization, for example a PowerNorm
:
[10]:
from matplotlib.colors import PowerNorm
fig, axes = plt.subplots(2, 2, figsize=(14, 10))
norms = ["lin", "log", "symlog", PowerNorm(0.5)]
for norm, ax in zip(norms, axes.flatten()):
disp = CameraDisplay(geom, image=image, ax=ax)
disp.norm = norm
disp.add_colorbar()
ax.set_title(str(norm))
axes[1, 1].set_title("PowerNorm(0.5)")
plt.show()
Overlays¶
Marking pixels¶
here we will mark pixels in the image mask. That will change their outline color
[11]:
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
disp = CameraDisplay(
geom, image=image, cmap="gray", ax=ax[0], title="Image mask in green"
)
disp.highlight_pixels(mask, alpha=0.8, linewidth=2, color="green")
disp = CameraDisplay(
geom, image=image, cmap="gray", ax=ax[1], title="Image mask in green (zoom)"
)
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)
[11]:
(-0.5, 0.5)
Drawing a Hillas-parameter ellipse¶
For this, we will first compute some Hillas Parameters in the current frame:
[12]:
clean_image = image.copy()
clean_image[~mask] = 0
hillas = hillas_parameters(geom, clean_image)
plt.figure(figsize=(6, 6))
disp = CameraDisplay(geom, image=image, cmap="gray_r")
disp.highlight_pixels(mask, alpha=0.5, color="dodgerblue")
disp.overlay_moments(hillas, color="red", linewidth=3, with_label=False)
Drawing a marker at a coordinate¶
This depends on the coordinate frame of the CameraGeometry
. Here we will sepcify the coordinate the EngineerngCameraFrame
, 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 automatically (if sufficient frame attributes are set).
Note that the parameter keep_old
is False by default, meaning adding a new point will clear the previous ones (useful for animations, but perhaps unexpected for a static plot). Set it to True
to plot multiple markers.
[13]:
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())
disp.overlay_coordinate(coord, markersize=20, marker="*")
disp.overlay_coordinate(
coord_in_another_frame, markersize=20, marker="*", keep_old=True
)
Generating an animation¶
Here we will make an animation of fake events by re-using a single display (much faster than generating a new one each time)
[14]:
from IPython import display
from matplotlib.animation import FuncAnimation
subarray = SubarrayDescription.read("dataset://gamma_prod5.simtel.zst")
geom = subarray.tel[1].camera.geometry
fov = 1.0
maxwid = 0.05
maxlen = 0.1
fig, ax = plt.subplots(1, 1, figsize=(8, 6))
disp = CameraDisplay(geom, ax=ax) # we only need one display (it can be re-used)
disp.cmap = "inferno"
disp.add_colorbar(ax=ax)
def update(frame):
"""this function will be called for each frame of the animation"""
x, y = np.random.uniform(-fov, fov, size=2)
width = np.random.uniform(0.01, maxwid)
length = np.random.uniform(width, maxlen)
angle = np.random.uniform(0, 180)
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,
psi=angle * u.deg,
)
image, _, _ = model.generate_image(
geom,
intensity=intens,
nsb_level_pe=5,
)
disp.image = image
# Create the animation and convert to a displayable video:
anim = FuncAnimation(fig, func=update, frames=10, interval=200)
plt.close(fig) # so it doesn't display here
video = anim.to_html5_video()
display.display(display.HTML(video))
Using CameraDisplays interactively¶
CameraDisplays
can be used interactivly whe displayed in a window, and also when using Jupyter notebooks/lab with appropriate backends.
When this is the case, the same CameraDisplay
object can be re-used. We can’t show this here in the documentation, but creating an animation when in a matplotlib window is quite easy! Try this in an interactive ipython session:
Running interactive displays in a matplotlib window¶
ipython -i --maplotlib=auto
That will open an ipython session with matplotlib graphics in a separate thread, meaning that you can type code and interact with plots simultaneneously.
In the ipython session try running the following code and you will see an animation (here in the documentation, it will of course be static)
First we load some real data so we have a nice image to view:
[15]:
import matplotlib.pyplot as plt
from ctapipe.io import EventSource
from ctapipe.visualization import CameraDisplay
import numpy as np
DATA = "dataset://gamma_20deg_0deg_run1___cta-prod5-lapalma_desert-2158m-LaPalma-dark_100evts.simtel.zst"
with EventSource(
DATA,
max_events=1,
focal_length_choice="EQUIVALENT",
) as source:
event = next(iter(source))
tel_id = list(event.r0.tel.keys())[0]
geom = source.subarray.tel[tel_id].camera.geometry
waveform = event.r0.tel[tel_id].waveform
n_chan, n_pix, n_samp = waveform.shape
Running the following the will bring up a window and animate the shower image as a function of time.
[16]:
disp = CameraDisplay(geom)
for ii in range(n_samp):
disp.image = waveform[0, :, ii]
plt.pause(0.1) # this lets matplotlib re-draw the scene
The output will be similar to the static animation created as follows:
[17]:
fig, ax = plt.subplots(1, 1)
disp = CameraDisplay(geom, ax=ax)
disp.add_colorbar()
disp.autoscale = False
def draw_sample(frame):
ax.set_title(f"sample: {frame}")
disp.set_limits_minmax(200, 400)
disp.image = waveform[0, :, frame]
anim = FuncAnimation(fig, func=draw_sample, frames=n_samp, interval=100)
plt.close(fig) # so it doesn't display here
video = anim.to_html5_video()
display.display(display.HTML(video))
Making it clickable¶
Also when running in a window, you can enable the disp.enable_pixel_picker()
option. This will then allow the user to click a pixel and a function will run. By default the function simply prints the pixel and value to stdout, however you can override the function on_pixel_clicked(pix_id)
to do anything you want by making a subclass
[18]:
class MyCameraDisplay(CameraDisplay):
def on_pixel_clicked(self, pix_id):
print(f"{pix_id=} has value {self.image[pix_id]:.2f}")
[19]:
disp = MyCameraDisplay(geom, image=image)
disp.enable_pixel_picker()
then, when a user clicks a pixel it would print:
pixel 5 has value 2.44
[ ]: