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: first_steps.ipynb | first_steps.py

Getting started with Gammapy¶

Introduction¶

This is a getting started tutorial for Gammapy.

In this tutorial we will use the Second Fermi-LAT Catalog of High-Energy Sources (3FHL) catalog, corresponding event list and images to learn how to work with some of the central Gammapy data structures.

We will cover the following topics:

Sky maps
- We will learn how to handle image based data with gammapy using a Fermi-LAT 3FHL example image. We will work with the following classes:
Event lists
- We will learn how to handle event lists with Gammapy. Important for this are the following classes:
  - gammapy.data.EventList
  - astropy.table.Table
Source catalogs
- We will show how to load source catalogs with Gammapy and explore the data using the following classes:
  - gammapy.catalog.SourceCatalog, specifically gammapy.catalog.SourceCatalog3FHL
  - astropy.table.Table
Spectral models and flux points
- We will pick an example source and show how to plot its spectral model and flux points. For this we will use the following classes:
  - gammapy.spectrum.SpectralModel, specifically the PowerLaw2 model.
  - gammapy.spectrum.FluxPoints
  - astropy.table.Table

Setup¶

Important: to run this tutorial the environment variable GAMMAPY_DATA must be defined and point to the directory on your machine where the datasets needed are placed. To check whether your setup is correct you can execute the following cell:

[1]:

import os

path = os.path.expandvars("$GAMMAPY_DATA")

if not os.path.exists(path):
    raise Exception("gammapy-data repository not found!")
else:
    print("Great your setup is correct!")

Great your setup is correct!

In case you encounter an error, you can un-comment and execute the following cell to continue. But we recommend to set up your enviroment correctly as decribed here after you are done with this notebook.

[2]:

# os.environ['GAMMAPY_DATA'] = os.path.join(os.getcwd(), '..')

Now we can continue with the usual IPython notebooks and Python imports:

[3]:

%matplotlib inline
import matplotlib.pyplot as plt

[4]:

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

Maps¶

The gammapy.maps package contains classes to work with sky images and cubes.

In this section, we will use a simple 2D sky image and will learn how to:

Read sky images from FITS files
Smooth images
Plot images
Cutout parts from images
Reproject images to different WCS

[5]:

from gammapy.maps import Map

gc_3fhl = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz")

The image is a WCSNDMap object:

[6]:

gc_3fhl

[6]:

WcsNDMap

    geom  : WcsGeom
    axes  : lon, lat
    shape : (400, 200)
    ndim  : 2
    unit  : ''
    dtype : >i8

The shape of the image is 400 x 200 pixel and it is defined using a cartesian projection in galactic coordinates.

The geom attribute is a WcsGeom object:

[7]:

gc_3fhl.geom

[7]:

WcsGeom

    axes       : lon, lat
    shape      : (400, 200)
    ndim       : 2
    coordsys   : GAL
    projection : CAR
    center     : 0.0 deg, 0.0 deg
    width      : 20.0 deg x 10.0 deg

Let’s take a closer look a the .data attribute:

[8]:

gc_3fhl.data

[8]:

array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ...,
       [0, 0, 0, ..., 0, 0, 1],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 1]])

That looks familiar! It just an ordinary 2 dimensional numpy array, which means you can apply any known numpy method to it:

[9]:

print("Total number of counts in the image: {:.0f}".format(gc_3fhl.data.sum()))

Total number of counts in the image: 32684

To show the image on the screen we can use the plot method. It basically calls plt.imshow, passing the gc_3fhl.data attribute but in addition handles axis with world coordinates using wcsaxes and defines some defaults for nicer plots (e.g. the colormap ‘afmhot’):

[10]:

gc_3fhl.plot(stretch="sqrt");

../_images/notebooks_first_steps_20_0.png

To make the structures in the image more visible we will smooth the data using a Gausian kernel with a radius of 0.5 deg. Again smooth() is a wrapper around existing functionality from the scientific Python libraries. In this case it is Scipy’s gaussian_filter method. For convenience the kernel shape can be specified with as string and the smoothing radius with a quantity. It returns again a map object, that we can plot directly the same way we did above:

[11]:

gc_3fhl_smoothed = gc_3fhl.smooth(kernel="gauss", width=0.2 * u.deg)

[12]:

gc_3fhl_smoothed.plot(stretch="sqrt");

../_images/notebooks_first_steps_23_0.png

The smoothed plot already looks much nicer, but still the image is rather large. As we are mostly interested in the inner part of the image, we will cut out a quadratic region of the size 9 deg x 9 deg around Vela. Therefore we use Map.cutout to make a cutout map:

[13]:

# define center and size of the cutout region
center = SkyCoord(0, 0, unit="deg", frame="galactic")
gc_3fhl_cutout = gc_3fhl_smoothed.cutout(center, 9 * u.deg)
gc_3fhl_cutout.plot(stretch="sqrt");

