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

Source catalogs

gammapy.catalog provides convenient access to common gamma-ray source catalogs. E.g. creating a spectral model and spectral points for a given Fermi-LAT catalog and source from the FITS table is tedious, gammapy.catalog has this implemented and makes it easy.

In this tutorial you will learn how to:

  • List available catalogs

  • Load a catalog

  • Select a source

  • Pretty-print the source information

  • Get source spectral and spatial models

  • Get flux points (if available)

  • Get lightcurves (if available)

  • Access the source catalog table data

In this tutorial we will show examples using the following catalogs:

All catalog and source classes work the same, as long as some information is available. E.g. trying to access a lightcurve from a catalog and source that doesn’t have that information will return None.

Further information is available at gammapy.catalog.

[1]:

%matplotlib inline
import matplotlib.pyplot as plt

[2]:

import astropy.units as u
from gammapy.catalog import SOURCE_CATALOGS

List available catalogs

gammapy.catalog contains a Python dictionary SOURCE_CATALOGS, which maps catalog names (e.g. “3fhl”) to catalog classes (e.g. SourceCatalog3FHL).

[3]:

SOURCE_CATALOGS

[3]:

{'gamma-cat': gammapy.catalog.gammacat.SourceCatalogGammaCat,
 'hgps': gammapy.catalog.hess.SourceCatalogHGPS,
 '2hwc': gammapy.catalog.hawc.SourceCatalog2HWC,
 '3fgl': gammapy.catalog.fermi.SourceCatalog3FGL,
 '4fgl': gammapy.catalog.fermi.SourceCatalog4FGL,
 '2fhl': gammapy.catalog.fermi.SourceCatalog2FHL,
 '3fhl': gammapy.catalog.fermi.SourceCatalog3FHL}

[4]:

list(SOURCE_CATALOGS)

[4]:

['gamma-cat', 'hgps', '2hwc', '3fgl', '4fgl', '2fhl', '3fhl']

Load catalogs

If you have run gammapy download datasets or gammapy download tutorials, you have a copy of the catalogs as FITS files in $GAMMAPY_DATA/catalogs, and that is the default location where gammapy.catalog loads from.

You can load a catalog by name via SOURCE_CATALOG[name]() (not the () to instantiate a catalog object from the catalog class - only this will load the catalog and be useful), or by importing the catalog class (e.g. SourceCatalog3FGL) directly. The two ways are equivalent, the result will be the same.

Note that $GAMMAPY_DATA/catalogs is just the default, you could pass a different filename when creating the catalog.

[5]:

!ls -1 $GAMMAPY_DATA/catalogs

2HWC.ecsv
2HWC.yaml
README.rst
fermi
gammacat
hgps_catalog_v1.fits.gz
make_2hwc.py

[6]:

# Catalog object - FITS file is loaded
catalog = SOURCE_CATALOGS["3fgl"]()
catalog

[6]:

<gammapy.catalog.fermi.SourceCatalog3FGL at 0x110466c88>

[7]:

from gammapy.catalog import SourceCatalog3FGL

catalog = SourceCatalog3FGL()
catalog

[7]:

<gammapy.catalog.fermi.SourceCatalog3FGL at 0x1104c5748>

[8]:

# Let's load the source catalogs we will use throughout this tutorial
catalog_gammacat = SOURCE_CATALOGS["gamma-cat"]
catalog_3fhl = SOURCE_CATALOGS["3fhl"]()
catalog_4fgl = SOURCE_CATALOGS["4fgl"]()
catalog_hgps = SOURCE_CATALOGS["hgps"]()

Select a source

To create a source object, index into the catalog using [], passing a catalog table row index (zero-based, first row is [0]), or a source name. If passing a name, catalog table columns with source names and association names (“ASSOC1” in the example below) are searched top to bottom. There is no name resolution web query.

[9]:

source = catalog_4fgl[42]
source

[9]:

<gammapy.catalog.fermi.SourceCatalogObject4FGL at 0x10282d588>

[10]:

source.row_index, source.name

[10]:

(42, '4FGL J0010.8-2154')

