This is a fixed-text formatted version of a Jupyter notebook.

Try online You can also contribute with your own notebooks in this GitHub repository. Source files: image_analysis.ipynb | image_analysis.py

# Image analysis with Gammapy¶

## Introduction¶

This tutorial shows how to make a significance image of the Crab nebula with Gammapy.

TODO: Refactor gammapy.scripts.image_fit into a class (simiar to gammapy.spectrum.SpectrumFit) and run it here to fit a Gauss and get the position / extension.

## Setup¶

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

In [2]:

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.convolution import Ring2DKernel, Tophat2DKernel
from astropy.visualization import simple_norm

from gammapy.data import DataStore
from gammapy.image import SkyImage, SkyImageList
from gammapy.detect import KernelBackgroundEstimator as KBE


## Dataset¶

We will use the gammapy.data.DataStore to access some example data.

These are observations of the Crab nebula with H.E.S.S. (preliminary, events are simulated for now).

In [3]:

data_store = DataStore.from_dir('\$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/')


## Counts image¶

Let’s make a counts image using the SkyMap class.

In [4]:

source_pos = SkyCoord(83.633083, 22.0145, unit='deg')
# If you have internet access, you could also use this to define the source_pos:
# source_pos = SkyCoord.from_name('crab')
print(source_pos)

<SkyCoord (ICRS): (ra, dec) in deg
( 83.633083,  22.0145)>

In [5]:

ref_image = SkyImage.empty(
nxpix=400, nypix=400, binsz=0.02,
xref=source_pos.ra.deg, yref=source_pos.dec.deg,
coordsys='CEL', proj='TAN',
)

In [6]:

# Make a counts image for a single observation
events = data_store.obs(obs_id=23523).events
counts_image = SkyImage.empty_like(ref_image)
counts_image.fill_events(events)

In [7]:

norm = simple_norm(counts_image.data, stretch='sqrt', min_cut=0, max_cut=0.3)

Out[7]:

(<matplotlib.figure.Figure at 0x107b99b70>,
<matplotlib.axes._subplots.WCSAxesSubplot at 0x108533d68>,
<matplotlib.colorbar.Colorbar at 0x102666b00>)

In [8]:

# Making a counts image for multiple observations is a bit inconvenient at the moment
# we'll make that better soon.
# For now, you can do it like this:
obs_ids = [23523, 23526]
counts_image2 = SkyImage.empty_like(ref_image)
for obs_id in obs_ids:
events = data_store.obs(obs_id=obs_id).events
counts_image2.fill_events(events)

In [9]:

norm = simple_norm(counts_image2.data, stretch='sqrt', min_cut=0, max_cut=0.5)

Out[9]:

(<matplotlib.figure.Figure at 0x107d9c0b8>,
<matplotlib.axes._subplots.WCSAxesSubplot at 0x108f72710>,
<matplotlib.colorbar.Colorbar at 0x1093fac50>)


# Background modeling¶

In Gammapy a few different methods to estimate the background are available.

Here we’ll use the gammapy.detect.KernelBackgroundEstimator to make a background image and the make a significance image.

In [10]:

source_kernel = Tophat2DKernel(radius=5)
source_kernel.normalize(mode='peak')
source_kernel = source_kernel.array

background_kernel.normalize(mode='peak')
background_kernel = background_kernel.array

In [11]:

plt.imshow(source_kernel, interpolation='nearest', cmap='gray')
plt.colorbar()
plt.grid('off')

In [12]:

plt.imshow(background_kernel, interpolation='nearest', cmap='gray')
plt.colorbar()
plt.grid('off')

In [13]:

# To use the KernelBackgroundEstimator you first have to set
# up a source and background kernel and put the counts image input
# into a container SkyImageList class.
images = SkyImageList()
images['counts'] = counts_image2

kbe = KBE(
kernel_src=source_kernel,
kernel_bkg=background_kernel,
significance_threshold=5,
)
# This takes about 10 seconds on my machine
result = kbe.run(images)

/Users/deil/code/gammapy/gammapy/stats/poisson.py:237: RuntimeWarning: invalid value encountered in greater_equal
/Users/deil/code/gammapy/gammapy/stats/poisson.py:254: RuntimeWarning: invalid value encountered in sqrt
term_b = sqrt(n_on * log(n_on / mu_bkg) - n_on + mu_bkg)
/Users/deil/code/gammapy/gammapy/detect/kernel.py:188: RuntimeWarning: invalid value encountered in less
mask = (significance.data < p['significance_threshold']) | np.isnan(significance)

In [14]:

# Let's have a look at the background image and the exclusion mask

# This doesn't work yet ... need to do SkyImage.plot fixes:
# fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(10, 3))
# background_image.plot(ax=axes[0])
# exclusion_image.plot(ax=axes[1])
# significance_image.plot(ax=axes[2])

In [15]:

background_image = result['background']
norm = simple_norm(background_image.data, stretch='sqrt', min_cut=0, max_cut=0.5)

Out[15]:

(<matplotlib.figure.Figure at 0x1096246a0>,
<matplotlib.axes._subplots.WCSAxesSubplot at 0x10963b470>,
<matplotlib.colorbar.Colorbar at 0x108d482b0>)

In [16]:

result['exclusion'].plot()

Out[16]:

(<matplotlib.figure.Figure at 0x109522a20>,
<matplotlib.axes._subplots.WCSAxesSubplot at 0x108d70780>,
None)

In [17]:

significance_image = result['significance']

Out[17]:

(<matplotlib.figure.Figure at 0x108dd2f98>,
<matplotlib.axes._subplots.WCSAxesSubplot at 0x109409518>,
<matplotlib.colorbar.Colorbar at 0x108d73b70>)


## Morphology fit¶

TODO

In [18]:

# from gammapy.scripts import image_fit


## Exercises¶

• Compute the counts, excess, background and significance at the Crab nebula position.
• Make an energy distribution of the events at the Crab nebula position.

## What next?¶

TODO: summarise

TODO: give links what to do next.