\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online [![Binder](https://mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.10?urlpath=lab/tree/image_analysis.ipynb)\n",
"- You can contribute with your own notebooks in this\n",
"[GitHub repository](https://github.com/gammapy/gammapy/tree/master/tutorials).\n",
"- **Source files:**\n",
"[image_analysis.ipynb](../_static/notebooks/image_analysis.ipynb) |\n",
"[image_analysis.py](../_static/notebooks/image_analysis.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fitting 2D images with Gammapy\n",
"\n",
"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: \n",
"\n",
"1. Using the sherpa pacakge, see: [image_fitting_with_sherpa.ipynb](image_fitting_with_sherpa.ipynb),\n",
"\n",
"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](analysis_3d.ipynb), taking care of a few things that we mention in this tutorial\n",
"\n",
"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.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"\n",
"from gammapy.data import DataStore\n",
"from gammapy.irf import make_mean_psf\n",
"from gammapy.maps import Map, MapAxis, WcsGeom\n",
"from gammapy.cube import MapMaker, PSFKernel, MapFit\n",
"from gammapy.cube.models import SkyModel\n",
"from gammapy.spectrum.models import PowerLaw2\n",
"from gammapy.image.models import SkyPointSource\n",
"\n",
"from regions import CircleSkyRegion"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Prepare modeling input data\n",
"\n",
"### The counts, exposure and the background maps\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Data store:\n",
"HDU index table:\n",
"BASE_DIR: /Users/jer/DATA/GAMMAPY/cta-1dc/index/gps\n",
"Rows: 24\n",
"OBS_ID: 110380 -- 111630\n",
"HDU_TYPE: ['aeff', 'bkg', 'edisp', 'events', 'gti', 'psf']\n",
"HDU_CLASS: ['aeff_2d', 'bkg_3d', 'edisp_2d', 'events', 'gti', 'psf_3gauss']\n",
"\n",
"Observation table:\n",
"Observatory name: 'N/A'\n",
"Number of observations: 4\n"
]
}
],
"source": [
"# Define which data to use and print some information\n",
"data_store = DataStore.from_dir(\"$GAMMAPY_DATA/cta-1dc/index/gps/\")\n",
"data_store.info()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Total observation time: 2.0 h\n"
]
}
],
"source": [
"print(\n",
" \"Total observation time: {}\".format(\n",
" data_store.obs_table[\"ONTIME\"].quantity.sum().to(\"hour\")\n",
" )\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Select some observations from these dataset by hand\n",
"obs_ids = [110380, 111140, 111159]\n",
"observations = data_store.get_observations(obs_ids)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"emin, emax = [0.1, 10] * u.TeV\n",
"energy_axis = MapAxis.from_bounds(\n",
" emin.value, emax.value, 10, unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
")\n",
"geom = WcsGeom.create(\n",
" skydir=(0, 0),\n",
" binsz=0.02,\n",
" width=(10, 8),\n",
" coordsys=\"GAL\",\n",
" proj=\"CAR\",\n",
" axes=[energy_axis],\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 11.6 s, sys: 2.31 s, total: 13.9 s\n",
"Wall time: 14.1 s\n"
]
}
],
"source": [
"%%time\n",
"maker = MapMaker(geom, offset_max=4.0 * u.deg)\n",
"maps = maker.run(observations)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"{'counts': WcsNDMap\n",
" \n",
" \tgeom : WcsGeom \n",
" \taxes : lon, lat, energy\n",
" \tshape : (500, 400, 10)\n",
" \tndim : 3\n",
" \tunit : '' \n",
" \tdtype : float32 , 'exposure': WcsNDMap\n",
" \n",
" \tgeom : WcsGeom \n",
" \taxes : lon, lat, energy\n",
" \tshape : (500, 400, 10)\n",
" \tndim : 3\n",
" \tunit : 'm2 s' \n",
" \tdtype : float32 , 'background': WcsNDMap\n",
" \n",
" \tgeom : WcsGeom \n",
" \taxes : lon, lat, energy\n",
" \tshape : (500, 400, 10)\n",
" \tndim : 3\n",
" \tunit : '' \n",
" \tdtype : float32 }"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"maps"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"spectrum = PowerLaw2(index=2)\n",
"maps2D = maker.make_images(spectrum=spectrum, keepdims=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"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."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# mean PSF\n",
"geom2d = maps2D[\"exposure\"].geom\n",
"src_pos = SkyCoord(0, 0, unit=\"deg\", frame=\"galactic\")\n",
"table_psf = make_mean_psf(observations, src_pos)\n",
"\n",
"table_psf_2d = table_psf.table_psf_in_energy_band(\n",
" (emin, emax), spectrum=spectrum\n",
")\n",
"\n",
"# PSF kernel used for the model convolution\n",
"psf_kernel = PSFKernel.from_table_psf(\n",
" table_psf_2d, geom2d, max_radius=\"0.3 deg\"\n",
")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, the analysis proceeds as usual. Just take care to use the proper geometry in this case."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Define a mask"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"mask = Map.from_geom(geom2d)\n",
"\n",
"region = CircleSkyRegion(center=src_pos, radius=0.6 * u.deg)\n",
"mask.data = mask.geom.region_mask([region])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modeling the source\n",
"\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"spatial_model = SkyPointSource(lon_0=\"0.01 deg\", lat_0=\"0.01 deg\")\n",
"spectral_model = PowerLaw2(\n",
" emin=emin, emax=emax, index=2.0, amplitude=\"3e-12 cm-2 s-1\"\n",
")\n",
"model = SkyModel(spatial_model=spatial_model, spectral_model=spectral_model)\n",
"model.parameters[\"index\"].frozen = True"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"fit = MapFit(\n",
" model=model,\n",
" counts=maps2D[\"counts\"],\n",
" exposure=maps2D[\"exposure\"],\n",
" background=maps2D[\"background\"],\n",
" mask=mask,\n",
" psf=psf_kernel,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 7.57 s, sys: 470 ms, total: 8.04 s\n",
"Wall time: 8.08 s\n"
]
}
],
"source": [
"%%time\n",
"result = fit.run()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To see the actual best-fit parameters, do a print on the result"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"SkyModel\n",
"\n",
"Parameters: \n",
"\n",
"\t name value error unit min max frozen\n",
"\t--------- ---------- ----- -------- --- --- ------\n",
"\t lon_0 -5.365e-02 nan deg nan nan False\n",
"\t lat_0 -5.059e-02 nan deg nan nan False\n",
"\tamplitude 4.264e-11 nan cm-2 s-1 nan nan False\n",
"\t index 2.000e+00 nan nan nan True\n",
"\t emin 1.000e-01 nan TeV nan nan True\n",
"\t emax 1.000e+01 nan TeV nan nan True\n"
]
}
],
"source": [
"print(model)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Todo: Demonstrate plotting a flux map"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercises\n",
"1. Plot residual maps as done in the `analysis_3d` notebook\n",
"2. Iteratively add and fit sources as explained in `image_fitting_with_sherpa` notebook\n"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.0"
},
"nbsphinx": {
"orphan": true
}
},
"nbformat": 4,
"nbformat_minor": 2
}