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: catalog.ipynb | catalog.py
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 0x1c197c6630>
[7]:
from gammapy.catalog import SourceCatalog3FGL
catalog = SourceCatalog3FGL()
catalog
[7]:
<gammapy.catalog.fermi.SourceCatalog3FGL at 0x1c197e75c0>
[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 0x1c1982b898>
[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 0x1c1984dd30>
[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 0x1c198212e8>, spectral_model=<gammapy.modeling.models.spectral.LogParabolaSpectralModel object at 0x1c19469e80>)temporal_model=None)
[19]:
print(model)
SkyModel
Name : 4FGL J2158.8-3013
Spectral model type : LogParabolaSpectralModel
Spatial model type : PointSpatialModel
Temporal model type :
Parameters:
lon_0 : 329.714 deg
lat_0 : -30.225 deg
amplitude : 1.38e-05 1 / (cm2 s TeV)
reference (frozen) : 0.001 TeV
alpha : 1.763
beta : 0.042
[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 0x1c1a339240>
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]:
| e_min | e_max | flux | flux_errn |
|---|---|---|---|
| MeV | MeV | 1 / (cm2 s) | 1 / (cm2 s) |
| float64 | float64 | float32 | float32 |
| 50.0 | 100.0 | 7.72853e-08 | 4.0041947e-08 |
| 100.0 | 300.0 | 6.633264e-08 | 3.3468415e-09 |
| 300.0 | 1000.0 | 3.2824694e-08 | 6.9634254e-10 |
| 1000.0 | 3000.0 | 1.2609229e-08 | 2.4615918e-10 |
| 3000.0 | 10000.0 | 5.407803e-09 | 1.3989285e-10 |
| 10000.0 | 30000.0 | 1.7188525e-09 | 7.4735836e-11 |
| 30000.0 | 300000.0 | 6.818578e-10 | 4.5728827e-11 |
[26]:
flux_points.plot()
[26]:
<matplotlib.axes._subplots.AxesSubplot at 0x1c1b7a1748>
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]:
| time_min | time_max | flux | flux_errp | flux_errn |
|---|---|---|---|---|
| 1 / (cm2 s) | 1 / (cm2 s) | 1 / (cm2 s) | ||
| float64 | float64 | float32 | float32 | float32 |
| 54682.655277777776 | 55047.603239293836 | 8.4104094e-08 | 4.119689e-09 | 4.119689e-09 |
| 55047.603239293836 | 55412.55121238397 | 3.6117186e-08 | 3.434389e-09 | 3.434389e-09 |
| 55412.55121238397 | 55777.4991854741 | 1.4224756e-08 | 2.7470255e-09 | 2.6878597e-09 |
[30]:
lightcurve.plot()
[30]:
<matplotlib.axes._subplots.AxesSubplot at 0x1c1badd550>
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]:
| Source_Name | RAJ2000 | DEJ2000 |
|---|---|---|
| deg | deg | |
| bytes18 | float32 | float32 |
| 3FHL J0001.2-0748 | 0.3107 | -7.8075 |
| 3FHL J0001.9-4155 | 0.4849 | -41.9303 |
| 3FHL J0002.1-6728 | 0.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 degandglat = 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.