[11]:

source = catalog_4fgl["4FGL J0010.8-2154"]
source

[11]:

<gammapy.catalog.fermi.SourceCatalogObject4FGL at 0x10a419c18>

[12]:

source.row_index, source.name

[12]:

(42, '4FGL J0010.8-2154')

[13]:

source.data["ASSOC1"]

[13]:

'PKS 0008-222                '

[14]:

source = catalog_4fgl["PKS 0008-222"]
source.row_index, source.name

[14]:

(42, '4FGL J0010.8-2154')

Pretty-print source information

A source object has a nice string representation that you can print. You can also call source.info() instead and pass an option what information to print.

[15]:

source = catalog_hgps["MSH 15-52"]
print(source)


*** Basic info ***

Catalog row index (zero-based) : 18
Source name          : HESS J1514-591
Analysis reference   : HGPS
Source class         : PWN
Identified object    : MSH 15-52
Gamma-Cat id         : 79


*** Info from map analysis ***

RA                   :  228.499 deg = 15h14m00s
DEC                  :  -59.161 deg = -59d09m41s
GLON                 :  320.315 +/- 0.008 deg
GLAT                 :   -1.188 +/- 0.007 deg
Position Error (68%) : 0.020 deg
Position Error (95%) : 0.033 deg
ROI number           : 13
Spatial model        : 3-Gaussian
Spatial components   : HGPSC 023, HGPSC 024, HGPSC 025
TS                   : 1763.4
sqrt(TS)             : 42.0
Size                 : 0.145 +/- 0.026 (UL: nan) deg
R70                  : 0.215 deg
RSpec                : 0.215 deg
Total model excess   : 3502.8
Excess in RSpec      : 2440.5
Model Excess in RSpec : 2414.3
Background in RSpec  : 1052.5
Livetime             : 41.4 hours
Energy threshold     : 0.61 TeV
Source flux (>1 TeV) : (6.434 +/- 0.211) x 10^-12 cm^-2 s^-1 = (28.47 +/- 0.94) % Crab

Fluxes in RSpec (> 1 TeV):
Map measurement                : 4.552 x 10^-12 cm^-2 s^-1 = 20.14 % Crab
Source model                   : 4.505 x 10^-12 cm^-2 s^-1 = 19.94 % Crab
Other component model          : 0.000 x 10^-12 cm^-2 s^-1 =  0.00 % Crab
Large scale component model    : 0.000 x 10^-12 cm^-2 s^-1 =  0.00 % Crab
Total model                    : 4.505 x 10^-12 cm^-2 s^-1 = 19.94 % Crab
Containment in RSpec                :  70.0 %
Contamination in RSpec              :   0.0 %
Flux correction (RSpec -> Total)    : 142.8 %
Flux correction (Total -> RSpec)    :  70.0 %

*** Info from spectral analysis ***

Livetime             : 13.7 hours
Energy range:        : 0.38 to 61.90 TeV
Background           : 1825.9
Excess               : 2061.1
Spectral model       : ECPL
TS ECPL over PL      : 10.2
Best-fit model flux(> 1 TeV) : (5.720 +/- 0.417) x 10^-12 cm^-2 s^-1  = (25.31 +/- 1.84) % Crab
Best-fit model energy flux(1 to 10 TeV) : (20.779 +/- 1.878) x 10^-12 erg cm^-2 s^-1
Pivot energy         : 1.54 TeV
Flux at pivot energy : (2.579 +/- 0.083) x 10^-12 cm^-2 s^-1 TeV^-1  = (11.41 +/- 0.24) % Crab
PL   Flux(> 1 TeV)   : (5.437 +/- 0.186) x 10^-12 cm^-2 s^-1  = (24.06 +/- 0.82) % Crab
PL   Flux(@ 1 TeV)   : (6.868 +/- 0.241) x 10^-12 cm^-2 s^-1 TeV^-1  = (30.39 +/- 0.69) % Crab
PL   Index           : 2.26 +/- 0.03
ECPL   Flux(@ 1 TeV) : (6.860 +/- 0.252) x 10^-12 cm^-2 s^-1 TeV^-1  = (30.35 +/- 0.72) % Crab
ECPL   Flux(> 1 TeV) : (5.720 +/- 0.417) x 10^-12 cm^-2 s^-1  = (25.31 +/- 1.84) % Crab
ECPL Index           : 2.05 +/- 0.06
ECPL Lambda          : 0.052 +/- 0.014 TeV^-1
ECPL E_cut           : 19.20 +/- 5.01 TeV

