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

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

Gamma-ray source catalogs with Gammapy

Introduction

This notebook shows how to work with gamma-ray source catalogs using gammapy.catalog

  • Load GeV and TeV gamma-ray source catalogs
  • Select a given source and plot it’s spectral model and spectral points
  • Select a subset of sources of interest
  • Plot parameter distributions
  • Plot a sky image with source position markers

You will mainly learn how to work with Astropy tables and how to make some plots to investigate the data.

References:

Setup

In [15]:
%matplotlib inline
import matplotlib.pyplot as plt
In [13]:
from pprint import pprint
import numpy as np
import astropy.units as u
In [2]:
import gammapy
print(gammapy.__version__)
0.6.dev4117

Access catalogs and sources

gamma-cat

In [1]:
from gammapy.catalog import SourceCatalogGammaCat
gammacat = SourceCatalogGammaCat()
gammacat_source = gammacat['Vela Junior']
In [2]:
type(gammacat)
Out[2]:
gammapy.catalog.gammacat.SourceCatalogGammaCat
In [3]:
type(gammacat.table)
Out[3]:
astropy.table.table.QTable
In [4]:
type(gammacat_source)
Out[4]:
gammapy.catalog.gammacat.SourceCatalogObjectGammaCat
In [5]:
type(gammacat_source.data)
Out[5]:
collections.OrderedDict
In [23]:
# You can print the catalog table and source data if you like
# by uncommenting one of the following commands and executing the cell
# gammacat.table.show_in_browser(jsviewer=True)
# gammacat.table.info()
# gammacat.table.info('stats')
# pprint(gammacat_source.data)
# gammacat_source.pprint()

Plot SED

Here’s an example how to plot a TeV SED.

TODO: move content from this section and http://docs.gammapy.org/dev/spectrum/plotting_fermi_spectra.html to a separate notebook about SED access and plotting

In [17]:
# gammacat_source.spectral_model.b
energy_range = [0.3, 40] * u.TeV
gammacat_source.spectral_model.plot(energy_range)
gammacat_source.flux_points.plot()
Out[17]:
<matplotlib.axes._subplots.AxesSubplot at 0x1101899e8>
../_images/notebooks_source_catalogs_14_1.png
In [44]:
from gammapy.catalog import SourceCatalog3FGL
fermi_3fgl = SourceCatalog3FGL()
fermi_3fgl_source = fermi_3fgl['Crab']
In [45]:
type(fermi_3fgl)
Out[45]:
gammapy.catalog.fermi.SourceCatalog3FGL
In [46]:
type(fermi_3fgl.table)
Out[46]:
astropy.table.table.Table
In [47]:
type(fermi_3fgl_source)
Out[47]:
gammapy.catalog.fermi.SourceCatalogObject3FGL
In [48]:
type(fermi_3fgl_source.data)
Out[48]:
collections.OrderedDict

Work in progress …

The rest of this notebook is outdated, it needs to be re-written.

Loading catalogs

The gammapy.datasets module gives you quick access to the most important gamma-ray source catalogs as Astropy table objects.

As you can see we store some parameters in the meta dict of the table that will make it easier to label and color the plots we’ll create below.