../_images/notebooks_first_steps_25_0.png

For a more detailed introdcution to ganmmapy.maps, take a look a the intro_maps.ipynb notebook.

Exercises¶

Add a marker and circle at the position of Sag A* (you can find examples in the WCSAxes documentation).

[ ]:

Event lists¶

Almost any high-level gamma-ray data analysis starts with the raw measured counts data, which is stored in event lists. In Gammapy event lists are represented by the gammapy.data.EventList class.

In this section we will learn how to:

Read event lists from FITS files
Access and work with the EventList attributes such as .table and .energy
Filter events lists using convenience methods

Let’s start with the import from the gammapy.data submodule:

[14]:

from gammapy.data import EventList

Very similar to the sky map class an event list can be created, by passing a filename to the .read() method:

[15]:

events_3fhl = EventList.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-events.fits.gz"
)

This time the actual data is stored as an astropy.table.Table object. It can be accessed with .table attribute:

[16]:

events_3fhl.table

[16]:

Table length=32843

ENERGY	RA	DEC	L	B	THETA	PHI	ZENITH_ANGLE	EARTH_AZIMUTH_ANGLE	TIME	EVENT_ID	RUN_ID	RECON_VERSION	CALIB_VERSION [3]	EVENT_CLASS [32]	EVENT_TYPE [32]	CONVERSION_TYPE	LIVETIME	DIFRSP0	DIFRSP1	DIFRSP2	DIFRSP3	DIFRSP4
MeV	deg	deg	deg	deg	deg	deg	deg	deg	s								s
float32	float32	float32	float32	float32	float32	float32	float32	float32	float64	int32	int32	int16	int16	bool	bool	int16	float64	float32	float32	float32	float32	float32
12186.642	260.45935	-33.553337	353.36273	1.7538676	71.977325	125.50694	59.22307	231.79672	239572401.29222104	1823040	239571670	0	0 .. 0	False .. True	False .. True	0	238.57837238907814	0.0	0.0	0.0	0.0	0.0
25496.598	261.37506	-34.395004	353.09607	0.6520652	42.49406	278.49347	41.092773	227.89838	239577842.16217342	550833	239577663	0	0 .. 0	False .. True	False .. False	1	176.16850754618645	0.0	0.0	0.0	0.0	0.0
15621.498	259.56973	-33.409416	353.05673	2.4450684	64.32412	234.22194	66.526794	232.75734	239578244.7997108	1353175	239577663	0	0 .. 0	False .. True	False .. False	1	9.392075657844543	0.0	0.0	0.0	0.0	0.0
12816.32	273.95883	-25.340391	6.45856	-4.0548873	43.292503	142.87392	13.232716	108.02273	239605914.66160735	9636241	239601276	0	0 .. 0	False .. True	False .. False	1	4.034786552190781	0.0	0.0	0.0	0.0	0.0
18988.387	260.8568	-36.355804	351.23734	-0.101912394	26.916113	290.39337	23.8726	212.91147	239611913.14460415	11233188	239606871	0	0 .. 0	False .. True	False .. True	0	131.60132896900177	0.0	0.0	0.0	0.0	0.0
11610.23	266.15518	-26.224436	2.1986027	1.6034819	35.77363	274.53387	23.537594	232.64166	239623554.55414733	14156811	239618329	0	0 .. 0	False .. True	False .. False	1	74.98110938072205	0.0	0.0	0.0	0.0	0.0
13960.802	271.44742	-29.615316	1.6267247	-4.1431155	25.917883	238.0368	15.037035	123.32094	239634549.1748726	14140569	239629788	0	0 .. 0	False .. True	False .. False	1	106.37336817383766	0.0	0.0	0.0	0.0	0.0
10477.372	266.3981	-28.96814	359.97003	-0.011748177	39.091587	275.5457	33.02354	229.59308	239635161.87982982	15688393	239629788	0	0 .. 0	False .. True	False .. True	0	214.62817406654358	0.0	0.0	0.0	0.0	0.0
13030.88	271.70428	-20.632627	9.59348	0.026241468	52.622505	161.3205	39.350845	91.9986	239639873.2076075	1736482	239639436	0	0 .. 0	False .. True	False .. True	0	94.68753063678741	0.0	0.0	0.0	0.0	0.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
387834.72	270.3779	-21.56711	8.171749	0.64531475	56.755512	221.84715	24.358454	86.67913	460185260.7970139	7595626	460181260	0	0 .. 0	False .. True	False .. True	0	34.214694023132324	0.0	0.0	0.0	0.0	0.0
20559.74	268.5538	-26.345692	3.200638	-0.30328986	49.523575	233.67285	12.370642	250.35716	460185803.2027966	8785590	460181260	0	0 .. 0	False .. True	False .. True	0	103.17629969120026	0.0	0.0	0.0	0.0	0.0
27209.146	266.59344	-30.52607	358.72775	-0.96718174	62.1856	140.27434	32.686306	109.74662	460190778.6372646	7274453	460186976	0	0 .. 0	False .. True	False .. True	0	43.22334712743759	0.0	0.0	0.0	0.0	0.0
13911.061	269.30997	-27.239439	2.7684028	-1.3365301	65.15399	224.52101	53.017742	242.62904	460197889.2652691	12049887	460192198	0	0 .. 0	False .. True	False .. True	0	95.46356403827667	0.0	0.0	0.0	0.0	0.0
13226.425	265.16287	-27.344238	0.7796942	1.7680178	59.38332	126.7019	32.12299	246.97205	460203215.10816145	11716888	460198235	0	0 .. 0	False .. True	False .. False	1	3.733097553253174	0.0	0.0	0.0	0.0	0.0
17445.463	266.63342	-28.807201	0.21464892	-0.1039705	55.48627	135.59155	14.227151	106.7812	460225372.9842249	1679082	460224933	0	0 .. 0	False .. True	False .. False	1	80.52235281467438	0.0	0.0	0.0	0.0	0.0
13133.864	270.42474	-22.651058	7.251185	0.071358204	48.704975	134.73102	2.489122	294.48605	460225688.52486295	2879335	460224933	0	0 .. 0	False .. True	False .. False	1	117.88173341751099	0.0	0.0	0.0	0.0	0.0
32095.705	266.0002	-29.77206	359.1034	-0.13615231	45.013103	236.72498	6.92107	212.86594	460231367.1387127	1113706	460231084	0	0 .. 0	False .. True	False .. True	0	108.92976492643356	0.0	0.0	0.0	0.0	0.0
18465.783	266.39728	-29.105953	359.85202	-0.08294058	55.97552	135.87787	18.909636	112.137924	459939497.057684	7689831	459935572	0	0 .. 0	False .. True	False .. True	0	70.72638684511185	0.0	0.0	0.0	0.0	0.0
14457.25	262.72217	-34.388405	353.7184	-0.26906812	45.683174	237.74162	25.728264	240.87035	459945845.4798405	10049680	459941302	0	0 .. 0	False .. True	False .. True	0	147.4274787902832	0.0	0.0	0.0	0.0	0.0

