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

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 (2FHL) 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:

If you’re not yet familiar with the listed Astropy classes, maybe check out the Astropy introduction for Gammapy users first.

Setup

Important: to run this tutorial the environment variable GAMMAPY_EXTRA must be defined and point to the directory, where the gammapy-extra is respository located on your machine. To check whether your setup is correct you can execute the following cell:

In [1]:

import os

path = os.path.expandvars('$GAMMAPY_EXTRA/datasets')

if not os.path.exists(path):
    raise Exception('gammapy-extra 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.

In [2]:

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

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

In [3]:

%matplotlib inline
import matplotlib.pyplot as plt

In [4]:

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

Sky images

The central data structure to work with image based data in Gammapy is the SkyImage class. It combines the raw data with world coordinate (WCS) information, FITS I/O functionality and convenience methods, that allow easy handling, processing and plotting of image based data.

In this section we will learn how to:

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

The SkyImage class is part of the gammapy.image submodule. So we will start by importing it from there:

In [5]:

from gammapy.image import SkyImage

As a first example, we will read a FITS file from a prepared Fermi-LAT 2FHL dataset:

In [6]:

vela_2fhl = SkyImage.read("$GAMMAPY_EXTRA/datasets/fermi_2fhl/fermi_2fhl_vela.fits.gz", hdu='COUNTS')

As the FITS file fermi_2fhl_vela.fits.gz contains mutiple image extensions and a SkyImage can only represent a single image, we explicitely specify to read the extension called ‘Counts’. Let’s print the image to get some basic information about it:

In [7]:

print(vela_2fhl)

Name: Counts
Data shape: (180, 320)
Data type: >f8
Data unit: ct
Data mean: 2.720e-02
WCS type: ['GLON-CAR', 'GLAT-CAR']

The shape of the image is 320x180 pixel, the data unit is counts (‘ct’) and it is defined using a cartesian projection in galactic coordinates.

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

In [8]:

vela_2fhl.data

Out[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.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.],
       [ 0.,  0.,  0., ...,  0.,  0.,  0.]])

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

In [9]:

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

Total number of counts in the image: 1567

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

In [10]:

vela_2fhl.show()
../_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 SkyImage.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 SkyImage object, that we can plot directly the same way we did above:

In [11]:

# smooth counts image with gaussian kernel of 0.5 deg
vela_2fhl_smoothed = vela_2fhl.smooth(kernel='gauss', radius=0.5 * u.deg)
vela_2fhl_smoothed.show()
../_images/notebooks_first_steps_22_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 SkyImage.cutout(), which wraps Astropy’s Cutout2D. Again it returns a SkyImage object:

In [12]:

# define center and size of the cutout region
center = SkyCoord(265.0, -2.0, unit='deg', frame='galactic')
size = 9.0 * u.deg

vela_2fhl_cutout = vela_2fhl_smoothed.cutout(center, size)
vela_2fhl_cutout.show()
../_images/notebooks_first_steps_24_0.png

To make this exercise a bit more scientifically useful, we will load a second image containing WMAP data from the same region:

In [13]:

vela_wmap = SkyImage.read("$GAMMAPY_EXTRA/datasets/images/Vela_region_WMAP_K.fits")

# define a norm to stretch the data, so it is better visible
norm = simple_norm(vela_wmap.data, stretch='sqrt', max_percent=99.9)
vela_wmap.show(cmap='viridis', norm=norm)
../_images/notebooks_first_steps_26_0.png

In order to make the structures in the data better visible we used the simple_norm() method, which allows to stretch the data for plotting. This is very similar to the methods that e.g. ds9 provides. In addition we used a different colormap called ‘viridis’. An overview of different colomaps can be found here.

Now let’s check the details of the WMAP image:

In [14]:

print(vela_wmap)

Name: None
Data shape: (300, 300)
Data type: >f4
Data unit: mK
Data mean: 1.827e+00
WCS type: ['RA---TAN', 'DEC--TAN']