In [2]:
from gammapy import datasets
tevcat = datasets.load_tevcat()
tevcat.meta['name'] = 'TeVCat'
tevcat.meta['color'] = 'red'
#hesscat = datasets.load_hess_galactic_catalog()
hesscat = tevcat.copy()
hesscat.meta['name'] = 'HESS Galactic Catalog'
hesscat.meta['color'] = 'green'
fermicat = datasets.fetch_fermi_catalog('2FGL', 'LAT_Point_Source_Catalog')
fermicat.meta['name'] = '2FGL'
fermicat.meta['color'] = 'blue'
WARNING: AstropyDeprecationWarning: Use of an 'eval' method when defining subclasses of FittableModel is deprecated; please rename this method to 'evaluate'.  Otherwise its semantics remain the same. [astropy.modeling.core]
WARNING:astropy:AstropyDeprecationWarning: Use of an 'eval' method when defining subclasses of FittableModel is deprecated; please rename this method to 'evaluate'.  Otherwise its semantics remain the same.
In [3]:
# Show the first three rows.
# In the IPython notebook this automatically gives an HTML table
tevcat[:3]
Out[3]:
canonical_namecatalog_idcatalog_id_namecatalog_namecoord_deccoord_gal_latcoord_gal_loncoord_racoord_typediscovererdiscovery_datedistance_valuedistance_typedistance_kpcethfluxidobservatory_nameother_namesownersize_xsize_ysource_typesource_type_namespec_idxsrc_rankvariability
degdegdegdegkpcTeVdegdeg
PKS 2005-4891Default CatalogTeV J2009-488-48.8311111111-32.605110192350.374122052302.3625-9962005-06-010.071z323761.88689400.00.03100H.E.S.S.HESS J2009-488-99nannan1HBL3.21-99
H 2356-3091Default CatalogTeV J2359-306-30.6229722222-78.041938334312.8562912922359.78925-9962006-04-010.165z800886.991458100.00.02101H.E.S.S.1H 2354-315, HESS J2359-306-99nannan1HBL3.091-99
PG 1553+1131Default CatalogTeV J1555+11111.194722222243.960237045821.9189101955238.93625-9962006-03-010.5z2874104.13541300.00.034102H.E.S.S.1ES 1553+113, HESS J1555+111-99nannan1HBL4.51-99
In [4]:
hesscat[:3]
Out[4]:
canonical_namecatalog_idcatalog_id_namecatalog_namecoord_deccoord_gal_latcoord_gal_loncoord_racoord_typediscovererdiscovery_datedistance_valuedistance_typedistance_kpcethfluxidobservatory_nameother_namesownersize_xsize_ysource_typesource_type_namespec_idxsrc_rankvariability
degdegdegdegkpcTeVdegdeg
PKS 2005-4891Default CatalogTeV J2009-488-48.8311111111-32.605110192350.374122052302.3625-9962005-06-010.071z323761.88689400.00.03100H.E.S.S.HESS J2009-488-99nannan1HBL3.21-99
H 2356-3091Default CatalogTeV J2359-306-30.6229722222-78.041938334312.8562912922359.78925-9962006-04-010.165z800886.991458100.00.02101H.E.S.S.1H 2354-315, HESS J2359-306-99nannan1HBL3.091-99
PG 1553+1131Default CatalogTeV J1555+11111.194722222243.960237045821.9189101955238.93625-9962006-03-010.5z2874104.13541300.00.034102H.E.S.S.1ES 1553+113, HESS J1555+111-99nannan1HBL4.51-99
In [5]:
# Show the first 8 columns and 3 rows
fermicat[fermicat.colnames[:8]][:3]
Out[5]:
Source_NameRAJ2000DEJ2000GLONGLATConf_68_SemiMajorConf_68_SemiMinorConf_68_PosAng
2FGL J0000.9-07480.233711-7.815588.8292-67.28050.1204250.1030447.85
2FGL J0001.7-41590.438849-41.9965334.076-71.99670.07536410.070615362.01
2FGL J0002.7+62200.67981262.3396117.3120.0007522920.05743850.05503698.6

Select a subset of sources of interest

Let’s say we’re interested only interested in sources in the inner Galaxy, namely the box |GLON| < 30 deg and |GLAT| < 5 deg … here’s how you make this selection.

Working with many catalogs is often a bit tedious because column names are different (e.g. GLON and coord_gal_lon) or because of other data cleaning that’s required (like wrapping GLON to the range (-180, 180] here).

In [6]:
def select_sources(catalog):
    """Select sources of interest to us."""
    print('Number of sources before selection: {}'.format(len(catalog)))

    glon = catalog['GLON']
    glat = catalog['GLAT']

    glon_min, glon_max = -30, 30
    glat_min, glat_max = -5, 5

    mask = (glon_min < glon) & (glon < glon_max)
    mask &= (glat_min < glat) & (glat < glat_max)

    catalog = catalog[mask]
    print('Number of sources after selection:  {}'.format(len(catalog)))

    return catalog