You can do len over event_3fhl.table to find the total number of events.

[17]:

print("Total number of events: {}".format(len(events_3fhl.table)))

Total number of events: 32843

And we can access any other attribute of the Table object as well:

[18]:

events_3fhl.table.colnames

[18]:

['ENERGY',
 'RA',
 'DEC',
 'L',
 'B',
 'THETA',
 'PHI',
 'ZENITH_ANGLE',
 'EARTH_AZIMUTH_ANGLE',
 'TIME',
 'EVENT_ID',
 'RUN_ID',
 'RECON_VERSION',
 'CALIB_VERSION',
 'EVENT_CLASS',
 'EVENT_TYPE',
 'CONVERSION_TYPE',
 'LIVETIME',
 'DIFRSP0',
 'DIFRSP1',
 'DIFRSP2',
 'DIFRSP3',
 'DIFRSP4']

For convenience we can access the most important event parameters as properties on the EventList objects. The attributes will return corresponding Astropy objects to represent the data, such as astropy.units.Quantity, astropy.coordinates.SkyCoord or astropy.time.Time objects:

[19]:

events_3fhl.energy.to("GeV")

[19]:

$$[12.186643,~25.496599,~15.621499,~\dots,~32.095707,~18.465784,~14.457251] \; \mathrm{GeV}$$

[20]:

events_3fhl.galactic
# events_3fhl.radec

[20]:

<SkyCoord (Galactic): (l, b) in deg
    [(353.36228879,  1.75408483), (353.09562941,  0.6522806 ),
     (353.05628243,  2.44528685), ..., (359.10295505, -0.1359316 ),
     (359.85157506, -0.08269984), (353.71795506, -0.26883694)]>

[21]:

events_3fhl.time

[21]:

<Time object: scale='tt' format='mjd' value=[54682.82946153 54682.89243456 54682.89709472 ... 57236.75267735
 57233.37455141 57233.44802852]>

In addition EventList provides convenience methods to filter the event lists. One possible use case is to find the highest energy event within a radius of 0.5 deg around the vela position:

[22]:

# select all events within a radius of 0.5 deg around center
from gammapy.utils.regions import SphericalCircleSkyRegion

region = SphericalCircleSkyRegion(center, radius=0.5 * u.deg)
events_gc_3fhl = events_3fhl.select_region(region)

# sort events by energy
events_gc_3fhl.table.sort("ENERGY")

# and show highest energy photon
events_gc_3fhl.energy[-1].to("GeV")

[22]:

$$1917.8592 \; \mathrm{GeV}$$

Exercises¶

Make a counts energy spectrum for the galactic center region, within a radius of 10 deg.

[ ]:

Source catalogs¶

Gammapy provides a convenient interface to access and work with catalog based data.

In this section we will learn how to:

Load builtins catalogs from gammapy.catalog
Sort and index the underlying Astropy tables
Access data from individual sources

Let’s start with importing the 3FHL catalog object from the gammapy.catalog submodule:

[23]:

from gammapy.catalog import SourceCatalog3FHL

First we initialize the Fermi-LAT 3FHL catalog and directly take a look at the .table attribute:

[24]:

fermi_3fhl = SourceCatalog3FHL()
fermi_3fhl.table

[24]:

Table masked=True length=1556

Source_Name	RAJ2000	DEJ2000	GLON	GLAT	Conf_95_SemiMajor	Conf_95_SemiMinor	Conf_95_PosAng	ROI_num	Signif_Avg	Pivot_Energy	Flux_Density	Unc_Flux_Density	Flux	Unc_Flux	Energy_Flux	Unc_Energy_Flux	Signif_Curve	SpectrumType	Spectral_Index	Unc_Spectral_Index	beta	Unc_beta	PowerLaw_Index	Unc_PowerLaw_Index	Flux_Band [5]	Unc_Flux_Band [5,2]	nuFnu [5]	Sqrt_TS_Band [5]	Npred	HEP_Energy	HEP_Prob	Variability_BayesBlocks	Extended_Source_Name	ASSOC_GAM	TEVCAT_FLAG	ASSOC_TEV	CLASS	ASSOC1	ASSOC2	ASSOC_PROB_BAY	ASSOC_PROB_LR	Redshift	NuPeak_obs
	deg	deg	deg	deg	deg	deg	deg			GeV	1 / (cm2 GeV s)	1 / (cm2 GeV s)	1 / (cm2 s)	1 / (cm2 s)	erg / (cm2 s)	erg / (cm2 s)									1 / (cm2 s)	1 / (cm2 s)	erg / (cm2 s)			GeV													Hz