*** Flux points info ***

Number of flux points: 6
Flux points table:

e_ref  e_min  e_max        dnde         dnde_errn       dnde_errp        dnde_ul     is_ul
 TeV    TeV    TeV   1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV) 1 / (cm2 s TeV)
------ ------ ------ --------------- --------------- --------------- --------------- -----
 0.562  0.383  0.825       2.439e-11       1.419e-12       1.509e-12       2.732e-11 False
 1.212  0.825  1.778       4.439e-12       2.489e-13       2.654e-13       4.970e-12 False
 2.738  1.778  4.217       7.295e-13       4.788e-14       4.898e-14       8.302e-13 False
 6.190  4.217  9.085       1.305e-13       1.220e-14       1.282e-14       1.571e-13 False
13.991  9.085 21.544       1.994e-14       2.723e-15       2.858e-15       2.588e-14 False
31.623 21.544 46.416       9.474e-16       3.480e-16       4.329e-16       1.919e-15 False

*** Gaussian component info ***

Number of components: 3
Spatial components   : HGPSC 023, HGPSC 024, HGPSC 025

Component HGPSC 023:
GLON                 :  320.303 +/- 0.005 deg
GLAT                 :   -1.124 +/- 0.007 deg
Size                 : 0.057 +/- 0.005 deg
Flux (>1 TeV)        : (2.01 +/- 0.23) x 10^-12 cm^-2 s^-1 = (8.9 +/- 1.0) % Crab

Component HGPSC 024:
GLON                 :  320.298 +/- 0.020 deg
GLAT                 :   -1.168 +/- 0.021 deg
Size                 : 0.206 +/- 0.020 deg
Flux (>1 TeV)        : (2.54 +/- 0.29) x 10^-12 cm^-2 s^-1 = (11.2 +/- 1.3) % Crab

Component HGPSC 025:
GLON                 :  320.351 +/- 0.005 deg
GLAT                 :   -1.284 +/- 0.007 deg
Size                 : 0.055 +/- 0.005 deg
Flux (>1 TeV)        : (1.88 +/- 0.22) x 10^-12 cm^-2 s^-1 = (8.3 +/- 1.0) % Crab


*** Source associations info ***

  Source_Name    Association_Catalog    Association_Name   Separation
                                                              deg
---------------- ------------------- --------------------- ----------
  HESS J1514-591                2FHL    2FHL J1514.0-5915e   0.098903
  HESS J1514-591                3FGL     3FGL J1513.9-5908   0.026914
  HESS J1514-591                3FGL    3FGL J1514.0-5915e   0.094834
  HESS J1514-591                COMP            G320.4-1.2   0.070483
  HESS J1514-591                 PSR              B1509-58   0.026891

*** Source identification info ***

Source_Name: HESS J1514-591
Identified_Object: MSH 15-52
Class: PWN
Evidence: Morphology
Reference: 2005A%26A...435L..17A
Distance_Reference: SNRCat
Distance: 5.199999809265137 kpc
Distance_Min: 3.799999952316284 kpc
Distance_Max: 6.599999904632568 kpc


[16]:

print(source.info("associations"))


*** Source associations info ***

  Source_Name    Association_Catalog    Association_Name   Separation
                                                              deg
---------------- ------------------- --------------------- ----------
  HESS J1514-591                2FHL    2FHL J1514.0-5915e   0.098903
  HESS J1514-591                3FGL     3FGL J1513.9-5908   0.026914
  HESS J1514-591                3FGL    3FGL J1514.0-5915e   0.094834
  HESS J1514-591                COMP            G320.4-1.2   0.070483
  HESS J1514-591                 PSR              B1509-58   0.026891