In [7]:
tevcat['GLON'] = tevcat['coord_gal_lon']
tevcat['GLAT'] = tevcat['coord_gal_lat']
tevcat['GLON'] = np.where(tevcat['GLON'] > 180, tevcat['GLON'] - 360, tevcat['GLON'])
tevcat2 = select_sources(tevcat)
Number of sources before selection: 173
Number of sources after selection:  45
In [8]:
#hesscat['GLON'] = np.where(hesscat['GLON'] > 180, hesscat['GLON'] - 360, hesscat['GLON'])
#hesscat2 = select_sources(hesscat)
In [9]:
fermicat['GLON'] = np.where(fermicat['GLON'] > 180, fermicat['GLON'] - 360, fermicat['GLON'])
fermicat2 = select_sources(fermicat)
Number of sources before selection: 1873
Number of sources after selection:  126

Plot parameter distributions

Let’s plot the GLON, GLAT parameter distributions.

For the source flux distributions we make separate plots for GeV and TeV sources, because their fluxes aren’t easily comparable.

In [10]:
catalogs = [fermicat2, tevcat2]
In [11]:
for catalog in catalogs:
    plt.hist(catalog['GLON'], bins=10, range=(-30, 30),
             label=catalog.meta['name'],
             color=catalog.meta['color'],
             alpha=0.3, histtype="stepfilled",)
plt.xlabel('Galactic longitude (deg)')
plt.ylabel('Number of sources')
plt.legend();
../_images/notebooks_source_catalogs_36_0.png
In [12]:
for catalog in catalogs:
    plt.hist(catalog['GLAT'], bins=10, range=(-5, 5),
             label=catalog.meta['name'],
             color=catalog.meta['color'],
             alpha=0.3, histtype="stepfilled",)
plt.xlabel('Galactic latitude (deg)')
plt.ylabel('Number of sources')
plt.legend();
../_images/notebooks_source_catalogs_37_0.png
In [13]:
flux = tevcat2['flux']
print(len(flux))
# We need to remove sources with missing flux info (a `NaN` value)
# or matplotlib will not plot it
flux = flux[~np.isnan(flux)]
print(len(flux))
plt.hist(flux, bins=np.logspace(-3, 1, 20),
         color=tevcat2.meta['color'],
         alpha=0.3, histtype="stepfilled",)
plt.xlabel('TeV Flux (% Crab)')
plt.ylabel('Number of sources')
plt.semilogx()
plt.legend();
45
33
/opt/local/Library/Frameworks/Python.framework/Versions/3.4/lib/python3.4/site-packages/matplotlib/axes/_axes.py:475: UserWarning: No labelled objects found. Use label='...' kwarg on individual plots.
  warnings.warn("No labelled objects found. "
../_images/notebooks_source_catalogs_38_2.png

Plot a sky image with source position markers

We use APLPy to plot an empty sky image showing the locations of the sources as markers.

This is a bit complicated because we don’t have an image or WCS to instantiate the aplpy.FITSFigure. Here we’ll use the gammapy.image.make_empty_image utility function to make an empty image to plot on.

In [14]:
from gammapy.image import make_empty_image
from aplpy import FITSFigure
In [15]:
image = make_empty_image(nxpix=62, nypix=12, binsz=1)
figure = FITSFigure(image)
# This call to `show_grayscale` is a workaround to get the Figure stretch correct.
figure.show_grayscale(vmin=-1, vmax=0)
for catalog in catalogs:
    figure.show_markers(catalog['GLON'], catalog['GLAT'],
                        facecolor=catalog.meta['color'],
                        alpha=0.5,
                        label=catalog.meta['name'],
                        )
figure.ticks.set_xspacing(5)
figure.ticks.set_yspacing(2)
figure.tick_labels.set_xformat('dd')
figure.tick_labels.set_yformat('dd')
figure.axis_labels.set_xtext('GLON (deg)')
figure.axis_labels.set_ytext('GLAT (deg)')
plt.legend();
../_images/notebooks_source_catalogs_42_0.png

Exercises

  • Choose another source, e.g. “Crab” or “PKS 2155-304”. Access it for the three catalogs 3fgl, 3fhl and gamma-cat, print the data on it, and make a plot showing the spectral model and spectral points from those catalogs (not necessarily in one plot, but you can try if you like).
  • Make a list of sources in the 3FGL catalog within one degree of the Galactic center, and compute their distance to the Galactic center. Make a plot with source markers and names and separation as labels. Tip: use gc = SkyCoord(0, 0, frame='galactic', unit='deg') to get a SkyCoord object for the Galactic center location.