bytes18	float32	float32	float32	float32	float32	float32	float32	int16	float32	float32	float32	float32	float32	float32	float32	float32	float32	bytes11	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	float32	int16	bytes18	bytes18	bytes1	bytes21	bytes7	bytes26	bytes26	float32	float32	float32	float32
3FHL J0001.2-0748	0.3107	-7.8075	89.0094	-67.3118	0.0424	0.0424	nan	64	5.362	23.73	5.3174e-13	2.0975e-13	2.9593e-11	1.1704e-11	1.6752e-12	1.0743e-12	1.02	PowerLaw	1.6724	0.8274	0.5916	0.7129	2.2226	0.4808	1.1127661e-11 .. 1.1422301e-22	-6.0763976e-12 .. 6.529277e-12	3.533989e-13 .. 1.1789072e-22	3.1458344 .. 0.0	7.63	86.975	0.9964	1		3FGL J0001.2-0748	N		bll	PMN J0001-0746		0.9974	0.9721	nan	306196370000000.0
3FHL J0001.9-4155	0.4849	-41.9303	334.1216	-72.0697	0.1018	0.1018	nan	429	5.638	28.42	5.4253e-13	1.6839e-13	4.3230e-11	1.3428e-11	3.4900e-12	1.8276e-12	0.45	PowerLaw	1.7819	0.4941	0.1187	0.2798	1.9418	0.3100	2.1003905e-11 .. 1.9287885e-18	-8.032091e-12 .. 5.8594097e-12	6.7452245e-13 .. 2.078675e-18	4.899907 .. 0.0	12.51	266.625	0.9622	1		3FGL J0002.2-4152	N		bcu	1RXS J000135.5-415519		0.9960	0.0000	nan	6309576500000000.0
3FHL J0002.1-6728	0.5283	-67.4825	310.0868	-48.9549	0.0357	0.0357	nan	386	8.470	20.82	1.2062e-12	3.2106e-13	5.0093e-11	1.3349e-11	2.3058e-12	9.5580e-13	1.53	PowerLaw	1.8109	0.6260	0.7933	0.5956	2.4285	0.3710	2.4550664e-11 .. 1.9009976e-21	-8.634195e-12 .. 4.8021903e-12	7.7340695e-13 .. 1.9026535e-21	5.900217 .. 0.0	17.11	52.152	0.9988	1		3FGL J0002.0-6722	N		bcu	SUMSS J000215-672653		0.0000	0.9395	nan	4466832000000000.0
3FHL J0003.3-5248	0.8300	-52.8150	318.9245	-62.7936	0.0425	0.0425	nan	145	7.229	23.66	7.5065e-13	2.3102e-13	4.1560e-11	1.2839e-11	2.2874e-12	1.1145e-12	1.70	PowerLaw	1.6010	0.5644	0.9972	0.1721	2.2481	0.3732	2.0886386e-11 .. 7.5867555e-23	-8.143967e-12 .. 5.31299e-12	6.6265456e-13 .. 7.800202e-23	5.298393 .. 0.0	13.02	67.310	0.9636	1		3FGL J0003.2-5246	N		bcu	RBS 0006		0.9996	0.9716	nan	7.079464e+16
3FHL J0007.0+7303	1.7647	73.0560	119.6625	10.4666	0.0101	0.0101	nan	277	75.265	12.80	1.7436e-10	7.5950e-12	1.5308e-09	6.1341e-11	3.6785e-11	1.5973e-12	3.24	LogParabola	3.1751	0.2103	0.9021	0.2659	3.8315	0.1141	1.3514667e-09 .. 3.839895e-18	-5.7581186e-11 .. 4.060418e-12	4.109739e-11 .. 2.9231144e-18	71.33829 .. 0.0	654.15	60.292	0.9972	1		3FGL J0007.0+7302	E	CTA 1	PSR	LAT PSR J0007+7303		1.0000	0.0000	nan	nan
3FHL J0007.9+4711	1.9931	47.1920	115.3093	-15.0354	0.0196	0.0196	nan	302	17.774	17.19	5.9778e-12	8.7683e-13	1.5131e-10	2.2181e-11	5.1444e-12	1.0540e-12	0.56	PowerLaw	2.6783	0.4196	0.1696	0.3282	2.8588	0.2685	1.0582407e-10 .. 1.9819723e-16	-1.7538379e-11 .. 4.823511e-12	3.278615e-12 .. 1.8668298e-16	15.209969 .. 0.0	50.95	68.152	0.9759	1		3FGL J0008.0+4713	N		bll	MG4 J000800+4712		1.0000	0.9873	0.2800	2511884200000000.0
3FHL J0008.4-2339	2.1243	-23.6514	50.2908	-79.7021	0.0366	0.0366	nan	517	9.679	16.96	3.0610e-12	7.3475e-13	7.4602e-11	1.7896e-11	2.4733e-12	8.1716e-13	0.34	PowerLaw	2.7388	0.7145	0.1737	0.5618	2.9070	0.4520	5.804992e-11 .. 1.1117311e-20	-1.4419374e-11 .. 6.10661e-12	1.7951775e-12 .. 1.0403958e-20	9.133706 .. 0.0	19.83	71.122	0.9968	1		3FGL J0008.6-2340	N		bll	RBS 0016		0.9996	0.9673	0.1470	524807800000000.0
3FHL J0009.1+0628	2.2874	6.4814	104.4637	-54.8669	0.0385	0.0385	nan	402	6.282	18.92	1.2691e-12	4.3696e-13	4.1597e-11	1.4317e-11	1.6903e-12	8.9372e-13	0.10	PowerLaw	2.5529	0.8363	0.0122	0.4477	2.5800	0.5391	2.4161059e-11 .. 6.6482124e-19	-9.546595e-12 .. 6.287476e-12	7.566492e-13 .. 6.5095056e-19	4.678369 .. 0.0	10.95	12.256	0.9721	1		3FGL J0009.1+0630	N		bll	CRATES J000903.95+062821.5		0.9993	0.9878	nan	663742400000000.0
3FHL J0009.4+5030	2.3504	50.5049	116.1257	-11.8105	0.0176	0.0176	nan	302	22.402	17.04	9.8252e-12	1.3192e-12	2.2191e-10	2.6212e-11	8.7336e-12	1.2488e-12	3.15	LogParabola	1.4305	0.3505	0.7965	0.3072	2.3610	0.1611	1.16274e-10 .. 9.252794e-17	-1.8225135e-11 .. 4.417993e-12	3.8564165e-12 .. 7.0436765e-17	15.780677 .. 0.0	78.50	72.762	0.9950	2		3FGL J0009.3+5030	C		bll	NVSS J000922+503028		1.0000	0.9698	nan	1412536400000000.0
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
3FHL J2347.9-1630	356.9978	-16.5106	65.5355	-71.8766	0.0288	0.0288	nan	450	9.297	16.28	3.1279e-12	8.0896e-13	6.7585e-11	1.7478e-11	2.0267e-12	6.5608e-13	0.07	PowerLaw	3.1259	0.7781	0.0104	0.5756	3.1324	0.5259	5.2519888e-11 .. 1.0747592e-20	-1.4230782e-11 .. 6.211813e-12	1.6103665e-12 .. 9.768252e-21	8.333468 .. 0.0	17.55	50.215	0.9869	3		3FGL J2348.0-1630	N		fsrq	PKS 2345-16		0.9994	0.9999	0.5760	9332549000000.0
3FHL J2350.5-3006	357.6354	-30.1070	16.7759	-76.3194	0.0491	0.0491	nan	70	6.497	21.20	1.0879e-12	3.2997e-13	4.7039e-11	1.4274e-11	2.2909e-12	1.1390e-12	0.63	PowerLaw	2.1012	0.6173	0.2880	0.4870	2.3678	0.4234	2.1939225e-11 .. 4.5892933e-16	-8.926376e-12 .. 6.097474e-12	6.927891e-13 .. 4.63469e-16	4.0536985 .. 0.0	12.84	49.286	0.9644	1		3FGL J2350.4-3004	N		bll	NVSS J235034-300603		0.9998	0.9218	0.2237	3981075200000000.0
3FHL J2351.5-7559	357.8926	-75.9890	307.6546	-40.5855	0.0650	0.0650	nan	55	6.067	26.82	4.9826e-13	1.6350e-13	3.5689e-11	1.1769e-11	2.3897e-12	1.2622e-12	0.61	PowerLaw	1.8474	0.5802	0.2003	0.3661	2.0816	0.3532	2.3730832e-11 .. 6.9375605e-17	-8.570627e-12 .. 4.7928705e-12	7.578736e-13 .. 7.316245e-17	5.2754674 .. 0.0	12.41	134.721	0.9892	1		3FGL J2351.9-7601	N		bll	SUMSS J235115-760012		0.0000	0.9625	nan	nan
3FHL J2352.1+1753	358.0415	17.8865	103.5764	-42.7466	0.0838	0.0838	nan	185	4.117	16.97	9.9227e-13	4.3475e-13	2.4254e-11	1.0640e-11	7.6327e-13	4.2356e-13	0.02	PowerLaw	3.0175	1.2164	0.0100	0.8524	3.0166	0.8270	1.5997077e-11 .. 2.9107688e-20	-7.581037e-12 .. 5.821708e-12	4.926488e-13 .. 2.6849966e-20	3.5496242 .. 0.0	6.73	43.107	0.9668	1		3FGL J2352.0+1752	N		bll	CLASS J2352+1749		0.9926	0.0000	nan	1737799900000000.0
3FHL J2356.2+4035	359.0746	40.5985	111.7521	-21.0732	0.0298	0.0298	nan	312	7.625	29.01	5.2427e-13	1.5104e-13	4.3400e-11	1.2511e-11	3.6677e-12	1.8547e-12	0.35	PowerLaw	2.0233	0.4242	-0.0706	0.1926	1.9095	0.2975	2.5777725e-11 .. 3.110794e-16	-8.514681e-12 .. 5.4134618e-12	8.2889175e-13 .. 3.3694582e-16	6.2127647 .. 0.0	13.81	417.861	0.9119	1		3FGL J2356.0+4037	N		bll	NVSS J235612+403648		0.9998	0.9199	0.1310	6309576500000000.0
3FHL J2357.4-1717	359.3690	-17.2996	68.4009	-74.1285	0.0327	0.0327	nan	450	6.961	29.52	5.4394e-13	1.7370e-13	4.6654e-11	1.4945e-11	3.7598e-12	1.9583e-12	1.11	PowerLaw	1.5762	0.5187	0.3513	0.3771	1.9430	0.3116	1.9003682e-11 .. 2.714288e-20	-8.131149e-12 .. 6.4742196e-12	6.1025685e-13 .. 2.92465e-20	4.552822 .. 0.0	12.30	146.757	0.9838	1		3FGL J2357.4-1716	N		bll	RBS 2066		0.9999	0.9631	nan	8.912525e+16
3FHL J2358.4-1808	359.6205	-18.1408	66.5520	-74.8501	0.0511	0.0511	nan	450	6.493	18.23	1.6335e-12	4.9686e-13	4.8680e-11	1.4811e-11	1.7825e-12	7.6480e-13	1.83	PowerLaw	2.0532	0.6673	0.9999	0.0134	2.7312	0.5024	2.6735683e-11 .. 6.0349635e-21	-9.960717e-12 .. 6.2551535e-12	8.323882e-13 .. 5.7844478e-21	4.3616037 .. 0.0	12.74	28.304	0.9845	1		3FGL J2358.6-1809	N					0.0000	0.0000	nan	nan
3FHL J2358.5+3829	359.6266	38.4963	111.6905	-23.2173	0.0584	0.0584	nan	312	5.797	18.24	1.4104e-12	4.4534e-13	4.2106e-11	1.3321e-11	1.7404e-12	9.8271e-13	0.44	PowerLaw	2.7466	0.6917	-0.1329	0.3013	2.5576	0.5781	2.824428e-11 .. 9.750687e-17	-9.458818e-12 .. 5.2791343e-12	8.852925e-13 .. 9.5778846e-17	5.7128677 .. 0.0	13.13	57.301	0.9782	1		3FGL J2358.5+3827	N		bcu	B3 2355+382		0.0000	0.9254	nan	nan
3FHL J2359.1-3038	359.7760	-30.6397	12.7909	-78.0268	0.0231	0.0231	nan	70	11.551	21.21	1.8903e-12	4.1965e-13	8.1774e-11	1.8149e-11	4.2849e-12	1.6806e-12	0.08	PowerLaw	2.2865	0.4632	0.0101	0.2434	2.2944	0.3092	5.5015617e-11 .. 6.037456e-17	-1.3604539e-11 .. 8.488618e-12	1.7422797e-12 .. 6.164239e-17	9.39347 .. 0.0	22.41	111.366	0.9607	1		3FGL J2359.3-3038	P	H 2356-309	bll	H 2356-309		0.9999	0.9975	0.1650	2.818388e+17
3FHL J2359.3-2049	359.8293	-20.8256	58.0522	-76.5411	0.0722	0.0722	nan	580	4.638	19.02	9.1911e-13	3.6043e-13	3.0559e-11	1.1979e-11	1.2593e-12	7.4704e-13	0.32	PowerLaw	2.3402	0.9445	0.1851	0.6600	2.5615	0.5838	2.3253791e-11 .. 8.3778735e-21	-8.939083e-12 .. 6.2386546e-12	7.2875863e-13 .. 8.224765e-21	4.8207045 .. 0.0	8.06	64.177	0.9859	1		3FGL J2359.5-2052	N		bll	TXS 2356-210		0.9894	0.9906	0.0960	4073799600000000.0