Source models

The gammapy.catalog.SourceCatalogObject classes have a sky_model() model which creates a gammapy.modeling.models.SkyModel object, with model parameter values and parameter errors from the catalog filled in.

In most cases, the spectral_model() method provides the gammapy.modeling.models.SpectralModel part of the sky model, and the spatial_model() method the gammapy.modeling.models.SpatialModel part individually.

We use the gammapy.catalog.SourceCatalog3FHL for the examples in this section.

[17]:

source = catalog_4fgl["PKS 2155-304"]

[18]:

model = source.sky_model()
model

[18]:

SkyModel(spatial_model=<gammapy.modeling.models.spatial.PointSpatialModel object at 0x111efcef0>, spectral_model=<gammapy.modeling.models.spectral.LogParabolaSpectralModel object at 0x111efc898>)

[19]:

print(model)

SkyModel

   name     value      error        unit         min        max    frozen
--------- ---------- --------- -------------- ---------- --------- ------
    lon_0  3.297e+02 3.735e-03            deg        nan       nan  False
    lat_0 -3.023e+01 3.227e-03            deg -9.000e+01 9.000e+01  False
amplitude  1.377e-05 1.801e-07 cm-2 s-1 TeV-1        nan       nan  False
reference  1.136e-03 0.000e+00            TeV        nan       nan   True
    alpha  1.763e+00 1.332e-02                       nan       nan  False
     beta  4.172e-02 4.680e-03                       nan       nan  False

[20]:

print(model.spatial_model)

PointSpatialModel

 name   value      error   unit    min        max    frozen
----- ---------- --------- ---- ---------- --------- ------
lon_0  3.297e+02 3.735e-03  deg        nan       nan  False
lat_0 -3.023e+01 3.227e-03  deg -9.000e+01 9.000e+01  False

[21]:

print(model.spectral_model)

LogParabolaSpectralModel

   name     value     error        unit      min max frozen
--------- --------- --------- -------------- --- --- ------
amplitude 1.377e-05 1.801e-07 cm-2 s-1 TeV-1 nan nan  False
reference 1.136e-03 0.000e+00            TeV nan nan   True
    alpha 1.763e+00 1.332e-02                nan nan  False
     beta 4.172e-02 4.680e-03                nan nan  False

[22]:

energy_range = (100 * u.MeV, 100 * u.GeV)
opts = dict(energy_power=2, flux_unit="erg-1 cm-2 s-1")
model.spectral_model.plot(energy_range, **opts)
model.spectral_model.plot_error(energy_range, **opts)

[22]:

<matplotlib.axes._subplots.AxesSubplot at 0x110f616a0>
../_images/notebooks_catalog_28_1.png

Flux points

The flux points are available via the flux_points property as a gammapy.spectrum.FluxPoints object.

[23]:

source = catalog_4fgl["PKS 2155-304"]
flux_points = source.flux_points

[24]:

flux_points

[24]:

FluxPoints(sed_type='flux', n_points=7)

[25]:

flux_points.table[["e_min", "e_max", "flux", "flux_errn"]]

[25]:
Table length=7
e_mine_maxfluxflux_errn
MeVMeV1 / (cm2 s)1 / (cm2 s)
float64float64float32float32
50.0100.07.72853e-084.0041947e-08
100.0300.06.633264e-083.3468415e-09
300.01000.03.2824694e-086.9634254e-10
1000.03000.01.2609229e-082.4615918e-10
3000.010000.05.407803e-091.3989285e-10
10000.030000.01.7188525e-097.4735836e-11
30000.0300000.06.818578e-104.5728827e-11
[26]:

flux_points.plot()

[26]:

<matplotlib.axes._subplots.AxesSubplot at 0x1123b9128>
../_images/notebooks_catalog_33_1.png

Lightcurves

The Fermi catalogs contain lightcurves for each source. It is available via the source.lightcurve property as a gammapy.time.LightCurve object.

[27]:

lightcurve = catalog_4fgl["4FGL J0349.8-2103"].lightcurve

[28]:

lightcurve

