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

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

Fitting 2D images with Gammapy

Gammapy does not have any special handling for 2D images, but treats them as a subset of maps. Thus, classical 2D image analysis can be done in 2 independent ways:

1. Using the sherpa pacakge, see: image_fitting_with_sherpa.ipynb,
2. Within gammapy, by assuming 2D analysis to be a sub-set of the generalised maps. Thus, analysis should proceeexactly as demonstrated in analysis_3d.ipynb, taking care of a few things that we mention in this tutorial

We consider 2D images to be a special case of 3D maps, ie, maps with only one energy bin. This is a major difference while analysing in sherpa, where the maps must not contain any energy axis. In this tutorial, we do a classical image analysis using three example observations of the Galactic center region with CTA - i.e., study the source flux and morphology.

Setup

%matplotlib inline

import astropy.units as u
from astropy.coordinates import SkyCoord

from gammapy.data import DataStore
from gammapy.irf import make_mean_psf
from gammapy.maps import Map, MapAxis, WcsGeom
from gammapy.cube import MapMaker, PSFKernel, MapFit
from gammapy.cube.models import SkyModel
from gammapy.spectrum.models import PowerLaw2
from gammapy.image.models import SkyPointSource

from regions import CircleSkyRegion

Prepare modeling input data

The counts, exposure and the background maps

This is the same drill - use DataStore object to access the CTA observations and retrieve a list of observations by passing the observations IDs to the .get_observations() method, then use MapMaker to make the maps.

# Define which data to use and print some information
data_store = DataStore.from_dir("$GAMMAPY_DATA/cta-1dc/index/gps/")
data_store.info()

Data store:
HDU index table:
BASE_DIR: /Users/jer/DATA/GAMMAPY/cta-1dc/index/gps
Rows: 24
OBS_ID: 110380 -- 111630
HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf']
HDU_CLASS: ['aeff_2d', 'bkg_3d', 'edisp_2d', 'events', 'gti', 'psf_3gauss']

Observation table:
Observatory name: 'N/A'
Number of observations: 4

print(
    "Total observation time: {}".format(
        data_store.obs_table["ONTIME"].quantity.sum().to("hour")
    )
)

Total observation time: 2.0 h

# Select some observations from these dataset by hand
obs_ids = [110380, 111140, 111159]
observations = data_store.get_observations(obs_ids)

emin, emax = [0.1, 10] * u.TeV
energy_axis = MapAxis.from_bounds(
    emin.value, emax.value, 10, unit="TeV", name="energy", interp="log"
)
geom = WcsGeom.create(
    skydir=(0, 0),
    binsz=0.02,
    width=(10, 8),
    coordsys="GAL",
    proj="CAR",
    axes=[energy_axis],
)

Note that even when doing a 2D analysis, it is better to use fine energy bins in the beginning and then sum them over. This is to ensure that the background shape can be approximated by a power law function in each energy bin.

%%time
maker = MapMaker(geom, offset_max=4.0 * u.deg)
maps = maker.run(observations)

CPU times: user 11.6 s, sys: 2.31 s, total: 13.9 s
Wall time: 14.1 s

maps

{'counts': WcsNDMap

    geom  : WcsGeom
    axes  : lon, lat, energy
    shape : (500, 400, 10)
    ndim  : 3
    unit  : ''
    dtype : float32 , 'exposure': WcsNDMap

    geom  : WcsGeom
    axes  : lon, lat, energy
    shape : (500, 400, 10)
    ndim  : 3
    unit  : 'm2 s'
    dtype : float32 , 'background': WcsNDMap

    geom  : WcsGeom
    axes  : lon, lat, energy
    shape : (500, 400, 10)
    ndim  : 3
    unit  : ''
    dtype : float32 }

As we can see, the maps now have multiple bins in energy. We need to squash them to have one bin, and this can be done by specifying keep_dims = True while calling make_images(). This will compute a summed counts and background maps, and a spectral weighted exposure map.

spectrum = PowerLaw2(index=2)
maps2D = maker.make_images(spectrum=spectrum, keepdims=True)

For a typical 2D analysis, using an energy dispersion usually does not make sense. A PSF map can be made as in the regular 3D case, taking care to weight it properly with the spectrum.

# mean PSF
geom2d = maps2D["exposure"].geom
src_pos = SkyCoord(0, 0, unit="deg", frame="galactic")
table_psf = make_mean_psf(observations, src_pos)

table_psf_2d = table_psf.table_psf_in_energy_band(
    (emin, emax), spectrum=spectrum
)

# PSF kernel used for the model convolution
psf_kernel = PSFKernel.from_table_psf(
    table_psf_2d, geom2d, max_radius="0.3 deg"
)

Now, the analysis proceeds as usual. Just take care to use the proper geometry in this case.

Define a mask

mask = Map.from_geom(geom2d)

region = CircleSkyRegion(center=src_pos, radius=0.6 * u.deg)
mask.data = mask.geom.region_mask([region])

Modeling the source

This is the important thing to note in this analysis. Since modelling and fitting in gammapy.maps needs to have a combination of spectral models, we have to use a dummy Powerlaw as for the spectral model and fix its index to 2. Since we are interested only in the integral flux, we will use the PowerLaw2 model which directly fits an integral flux.

spatial_model = SkyPointSource(lon_0="0.01 deg", lat_0="0.01 deg")
spectral_model = PowerLaw2(
    emin=emin, emax=emax, index=2.0, amplitude="3e-12 cm-2 s-1"
)
model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)
model.parameters["index"].frozen = True

fit = MapFit(
    model=model,
    counts=maps2D["counts"],
    exposure=maps2D["exposure"],
    background=maps2D["background"],
    mask=mask,
    psf=psf_kernel,
)

%%time
result = fit.run()

CPU times: user 7.57 s, sys: 470 ms, total: 8.04 s
Wall time: 8.08 s

To see the actual best-fit parameters, do a print on the result

print(model)

SkyModel

Parameters:

           name     value    error   unit   min max frozen
        --------- ---------- ----- -------- --- --- ------
            lon_0 -5.365e-02   nan      deg nan nan  False
            lat_0 -5.059e-02   nan      deg nan nan  False
        amplitude  4.264e-11   nan cm-2 s-1 nan nan  False
            index  2.000e+00   nan          nan nan   True
             emin  1.000e-01   nan      TeV nan nan   True
             emax  1.000e+01   nan      TeV nan nan   True

Todo: Demonstrate plotting a flux map

Exercises

1. Plot residual maps as done in the analysis_3d notebook
2. Iteratively add and fit sources as explained in image_fitting_with_sherpa notebook