This looks very familiar again. The data is just stored as an astropy.table.Table object. We have all the methods and attributes of the Table object available. E.g. we can sort the underlying table by Signif_Avg to find the top 5 most significant sources:

[25]:

# sort table by significance
fermi_3fhl.table.sort("Signif_Avg")

# invert the order to find the highest values and take the top 5
top_five_TS_3fhl = fermi_3fhl.table[::-1][:5]

# print the top five significant sources with association and source class
top_five_TS_3fhl[["Source_Name", "ASSOC1", "ASSOC2", "CLASS", "Signif_Avg"]]

[25]:

Table masked=True length=5

Source_Name	ASSOC1	ASSOC2	CLASS	Signif_Avg
bytes18	bytes26	bytes26	bytes7	float32
3FHL J0534.5+2201	Crab Nebula		PWN	168.641
3FHL J1104.4+3812	Mkn 421		BLL	144.406
3FHL J0835.3-4510	PSR J0835-4510	Vela X field	PSR	138.801
3FHL J0633.9+1746	PSR J0633+1746		PSR	99.734
3FHL J1555.7+1111	PG 1553+113		BLL	94.411

If you are interested in the data of an individual source you can access the information from catalog using the name of the source or any alias source name that is defined in the catalog:

[26]:

mkn_421_3fhl = fermi_3fhl["3FHL J1104.4+3812"]