[28]:

LightCurve(len=8)

[29]:

lightcurve.table[:3]

[29]:
Table length=3
time_mintime_maxfluxflux_errpflux_errn
1 / (cm2 s)1 / (cm2 s)1 / (cm2 s)
float64float64float32float32float32
54682.65527777777655047.6032392938368.4104094e-084.119689e-094.119689e-09
55047.60323929383655412.551212383973.6117186e-083.434389e-093.434389e-09
55412.5512123839755777.49918547411.4224756e-082.7470255e-092.6878597e-09
[30]:

lightcurve.plot()

[30]:

<matplotlib.axes._subplots.AxesSubplot at 0x11202a550>
../_images/notebooks_catalog_38_1.png

Catalog table and source dictionary

Source catalogs are given as FITS files that contain one or multiple tables. Above we showed how to get spectra, light curves and other information as Gammapy objects.

However, you can also access the underlying ~astropy.table.Table for a catalog, and the row data as a Python dict. This can be useful if you want to do something that is not pre-scripted by the gammapy.catalog classes, such as e.g. selecting sources by sky position or association class, or accessing special source information (like e.g. Npred in the example below).

Note that you can also do a for source in catalog loop, to find or process sources of interest.

[31]:

type(catalog_3fhl.table)

[31]:

astropy.table.table.Table

[32]:

len(catalog_3fhl.table)

[32]:

1556

[33]:

catalog_3fhl.table[:3][["Source_Name", "RAJ2000", "DEJ2000"]]

[33]:
Table masked=True length=3
Source_NameRAJ2000DEJ2000
degdeg
bytes18float32float32
3FHL J0001.2-07480.3107-7.8075
3FHL J0001.9-41550.4849-41.9303
3FHL J0002.1-67280.5283-67.4825
[34]:

source = catalog_3fhl["PKS 2155-304"]

[35]:

source.data["Source_Name"]

[35]:

'3FHL J2158.8-3013 '

[36]:

source.data["Npred"]

[36]:

678.95544

[37]:

source.position

[37]:

<SkyCoord (ICRS): (ra, dec) in deg
    (329.71694946, -30.22418594)>

[38]:

# Find the brightest sources in the 100 to 200 GeV energy band
for source in catalog_3fhl:
    flux = (
        source.spectral_model()
        .integral(100 * u.GeV, 200 * u.GeV)
        .to("cm-2 s-1")
    )
    if flux > 1e-10 * u.Unit("cm-2 s-1"):
        print(f"{source.row_index:<7d} {source.name:20s} {flux:.3g}")

352     3FHL J0534.5+2201    2.99e-10 1 / (cm2 s)
553     3FHL J0851.9-4620e   1.24e-10 1 / (cm2 s)
654     3FHL J1036.3-5833e   1.57e-10 1 / (cm2 s)
691     3FHL J1104.4+3812    3.34e-10 1 / (cm2 s)
1111    3FHL J1653.8+3945    1.27e-10 1 / (cm2 s)
1219    3FHL J1824.5-1351e   1.77e-10 1 / (cm2 s)
1361    3FHL J2028.6+4110e   1.75e-10 1 / (cm2 s)

Exercises

  • How many sources are in the 4FGL catalog? (try len(catalog.table)

  • What is the name of the source with row index 42?

  • What is the row index of the source with name “4FGL J0536.1-1205”?

  • What is the integral flux of “4FGL J0536.1-1205” in the energy range 100 GeV to 1 TeV according to the best-fit spectral model?

  • Which source in the HGPS catalog is closest to Galactic position glon = 42 deg and glat = 0 deg?

[39]:

# Start coding here ...

Next steps

gammapy.catalog is mostly independent from the rest of Gammapy. Typically you use it to compare new analyses against catalog results, e.g. overplot the spectral model, or compare the source position.

You can also use gammapy.catalog in your scripts to create initial source models for your analyses. This is very common for Fermi-LAT, to start with a catalog model. For TeV analysis, especially in crowded Galactic regions, using the HGPS, gamma-cat or 2HWC catalog in this way can also be useful.