As you can see it is defined using a tangential projection and ICRS coordinates, which is different from the projection used for the vela_2fhl image. To compare both images we have to reproject one image to the WCS of the other. This can be achieved with the SkyImage.reproject() method:

In [15]:

# reproject WMAP image
vela_wmap_reprojected = vela_wmap.reproject(vela_2fhl)

# cutout part we're interested in
vela_wmap_reprojected_cutout = vela_wmap_reprojected.cutout(center, size)
vela_wmap_reprojected_cutout.show(cmap='viridis', norm=norm)
../_images/notebooks_first_steps_30_0.png

Finally we will combine both images in single plot, by plotting WMAP contours on top of the smoothed Fermi-LAT 2FHL image:

In [16]:

fig, ax, _ = vela_2fhl_cutout.plot()
ax.contour(vela_wmap_reprojected_cutout.data, cmap='Blues')

Out[16]:

<matplotlib.contour.QuadContourSet at 0x11639b278>
../_images/notebooks_first_steps_32_1.png

Exercises

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:

In [17]:

from gammapy.data import EventList

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

In [18]:

events_2fhl = EventList.read('$GAMMAPY_EXTRA/datasets/fermi_2fhl/2fhl_events.fits.gz')

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

In [19]:

events_2fhl.table

Out[19]:
<Table length=60978>
ENERGYRADECLBTHETAPHIZENITH_ANGLEEARTH_AZIMUTH_ANGLETIMEEVENT_IDRUN_IDRECON_VERSIONCALIB_VERSION [3]EVENT_CLASS [32]EVENT_TYPE [32]CONVERSION_TYPELIVETIMEDIFRSP0DIFRSP1DIFRSP2DIFRSP3DIFRSP4
MeVdegdegdegdegdegdegdegdegss
float32float32float32float32float32float32float32float32float32float64int32int32int16int16boolboolint16float64float32float32float32float32float32
145927.0291.66242.234174.543711.867838.045583.535855.6387314.034239561457.866485143723955956500 .. 0False .. TrueFalse .. True0275.6980889742.69657e-113238.250.00.00.0
221273.046.9858-40.6389247.489-58.873934.1051224.20968.2524198.319239562739.085752143223955956500 .. 0False .. TrueFalse .. True064.79749315981.53573e-122774.20.00.00.0
57709.2121.84149.2288169.86832.301771.563634.292536.717325.5439239563180.302869069323955956500 .. 0False .. TrueFalse .. True030.572186478.11096e-12253.2210.00.00.0
221224.083.5626-4.21744207.783-19.077120.508992.160532.3033239.141239563382.213920842423955956500 .. 0False .. TrueFalse .. True027.41250959041.66075e-112980.120.00.00.0
698997.0320.895-1.3278951.2218-33.971835.3621158.74112.086772.2029239566572.951248048323956564500 .. 0False .. TrueFalse .. True0106.4754811232.26543e-132706.140.00.00.0
119159.0318.81112.302862.6361-24.41626.5896213.89417.715623.9409239572348.06172527623957167000 .. 0False .. TrueFalse .. True0185.3464272922.75677e-113439.290.00.00.0
56175.6279.25147.883576.691522.073929.103461.004862.1731321.104239572763.4317601723957273600 .. 0False .. TrueFalse .. True024.45073395972.08265e-104071.250.00.00.0
1.41812e+06100.311-47.4481256.468-21.264161.2256294.1890.4753144.149239573788.813178156923957273600 .. 0False .. TrueFalse .. True068.2716146413.20022e-141086.230.00.00.0
62164.9331.492-41.2264359.42-53.404928.1408229.92752.0142189.054239578601.168200070023957766300 .. 0False .. TrueFalse .. True090.3322626651.21479e-103566.070.00.00.0
.....................................................................
51296.879.031178.9327133.86822.239951.1726267.81499.7602357.79444418973.822128485344441859000 .. 0False .. TrueFalse .. True0181.68898352.14127e-102586.880.00.00.0
60315.7243.681-50.5669332.430.2826824.8501255.68374.8804195.14444421777.761658294944441859000 .. 0False .. TrueFalse .. False157.33589029312.35854e-083541.040.00.00.0
90000.1282.698-33.19872.66219-14.33398.49144343.55947.3965191.033444422182.738750746644441859000 .. 0False .. TrueFalse .. True0187.6024691463.62123e-103792.570.00.00.0
61988.7247.98-48.1091336.1460.022033933.6499144.27180.3966221.396444422758.145869914744441859000 .. 0False .. TrueFalse .. True039.20983272792.26421e-083395.070.00.00.0
54282.9159.924-58.3628286.3620.20555424.1196287.13865.0847177.524444425794.794252223744442461900 .. 0False .. TrueFalse .. True0196.2096688754.8258e-093672.780.00.00.0
146728.0244.848-46.5876335.7542.6032645.25395.3364491.871137.032444425911.819270521044442461900 .. 0False .. TrueFalse .. False165.63278573754.4595e-102899.330.00.00.0
135433.083.527827.9053179.529-2.6910617.491152.085853.306313.822444430968.96893918844443059900 .. 0False .. TrueFalse .. False197.64468312261.57989e-103515.60.00.00.0
61592.1231.214-5.45521357.43540.647346.6356141.04762.7584256.631444433607.906544394444443059900 .. 0False .. TrueFalse .. False128.72906577592.82028e-102793.450.00.00.0
80480.8228.244-45.044327.56110.952837.3149193.71483.3882221.57444433702.069561244344443059900 .. 0False .. TrueFalse .. False1122.8918192982.26095e-103176.430.00.00.0
124449.0238.008-51.0371329.4292.300632.4522199.50480.9768214.48444433764.433572336144443059900 .. 0False .. TrueFalse .. True0185.2554872048.54908e-103366.30.00.00.0