# or use any alias source name that is defined in the catalog
mkn_421_3fhl = fermi_3fhl["Mkn 421"]
print(mkn_421_3fhl.data["Signif_Avg"])

144.40611

Exercises¶

Try to load the Fermi-LAT 2FHL catalog and check the total number of sources it contains.
Select all the sources from the 2FHL catalog which are contained in the Galactic Center region. The methods `WcsGeom.contains() <../api/gammapy.maps.WcsGeom.rst#gammapy.maps.WcsGeom.contains>`__ and `SourceCatalog.positions <../api/gammapy.catalog.SourceCatalog.rst#gammapy.catalog.SourceCatalog.positions>`__ might be helpful for this. Add markers for all these sources and try to add labels with the source names. The function ax.text() might be also helpful.
Try to find the source class of the object at position ra=68.6803, dec=9.3331

[ ]:

Spectral models and flux points¶

In the previous section we learned how access basic data from individual sources in the catalog. Now we will go one step further and explore the full spectral information of sources. We will learn how to:

Plot spectral models
Compute integral and energy fluxes
Read and plot flux points

As a first example we will start with the Crab Nebula:

[27]:

crab_3fhl = fermi_3fhl["Crab Nebula"]
print(crab_3fhl.spectral_model)

PowerLaw

Parameters:

           name     value     error        unit      min max frozen
        --------- --------- --------- -------------- --- --- ------
            index 2.220e+00 2.498e-02                nan nan  False
        amplitude 1.713e-10 3.389e-12 cm-2 GeV-1 s-1 nan nan  False
        reference 2.273e+01 0.000e+00            GeV nan nan   True

Covariance:

           name     index   amplitude reference
        --------- --------- --------- ---------
            index 6.241e-04 0.000e+00 0.000e+00
        amplitude 0.000e+00 1.148e-23 0.000e+00
        reference 0.000e+00 0.000e+00 0.000e+00

