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

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]:

fromgammapy.imageimportSkyImage

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

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:

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()

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:

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 regioncenter=SkyCoord(265.0,-2.0,unit='deg',frame='galactic')size=9.0*u.degvela_2fhl_cutout=vela_2fhl_smoothed.cutout(center,size)vela_2fhl_cutout.show()

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 visiblenorm=simple_norm(vela_wmap.data,stretch='sqrt',max_percent=99.9)vela_wmap.show(cmap='viridis',norm=norm)

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.

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 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 centerevents_vela_2fhl=events_2fhl.select_sky_cone(center=center,radius=0.5*u.deg)# sort events by energyevents_vela_2fhl.table.sort('ENERGY')# and show highest energy photonevents_vela_2fhl.energy[-1].to('GeV')

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 TSfermi_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 classtop_five_TS_2fhl[['Source_Name','ASSOC','CLASS']]

Out[29]:

<Table length=5>

Source_Name

ASSOC

CLASS

str18

str25

str8

2FHL J1104.4+3812

Mkn 421

bll

2FHL J0534.5+2201

Crab

pwn

2FHL J1653.9+3945

Mkn 501

bll

2FHL J1555.7+1111

PG 1553+113

bll

2FHL J2158.8-3013

PKS 2155-304

bll

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 catalogmkn_421_2fhl=fermi_2fhl['Mkn 421']print(mkn_421_2fhl.data['TS'])

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

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 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.