Let’s try to find the total number of events contained int the list. This doesn’t work:

In [20]:

print('Total number of events: {}'.format(len(events_2fhl.table)))

Total number of events: 60978

Because Gammapy objects don’t redefine properties, that are accessible from the underlying Astropy of Numpy data object.

So in this case of course we can directly use the .table attribute to find the total number of events:

In [21]:

print('Total number of events: {}'.format(len(events_2fhl.table)))

Total number of events: 60978

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

In [22]:

events_2fhl.table.colnames

Out[22]:

['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:

In [23]:

events_2fhl.energy.to('GeV')

Out[23]:
$[145.92725,~221.27338,~57.709244,~\dots, 61.592117,~80.480789,~124.44917] \; \mathrm{GeV}$
In [24]:

events_2fhl.galactic
#events_2fhl.radec

Out[24]:

<SkyCoord (Galactic): (l, b) in deg
    [(  74.54326949,  11.86794629), ( 247.48874019, -58.87409431),
     ( 169.86765712,  32.30144389), ..., ( 357.43496217,  40.64753602),
     ( 327.56037412,  10.95297324), ( 329.42901812,   2.30075607)]>

In [25]:

events_2fhl.time

Out[25]:

<Time object: scale='tt' format='mjd' value=[ 54682.7028015   54682.71763043  54682.72273711 ...,  57053.90824179
  57053.90933163  57053.91005343]>

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:

In [26]:

# select all events within a radius of 0.5 deg around center
events_vela_2fhl = events_2fhl.select_sky_cone(center=center, radius=0.5 * u.deg)

# sort events by energy
events_vela_2fhl.table.sort('ENERGY')

# and show highest energy photon
events_vela_2fhl.energy[-1].to('GeV')

Out[26]:
$342.76625 \; \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 2FHL catalog object from the gammapy.catalog submodule:

In [27]:

from gammapy.catalog import SourceCatalog2FHL

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

In [28]:

fermi_2fhl = SourceCatalog2FHL('$GAMMAPY_EXTRA/datasets/catalogs/fermi/gll_psch_v08.fit.gz')
fermi_2fhl.table

Out[28]:
<Table length=360>
Source_NameRAJ2000DEJ2000GLONGLATPos_err_68TSSpectral_IndexUnc_Spectral_IndexIntr_Spectral_Index_D11Unc_Intr_Spectral_Index_D11Intr_Spectral_Index_G12Unc_Intr_Spectral_Index_G12Flux50Unc_Flux50Energy_Flux50Unc_Energy_Flux50Flux50_171GeVUnc_Flux50_171GeV [2]Sqrt_TS50_171GeVFlux171_585GeVUnc_Flux171_585GeV [2]Sqrt_TS171_585GeVFlux585_2000GeVUnc_Flux585_2000GeV [2]Sqrt_TS585_2000GeVNpredHEP_EnergyHEP_ProbROIASSOCASSOC_PROB_BAYASSOC_PROB_LRCLASSRedshiftNuPeak_obs3FGL_Name1FHL_NameTeVCat_Name
str18float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32int16str25float32float32str8float32float32str18str18str18
2FHL J0008.1+47092.043747.1642115.339-15.06880.06111428.646.242.753.963.192.164.211.23e-116.71e-121.21e-126.71e-133.36344e-12-1.53668e-12 .. 2.16095e-125.353544.34947e-18nan .. 4.75389e-120.08.44871e-18nan .. 7.29101e-120.04.068.150.991MG4 J000800+47120.997220.834827bll2.12.51188e+153FGL J0008.0+47131FHL J0007.7+4709
2FHL J0009.3+50312.343550.5217116.124-11.79320.045439253.975.081.66nannannannan1.91e-117.82e-122.03e-128.79e-138.36282e-12-2.98962e-12 .. 3.89512e-127.351322.91458e-17nan .. 5.10003e-120.03.50875e-16nan .. 4.87458e-120.06.472.761.01NVSS J000922+5030280.9997240.734808bll0.01.41254e+153FGL J0009.3+50301FHL J0009.2+5032
2FHL J0018.5+29474.635529.7879114.463-32.54240.037093630.892.580.992.411.042.41.041.06e-116.15e-122.05e-121.72e-129.65438e-12-5.55342e-12 .. 5.55342e-125.765791.179e-15nan .. 5.37865e-120.01.60465e-16nan .. 6.12008e-120.03.0127.321.03RBS 00420.9998680.978522bll0.15.91561e+163FGL J0018.4+29471FHL J0018.6+2946
2FHL J0022.0+00065.50010.1059107.172-61.86180.051185229.961.860.570.950.720.880.711.97e-119.56e-126.86e-125.29e-121.61197e-11-7.17255e-12 .. 9.95349e-125.38223.63922e-12-4.38636e-12 .. 4.38636e-121.385398.42356e-16nan .. 7.3424e-120.04.8180.130.8625BZGJ0022+00060.999280.900089bll-g0.3064.31519e+16
2FHL J0033.6-19218.4115-19.357594.28-81.22240.0348384148.313.320.692.560.882.330.925.46e-111.5e-117.62e-122.69e-124.00161e-11-1.01615e-11 .. 1.22378e-1112.17252.13901e-12nan .. 8.25812e-120.9079582.43955e-17nan .. 6.84226e-120.013.8170.010.992KUV 00311-19380.9997590.981424bll0.618.31764e+153FGL J0033.6-19211FHL J0033.6-1921TeV J0033-1921
2FHL J0035.8+59498.962559.8312120.972-2.981220.0316026402.42.230.21nannannannan1.25e-101.9e-113.11e-117.23e-129.42294e-11-1.52566e-11 .. 1.71572e-1117.31242.66991e-11-7.62818e-12 .. 9.45084e-1210.39093.72609e-16nan .. 4.67444e-120.046.5247.620.9641ES 0033+5950.9999580.992868bll0.01.31826e+173FGL J0035.9+59491FHL J0035.9+5950TeV J0035+5950
2FHL J0040.3+404910.094940.8315120.676-21.99180.035514826.762.120.81nannannannan1.05e-116.3e-122.84e-122.67e-127.41577e-12-4.0778e-12 .. 6.45182e-124.954423.00308e-12-2.31775e-12 .. 4.46644e-121.708354.73166e-16nan .. 5.79809e-120.03.2258.770.853B3 0037+4050.9987380.9347bcu I0.01.03FGL J0040.3+40491FHL J0040.3+4049
2FHL J0043.9+342410.975534.4109121.164-28.4350.059241739.54.571.613.461.852.951.911.83e-118.24e-122.03e-129.93e-139.65762e-12-4.31687e-12 .. 4.31687e-126.311528.75861e-17nan .. 5.24335e-120.09.43812e-16nan .. 5.68787e-120.05.4109.970.93GB6 J0043+34260.9985940.838844fsrq0.9666.45655e+133FGL J0043.8+34251FHL J0043.7+3425
2FHL J0045.2+212711.316121.4555121.017-41.39330.0378046110.433.070.64nannannannan4.14e-111.26e-116.29e-122.45e-123.11532e-11-8.93008e-12 .. 1.10103e-1110.24042.99712e-12-3.03054e-12 .. 3.03054e-122.70894.68348e-17nan .. 6.42314e-120.011.4246.750.995GB6 J0045+21270.999610.965127bll0.01.0023e+163FGL J0045.3+21261FHL J0045.2+2126
2FHL J0048.0+544912.002354.8281122.433-8.040330.057205535.151.30.51nannannannan1.21e-116.12e-127.64e-125.72e-126.43137e-12-3.56228e-12 .. 5.62894e-124.64012.94345e-12-3.03542e-12 .. 3.03542e-122.275613.59026e-12-2.51846e-12 .. 4.90533e-122.862034.1605.90.9941RXS J004754.5+544758nan0.906161bcu II0.07.87046e+153FGL J0047.9+5447
.....................................................................................................................
2FHL J2322.5+3436350.62834.613102.876-24.77180.071799828.772.330.82.130.852.110.851.42e-117.46e-123.24e-122.48e-121.37062e-11-5.99053e-12 .. 8.27325e-125.511091.50445e-15nan .. 5.25568e-120.03.34396e-16nan .. 6.35236e-120.04.1166.20.99133TXS 2320+3430.9966390.877901bll-g0.0981.04713e+153FGL J2322.5+34361FHL J2322.5+3436
2FHL J2322.7-4916350.698-49.2706334.608-62.04920.03894936.844.581.69nannannannan1.75e-118.3e-121.94e-129.96e-139.13009e-12-3.73356e-12 .. 5.01738e-126.052236.65757e-17nan .. 5.27137e-120.06.34713e-16nan .. 5.77243e-120.05.1106.330.9138SUMSS J232254-4916290.9994340.937726bll0.04.46683e+153FGL J2322.9-4917
2FHL J2323.4+5848350.87458.808111.744-2.140310.0348384183.622.450.3nannannannan7.81e-111.53e-111.64e-114.68e-125.93689e-11-1.21101e-11 .. 1.38986e-1111.62381.03761e-11-4.37845e-12 .. 6.14398e-125.814352.90309e-12-2.02416e-12 .. 3.98202e-123.5083428.8662.121.0135Cas Anan0.999833snrnannan3FGL J2323.4+58491FHL J2323.3+5849TeV J2323+588
2FHL J2323.8+4209350.97242.1629106.059-17.82240.045601670.614.211.164.141.154.131.153.05e-111.05e-113.55e-121.35e-121.62017e-11-5.24708e-12 .. 6.61912e-128.119382.00018e-12-1.41271e-12 .. 2.77927e-122.568492.18906e-18nan .. 5.54082e-120.09.3230.90.921331ES 2321+4190.999560.961506bll0.0594.84173e+153FGL J2323.9+42111FHL J2323.8+4210
2FHL J2324.7-4041351.195-40.6934350.157-67.58360.035865880.933.180.832.940.872.90.872.97e-111.08e-114.34e-122.03e-122.13391e-11-7.09939e-12 .. 9.10435e-128.853832.71834e-12-2.11714e-12 .. 4.11668e-121.657923.54717e-17nan .. 6.3514e-120.08.0118.481.01381ES 2322-4090.9994670.985516bll0.1735878.31764e+153FGL J2324.7-40401FHL J2324.6-4041
2FHL J2329.2+3754352.30837.9157105.532-22.16270.03853939.43.291.192.991.282.921.281.43e-117.17e-122.02e-121.25e-121.11304e-11-4.6962e-12 .. 6.5288e-126.311494.40415e-16nan .. 5.47737e-120.02.64065e-17nan .. 5.92547e-120.04.3142.130.98133NVSS J232914+3754140.9997920.956363bll0.2648.91251e+153FGL J2329.2+37541FHL J2329.1+3754
2FHL J2340.8+8014355.20380.2447119.83817.79990.034893469.794.21.284.01.363.941.371.95e-117.24e-122.26e-129.72e-131.16398e-11-3.81183e-12 .. 4.82791e-128.387451.25703e-16nan .. 3.71812e-120.01.13621e-15nan .. 4.04723e-120.07.980.850.861251RXS J234051.4+8015130.9998740.979435bll0.2743.42768e+153FGL J2340.7+80161FHL J2341.3+8016
2FHL J2343.5+3438355.89834.6484107.422-26.17170.071505627.434.992.364.882.434.792.441.33e-117.32e-121.41e-128.56e-135.94254e-12-2.73222e-12 .. 3.85247e-125.23475.08627e-17nan .. 8.36317e-120.03.14709e-16nan .. 5.6852e-120.03.972.661.01331RXS J234332.5+3439570.9991480.955429bll0.3661.0292e+153FGL J2343.7+3437
2FHL J2347.1+5142356.7551.7081112.88-9.901590.0335464209.841.850.251.740.261.730.267.48e-111.54e-112.62e-118.56e-124.78884e-11-1.15596e-11 .. 1.339e-1111.2822.42894e-11-7.68104e-12 .. 9.87777e-128.85773.51388e-12-3.55282e-12 .. 3.55282e-122.6657225.2652.310.9911ES 2344+5140.999860.988082bll0.0447.07946e+153FGL J2347.0+51421FHL J2347.0+5142TeV J2346+5142
2FHL J2352.0-7558358.019-75.9718307.623-40.61110.073373328.493.371.24nannannannan1.18e-116.14e-121.63e-121.03e-129.00308e-12-3.9133e-12 .. 5.47729e-125.389353.90782e-16nan .. 4.95279e-120.03.23398e-17nan .. 5.16722e-120.03.9134.721.0137nannanunknannan3FGL J2351.9-7601

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 TS to find the top 5 most significant sources:

In [29]:

# sort table by TS
fermi_2fhl.table.sort('TS')

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

# print the top five significant sources with association and source class
top_five_TS_2fhl[['Source_Name', 'ASSOC', 'CLASS']]

Out[29]:
<Table length=5>
Source_NameASSOCCLASS
str18str25str8
2FHL J1104.4+3812Mkn 421bll
2FHL J0534.5+2201Crabpwn
2FHL J1653.9+3945Mkn 501bll
2FHL J1555.7+1111PG 1553+113bll
2FHL J2158.8-3013PKS 2155-304bll

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:

In [30]:

mkn_421_2fhl = fermi_2fhl['2FHL J1104.4+3812']

# or use any alias source name that is defined in the catalog
mkn_421_2fhl = fermi_2fhl['Mkn 421']
print(mkn_421_2fhl.data['TS'])

5600.98

Exercises

  • Try to load the Fermi-LAT 3FHL catalog and check the total number of sources it contains.
  • Select all the sources from the 2FHL catalog which are contained in the Vela region. Add markers for all these sources and try to add labels with the source names. The methods SkyImage.contains() and ax.text() might be 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:

In [31]:

crab_2fhl = fermi_2fhl['Crab']
print(crab_2fhl.spectral_model)

PowerLaw2
ParameterList
Parameter(name='amplitude', value=1.3099999973675835e-09, unit=Unit("1 / (cm2 s)"), min=0, max=None, frozen=False)
Parameter(name='index', value=2.130000114440918, unit=Unit(dimensionless), min=0, max=None, frozen=False)
Parameter(name='emin', value=0.05, unit=Unit("TeV"), min=None, max=None, frozen=False)
Parameter(name='emax', value=2.0, unit=Unit("TeV"), min=None, max=None, frozen=False)

Covariance: [[  4.66488967e-21   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   4.90000004e-03   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]]

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

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

In [32]:

ax_crab_2fhl = crab_2fhl.spectral_model.plot(
    energy_range=[50, 2000] * u.GeV, energy_power=0)
../_images/notebooks_first_steps_65_0.png

We assign the return axes object to variable called ax_crab_2fhl, 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:

In [33]:

crab_2fhl.spectral_model(100 * u.GeV).to('cm-2 s-1 GeV-1')

Out[33]:
$6.8700477 \times 10^{-12} \; \mathrm{\frac{1}{GeV\,s\,cm^{2}}}$

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

In [34]:

crab_2fhl.spectral_model.integral(
    emin=50 * u.GeV, emax=2000 * u.GeV,
).to('cm-2 s-1')

Out[34]:
$1.31 \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 2FHL catalog:

In [35]:

crab_2fhl.data['Flux50']

Out[35]:

1.31e-09

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

In [36]:

crab_2fhl.spectral_model.energy_flux(
    emin=50 * u.GeV, emax=2000 * u.GeV,
).to('erg cm-2 s-1')

Out[36]:
$3.5295398 \times 10^{-10} \; \mathrm{\frac{erg}{s\,cm^{2}}}$

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

In [37]:

print(crab_2fhl.flux_points)

Flux points of type 'dnde'

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:

In [38]:

crab_2fhl.flux_points.table

Out[38]:
<Table length=3>
e_mine_maxfluxflux_errnflux_errpis_ulflux_ule_refdndednde_errndnde_errpdnde_ul
GeVGeV1 / (cm2 s)1 / (cm2 s)1 / (cm2 s)1 / (cm2 s)TeVGeV(9/8) / (cm2 s TeV(17/8))GeV(9/8) / (cm2 s TeV(17/8))GeV(9/8) / (cm2 s TeV(17/8))GeV(9/8) / (cm2 s TeV(17/8))
float64float64float32float32float32boolfloat64float64float64float64float64float64
50.0171.09.94072e-105.7595e-115.94726e-11Falsenan0.09173269699642.01666069866e-051.16842211486e-061.20651355523e-06nan
171.0585.02.43875e-102.82542e-113.04354e-11Falsenan0.3137728834461.44599896753e-061.6752658163e-071.80459220037e-07nan
585.02000.05.27536e-111.33401e-111.60478e-11Falsenan1.073089512869.15157602841e-082.31421507433e-082.78392984262e-08nan

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:

In [39]:

ax = crab_2fhl.spectral_model.plot(
    energy_range=[50, 2000] * u.GeV, energy_power=2,
)
crab_2fhl.flux_points.plot(ax=ax, energy_power=2)

Out[39]:

<matplotlib.axes._subplots.AxesSubplot at 0x116a8f7f0>
../_images/notebooks_first_steps_79_1.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.