The crab_3fhl.spectral_model is an instance of the gammapy.spectrum.models.PowerLaw2 model, with the parameter values and errors taken from the 3FHL catalog.

Let’s plot the spectral model in the energy range between 10 GeV and 2000 GeV:

[28]:

ax_crab_3fhl = crab_3fhl.spectral_model.plot(
    energy_range=[10, 2000] * u.GeV, energy_power=0
)

../_images/notebooks_first_steps_59_0.png

We assign the return axes object to variable called ax_crab_3fhl, because we will re-use it later to plot the flux points on top.

To compute the differential flux at 100 GeV we can simply call the model like normal Python function and convert to the desired units:

[29]:

crab_3fhl.spectral_model(100 * u.GeV).to("cm-2 s-1 GeV-1")

[29]:

$$6.3848913 \times 10^{-12} \; \mathrm{\frac{1}{GeV\,s\,cm^{2}}}$$

Next we can compute the integral flux of the Crab between 10 GeV and 2000 GeV:

[30]:

crab_3fhl.spectral_model.integral(emin=10 * u.GeV, emax=2000 * u.GeV).to(
    "cm-2 s-1"
)

[30]:

$$8.6745734 \times 10^{-9} \; \mathrm{\frac{1}{s\,cm^{2}}}$$

We can easily convince ourself, that it corresponds to the value given in the Fermi-LAT 3FHL catalog:

[31]:

crab_3fhl.data["Flux"]

[31]:

$$8.6589091 \times 10^{-9} \; \mathrm{\frac{1}{s\,cm^{2}}}$$

In addition we can compute the energy flux between 10 GeV and 2000 GeV:

[32]:

crab_3fhl.spectral_model.energy_flux(emin=10 * u.GeV, emax=2000 * u.GeV).to(
    "erg cm-2 s-1"
)

[32]:

$$5.3114891 \times 10^{-10} \; \mathrm{\frac{erg}{s\,cm^{2}}}$$

Next we will access the flux points data of the Crab:

[33]:

print(crab_3fhl.flux_points)

FluxPoints(sed_type="flux", n_points=5)

If you want to learn more about the different flux point formats you can read the specification here.

No we can check again the underlying astropy data structure by accessing the .table attribute:

[34]:

crab_3fhl.flux_points.table

[34]:

Table length=5

e_min	e_max	flux	flux_errn	flux_errp	e2dnde	e2dnde_errn	e2dnde_errp	is_ul	flux_ul	e2dnde_ul	sqrt_ts
GeV	GeV	1 / (cm2 s)	1 / (cm2 s)	1 / (cm2 s)	erg / (cm2 s)	erg / (cm2 s)	erg / (cm2 s)		1 / (cm2 s)	erg / (cm2 s)
float64	float64	float32	float32	float32	float32	float32	float32	bool	float64	float64	float32
10.0	20.0	5.1698894e-09	1.3338798e-10	1.3338798e-10	1.6420408e-10	4.236619e-12	4.236619e-12	False	nan	nan	125.15745
20.0	50.0	2.2450237e-09	8.671894e-11	8.671894e-11	1.1808442e-10	4.561268e-12	4.561268e-12	False	nan	nan	88.71535
50.0	150.0	9.2431746e-10	5.497474e-11	5.497474e-11	1.08686914e-10	6.464267e-12	6.464267e-12	False	nan	nan	59.087498
150.0	500.0	2.7589558e-10	2.9164688e-11	3.1360175e-11	9.2301645e-11	9.757128e-12	1.0491635e-11	False	nan	nan	33.076164
500.0	2000.0	6.6840214e-11	1.4629655e-11	1.6916293e-11	6.9011685e-11	1.5104936e-11	1.746586e-11	False	nan	nan	15.573053

Finally let’s combine spectral model and flux points in a single plot and scale with energy_power=2 to obtain the spectral energy distribution:

[35]:

ax = crab_3fhl.spectral_model.plot(
    energy_range=[10, 2000] * u.GeV, energy_power=2
)
crab_3fhl.flux_points.to_sed_type("dnde").plot(ax=ax, energy_power=2);

../_images/notebooks_first_steps_73_0.png

Exercises¶

Plot the spectral model and flux points for PKS 2155-304 for the 3FGL and 2FHL catalogs. Try to plot the error of the model (aka “Butterfly”) as well. Note this requires the uncertainties package to be installed on your machine.

[ ]:

What next?¶

This was a quick introduction to some of the high-level classes in Astropy and Gammapy.

To learn more about those classes, go to the API docs (links are in the introduction at the top).
To learn more about other parts of Gammapy (e.g. Fermi-LAT and TeV data analysis), check out the other tutorial notebooks.
To see what’s available in Gammapy, browse the Gammapy docs or use the full-text search.
If you have any questions, ask on the mailing list.