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: astropy_introduction.ipynb | astropy_introduction.py
Astropy introduction for Gammapy users¶
Introduction¶
To become efficient at using Gammapy, you have to learn to use some parts of Astropy, especially FITS I/O and how to work with Table, Quantity, SkyCoord, Angle and Time objects.
Gammapy is built on Astropy, meaning that data in Gammapy is often stored in Astropy objects, and the methods on those objects are part of the public Gammapy API.
This tutorial is a quick introduction to the parts of Astropy you should know become familiar with to use Gammapy (or also when not using Gammapy, just doing astronomy from Python scripts). The largest part is devoted to tables, which are a the most important building block for Gammapy (event lists, flux points, light curves, … many other thing are store in Table objects).
We will:
- open and write fits files with io.fits
- manipulate coordinates: SkyCoord and Angle classes
- use units and Quantities. See also this tutorial
- manipulate Times and Dates
- use tables with the Table class with the Fermi catalog
- define regions in the sky with the region package
Setup¶
[1]:
# to make plots appear in the notebook
%matplotlib inline
import matplotlib.pyplot as plt
[2]:
# If something below doesn't work, here's how you
# can check what version of Numpy and Astropy you have
# All examples should work with Astropy 1.3,
# most even with Astropy 1.0
from pathlib import Path
import numpy as np
import astropy
import os
print("numpy:", np.__version__)
print("astropy:", astropy.__version__)
numpy: 1.16.2
astropy: 3.1.1
[3]:
# Units, Quantities and constants
import astropy.units as u
from astropy.units import Quantity
import astropy.constants as cst
from astropy.io import fits
from astropy.table import Table
from astropy.coordinates import SkyCoord, Angle
from astropy.time import Time
Units and constants¶
Basic usage¶
[4]:
# One can create a Quantity like this
L = Quantity(1e35, unit="erg/s")
# or like this
d = 8 * u.kpc
# then one can produce new Quantities
flux = L / (4 * np.pi * d ** 2)
[5]:
# And convert its value to an equivalent unit
flux.to("erg cm-2 s-1")
[5]:
$$1.3058974 \times 10^{-11} \; \mathrm{\frac{erg}{s\,cm^{2}}}$$
[6]:
flux.to("W/m2")
[6]:
$$1.3058974 \times 10^{-14} \; \mathrm{\frac{W}{m^{2}}}$$
More generally a Quantity is a numpy array with a unit.
[7]:
E = np.logspace(1, 4, 10) * u.GeV
E.to("TeV")
[7]:
$$[0.01,~0.021544347,~0.046415888,~0.1,~0.21544347,~0.46415888,~1,~2.1544347,~4.6415888,~10] \; \mathrm{TeV}$$
Here we compute the interaction time of protons.
[8]:
x_eff = 30 * u.mbarn
density = 1 * u.cm ** -3
interaction_time = (density * x_eff * cst.c) ** -1
interaction_time.to("Myr")
[8]:
$$35.233361 \; \mathrm{Myr}$$
Use Quantities in functions¶
We compute here the energy loss rate of an electron of kinetic energy E in magnetic field B. See formula (5B10) in this lecture
[9]:
def electron_energy_loss_rate(B, E):
""" energy loss rate of an electron of kinetic energy E in magnetic field B
"""
U_B = B ** 2 / (2 * cst.mu0)
gamma = (
E / (cst.m_e * cst.c ** 2) + 1
) # note that this works only because E/(cst.m_e*cst.c**2) is dimensionless
beta = np.sqrt(1 - 1 / gamma ** 2)
return 4.0 / 3.0 * cst.sigma_T * cst.c * gamma ** 2 * beta ** 2 * U_B
print(electron_energy_loss_rate(1e-5 * u.G, 1 * u.TeV).to("erg/s"))
4.0519325060487625e-13 erg / s
[10]:
# Now plot it
E_elec = np.logspace(-1.0, 6, 100) * u.MeV
B = 1 * u.G
y = (E_elec / electron_energy_loss_rate(B, E_elec)).to("yr")
plt.loglog(E_elec, y);
A frequent issue is homogeneity. One can use decorators to ensure it.
[11]:
# This ensures that B and E are homogeneous to magnetic field strength and energy
# If not will raise a UnitError exception
@u.quantity_input(B=u.T, E=u.J)
def electron_energy_loss_rate(B, E):
""" energy loss rate of an electron of kinetic energy E in magnetic field B
"""
U_B = B ** 2 / (2 * cst.mu0)
gamma = (
E / (cst.m_e * cst.c ** 2) + 1
) # note that this works only because E/(cst.m_e*cst.c**2) is dimensionless
beta = np.sqrt(1 - 1 / gamma ** 2)
return 4.0 / 3.0 * cst.sigma_T * cst.c * gamma ** 2 * beta ** 2 * U_B
# Now try it
try:
print(electron_energy_loss_rate(1e-5 * u.G, 1 * u.Hz).to("erg/s"))
except u.UnitsError as message:
print("Incorrect unit: " + str(message))
Incorrect unit: Argument 'E' to function 'electron_energy_loss_rate' must be in units convertible to 'J'.
Coordinates¶
Note that SkyCoord are arrays of coordinates. We will see that in more detail in the next section.
[12]:
# Different ways to create a SkyCoord
c1 = SkyCoord(10.625, 41.2, frame="icrs", unit="deg")
c1 = SkyCoord("00h42m30s", "+41d12m00s", frame="icrs")
c2 = SkyCoord(83.633083, 22.0145, unit="deg")
# If you have internet access, you could also use this to define the `source_pos`:
# c2 = SkyCoord.from_name("Crab") # Get the name from CDS
print(c1.ra, c2.dec)
# separation returns an Angle object
print("Distance to Crab: ", c1.separation(c2))
print("Distance to Crab: ", c1.separation(c2).degree)
10d37m30s 22d00m52.2s
Distance to Crab: 63d12m28.2297s
Distance to Crab: 63.2078415848386
Coordinate transformations¶
How to change between coordinate frames. The Crab in Galactic coordinates.
[13]:
c2b = c2.galactic
print(c2b)
print(c2b.l, c2b.b)
<SkyCoord (Galactic): (l, b) in deg
(184.55745771, -5.78435696)>
184d33m26.8478s -5d47m03.6851s
Time¶
Is the Crab visible now?
[14]:
now = Time.now()
print(now)
print(now.mjd)
2019-03-29 16:41:51.788169
58571.695738288996
[15]:
# define the location for the AltAz system
from astropy.coordinates import EarthLocation, AltAz
paris = EarthLocation(lat=48.8567 * u.deg, lon=2.3508 * u.deg)
# calculate the horizontal coordinates
crab_altaz = c2.transform_to(AltAz(obstime=now, location=paris))
print(crab_altaz)
<SkyCoord (AltAz: obstime=2019-03-29 16:41:51.788169, location=(4200910.64325784, 172456.78503911, 4780088.65877593) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
(171.19416159, 62.95038106)>
Table: Manipulating the 3FGL catalog¶
Here we are going to do some selections with the 3FGL catalog. To do so we use the Table class from astropy.
Accessing the table¶
First, we need to open the catalog in a Table.
[16]:
# Open Fermi 3FGL from the repo
filename = os.environ["GAMMAPY_DATA"] / Path(
"catalogs/fermi/gll_psc_v16.fit.gz"
)
table = Table.read(str(filename), hdu=1)
# Alternatively, one can grab it from the server.
# table = Table.read("http://fermi.gsfc.nasa.gov/ssc/data/access/lat/4yr_catalog/gll_psc_v16.fit")
WARNING: UnitsWarning: 'photon/cm**2/MeV/s' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'photon/cm**2/s' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
WARNING: UnitsWarning: 'erg/cm**2/s' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]
[17]:
# Note that a single FITS file might contain different tables in different HDUs
filename = os.environ["GAMMAPY_DATA"] / Path(
"catalogs/fermi/gll_psc_v16.fit.gz"
)
# You can load a `fits.HDUList` and check the extension names
print([_.name for _ in fits.open(str(filename))])
# Then you can load by name or integer index via the `hdu` option
extended_source_table = Table.read(str(filename), hdu="ExtendedSources")
['PRIMARY', 'LAT_Point_Source_Catalog', 'ROIs', 'Hist_Start', 'GTI', 'ExtendedSources']
General informations on the Table¶
[18]:
table.info()
<Table length=3034>
name dtype shape unit format n_bad
-------------------- ------- ------- ---------------- -------- -----
Source_Name bytes18 0
RAJ2000 float32 deg {:8.4f} 0
DEJ2000 float32 deg {:8.4f} 0
GLON float32 deg {:8.4f} 0
GLAT float32 deg {:8.4f} 0
Conf_68_SemiMajor float32 deg {:8.4f} 28
Conf_68_SemiMinor float32 deg {:8.4f} 28
Conf_68_PosAng float32 deg {:8.3f} 28
Conf_95_SemiMajor float32 deg {:8.4f} 28
Conf_95_SemiMinor float32 deg {:8.4f} 28
Conf_95_PosAng float32 deg {:8.3f} 28
ROI_num int16 0
Signif_Avg float32 {:8.3f} 5
Pivot_Energy float32 MeV {:10.2f} 0
Flux_Density float32 ph / (cm2 MeV s) {:10.4e} 0
Unc_Flux_Density float32 ph / (cm2 MeV s) {:10.4e} 5
Flux1000 float32 ph / (cm2 s) {:10.4e} 0
Unc_Flux1000 float32 ph / (cm2 s) 5
Energy_Flux100 float32 erg / (cm2 s) {:10.4e} 0
Unc_Energy_Flux100 float32 erg / (cm2 s) {:10.4e} 5
Signif_Curve float32 {:10.2f} 0
SpectrumType bytes16 0
Spectral_Index float32 {:8.4f} 0
Unc_Spectral_Index float32 {:8.4f} 9
beta float32 {:8.4f} 2639
Unc_beta float32 {:8.4f} 2680
Cutoff float32 MeV {:10.2f} 2918
Unc_Cutoff float32 MeV {:10.2f} 2918
Exp_Index float32 {:8.4f} 2918
Unc_Exp_Index float32 {:8.4f} 3028
PowerLaw_Index float32 {:8.4f} 0
Flux30_100 float32 ph / (cm2 s) 3034
Unc_Flux30_100 float32 ph / (cm2 s) 3034
nuFnu30_100 float32 erg / (cm2 s) 3034
Sqrt_TS30_100 float32 {:8.3f} 3034
Flux100_300 float32 ph / (cm2 s) 0
Unc_Flux100_300 float32 (2,) ph / (cm2 s) 836
nuFnu100_300 float32 erg / (cm2 s) 0
Sqrt_TS100_300 float32 {:8.3f} 5
Flux300_1000 float32 ph / (cm2 s) 0
Unc_Flux300_1000 float32 (2,) ph / (cm2 s) 327
nuFnu300_1000 float32 erg / (cm2 s) 0
Sqrt_TS300_1000 float32 {:8.3f} 5
Flux1000_3000 float32 ph / (cm2 s) 0
Unc_Flux1000_3000 float32 (2,) ph / (cm2 s) 111
nuFnu1000_3000 float32 erg / (cm2 s) 0
Sqrt_TS1000_3000 float32 {:8.3f} 0
Flux3000_10000 float32 ph / (cm2 s) 0
Unc_Flux3000_10000 float32 (2,) ph / (cm2 s) 137
nuFnu3000_10000 float32 erg / (cm2 s) 0
Sqrt_TS3000_10000 float32 {:8.3f} 3
Flux10000_100000 float32 ph / (cm2 s) 0
Unc_Flux10000_100000 float32 (2,) ph / (cm2 s) 670
nuFnu10000_100000 float32 erg / (cm2 s) 0
Sqrt_TS10000_100000 float32 {:8.3f} 2
Variability_Index float32 0
Signif_Peak float32 2387
Flux_Peak float32 ph / (cm2 s) 2387
Unc_Flux_Peak float32 ph / (cm2 s) 2387
Time_Peak float64 s 2387
Peak_Interval float32 s 2387
Flux_History float32 (48,) ph / (cm2 s) 0
Unc_Flux_History float32 (48, 2) ph / (cm2 s) 64623
Extended_Source_Name bytes18 0
0FGL_Name bytes17 0
1FGL_Name bytes18 0
2FGL_Name bytes18 0
1FHL_Name bytes18 0
ASSOC_GAM1 bytes15 0
ASSOC_GAM2 bytes14 0
ASSOC_GAM3 bytes15 0
TEVCAT_FLAG bytes1 0
ASSOC_TEV bytes21 0
CLASS1 bytes5 0
ASSOC1 bytes26 0
ASSOC2 bytes26 0
Flags int16 0
[19]:
# Statistics on each column
table.info("stats")
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
<Table length=3034>
name mean std min max n_bad
-------------------- ----------------- ------------------ ------------- -------------- -----
Source_Name -- -- -- -- 0
RAJ2000 185.63785 100.58523 0.0377 359.8812 0
DEJ2000 -2.235143 41.53712 -87.6185 88.7739 0
GLON 186.2728 107.57225 0.12616104 359.9692 0
GLAT 1.2272784 36.499836 -87.95155 86.36791 0
Conf_68_SemiMajor -inf nan -inf 0.6587264 28
Conf_68_SemiMinor -inf nan -inf 0.32094488 28
Conf_68_PosAng -0.21967074 51.926346 -89.94 89.92 28
Conf_95_SemiMajor -inf nan -inf 1.0681 28
Conf_95_SemiMinor -inf nan -inf 0.5204 28
Conf_95_PosAng -0.21967074 51.926346 -89.94 89.92 28
ROI_num 424.7771918259723 238.84623199551552 1 840 0
Signif_Avg -inf nan -inf 1048.9587 5
Pivot_Energy 2104.4626 2578.198 100.79754 37201.97 0
Flux_Density 1.2938883e-11 2.1723061e-10 4.2151523e-16 1.1447582e-08 0
Unc_Flux_Density -inf nan -inf 6.864756e-10 5
Flux1000 2.7506988e-09 2.764754e-08 4.905819e-12 1.2976768e-06 0
Unc_Flux1000 -inf nan -inf 2.8939553e-09 5
Energy_Flux100 2.5601129e-11 1.8885869e-10 1.6522738e-12 8.9300825e-09 0
Unc_Energy_Flux100 -inf nan -inf 2.5295337e-11 5
Signif_Curve 2.431611 3.8773124 0.0 84.97297 0
SpectrumType -- -- -- -- 0
Spectral_Index 2.2145069 0.36053884 0.50000006 5.715894 0
Unc_Spectral_Index -inf nan -inf 1.0611331 9
beta -inf nan -inf 1.0 2639
Unc_beta -inf nan -inf 0.3376085 2680
Cutoff -inf nan -inf 13326.353 2918
Unc_Cutoff -inf nan -inf 3750.5642 2918
Exp_Index -inf nan -inf 1.0 2918
Unc_Exp_Index -inf nan -inf 0.044299893 3028
PowerLaw_Index 2.249819 0.30635425 1.0975342 5.715894 0
Flux30_100 nan nan nan nan 3034
Unc_Flux30_100 nan nan nan nan 3034
nuFnu30_100 nan nan nan nan 3034
Sqrt_TS30_100 nan nan nan nan 3034
Flux100_300 2.3779466e-08 1.2034207e-07 1.0022202e-14 5.37625e-06 0
Unc_Flux100_300 -inf nan -inf 9.567011e-08 836
nuFnu100_300 5.7192215e-12 3.026846e-11 2.4912323e-18 1.3717709e-09 0
Sqrt_TS100_300 3.8165362 8.915827 0.0 215.4162 5
Flux300_1000 8.619365e-09 6.732932e-08 3.9778793e-15 3.2167186e-06 0
Unc_Flux300_1000 -inf nan -inf 1.897159e-08 327
nuFnu300_1000 5.904836e-12 4.7984526e-11 2.9117903e-18 2.2924256e-09 0
Sqrt_TS300_1000 7.471796 19.89813 0.0 568.66296 5
Flux1000_3000 2.1950213e-09 2.2779465e-08 1.5087699e-15 1.0701934e-06 0
Unc_Flux1000_3000 -inf nan -inf 6.0994774e-09 111
nuFnu1000_3000 5.1461344e-12 5.3878357e-11 3.765228e-18 2.5253446e-09 0
Sqrt_TS1000_3000 8.607782 22.247915 0.0 711.6319 0
Flux3000_10000 4.836294e-10 4.6478523e-09 9.338969e-16 2.1617714e-07 0
Unc_Flux3000_10000 -inf nan -inf 2.9860603e-09 137
nuFnu3000_10000 3.0733539e-12 2.8322539e-11 5.6769025e-18 1.3210241e-09 0
Sqrt_TS3000_10000 6.98671 13.919392 0.0 465.89935 3
Flux10000_100000 7.399214e-11 3.2253794e-10 1.0387837e-17 1.13814576e-08 0
Unc_Flux10000_100000 -inf nan -inf 7.880918e-10 670
nuFnu10000_100000 1.1641434e-12 4.635748e-12 9.385276e-20 1.1601801e-10 0
Sqrt_TS10000_100000 4.019383 5.470039 0.0 111.192406 2
Variability_Index 154.12613 1230.7148 0.0 60733.92 0
Signif_Peak -inf nan -inf 255.07051 2387
Flux_Peak -inf nan -inf 2.1237443e-05 2387
Unc_Flux_Peak -inf nan -inf 1.4785873e-07 2387
Time_Peak -inf nan -inf 364203776.0 2387
Peak_Interval -inf nan -inf 2630000.0 2387
Flux_History 3.5440095e-08 2.2605154e-07 0.0 2.1237443e-05 0
Unc_Flux_History -inf nan -inf 2.3038561e-07 64623
Extended_Source_Name -- -- -- -- 0
0FGL_Name -- -- -- -- 0
1FGL_Name -- -- -- -- 0
2FGL_Name -- -- -- -- 0
1FHL_Name -- -- -- -- 0
ASSOC_GAM1 -- -- -- -- 0
ASSOC_GAM2 -- -- -- -- 0
ASSOC_GAM3 -- -- -- -- 0
TEVCAT_FLAG -- -- -- -- 0
ASSOC_TEV -- -- -- -- 0
CLASS1 -- -- -- -- 0
ASSOC1 -- -- -- -- 0
ASSOC2 -- -- -- -- 0
Flags 45.77356624917601 265.7116050036644 0 2565 0
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/adonath/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/numpy/lib/nanfunctions.py:1503: RuntimeWarning: invalid value encountered in subtract
np.subtract(arr, avg, out=arr, casting='unsafe')
[20]:
### list of column names
table.colnames
[20]:
['Source_Name',
'RAJ2000',
'DEJ2000',
'GLON',
'GLAT',
'Conf_68_SemiMajor',
'Conf_68_SemiMinor',
'Conf_68_PosAng',
'Conf_95_SemiMajor',
'Conf_95_SemiMinor',
'Conf_95_PosAng',
'ROI_num',
'Signif_Avg',
'Pivot_Energy',
'Flux_Density',
'Unc_Flux_Density',
'Flux1000',
'Unc_Flux1000',
'Energy_Flux100',
'Unc_Energy_Flux100',
'Signif_Curve',
'SpectrumType',
'Spectral_Index',
'Unc_Spectral_Index',
'beta',
'Unc_beta',
'Cutoff',
'Unc_Cutoff',
'Exp_Index',
'Unc_Exp_Index',
'PowerLaw_Index',
'Flux30_100',
'Unc_Flux30_100',
'nuFnu30_100',
'Sqrt_TS30_100',
'Flux100_300',
'Unc_Flux100_300',
'nuFnu100_300',
'Sqrt_TS100_300',
'Flux300_1000',
'Unc_Flux300_1000',
'nuFnu300_1000',
'Sqrt_TS300_1000',
'Flux1000_3000',
'Unc_Flux1000_3000',
'nuFnu1000_3000',
'Sqrt_TS1000_3000',
'Flux3000_10000',
'Unc_Flux3000_10000',
'nuFnu3000_10000',
'Sqrt_TS3000_10000',
'Flux10000_100000',
'Unc_Flux10000_100000',
'nuFnu10000_100000',
'Sqrt_TS10000_100000',
'Variability_Index',
'Signif_Peak',
'Flux_Peak',
'Unc_Flux_Peak',
'Time_Peak',
'Peak_Interval',
'Flux_History',
'Unc_Flux_History',
'Extended_Source_Name',
'0FGL_Name',
'1FGL_Name',
'2FGL_Name',
'1FHL_Name',
'ASSOC_GAM1',
'ASSOC_GAM2',
'ASSOC_GAM3',
'TEVCAT_FLAG',
'ASSOC_TEV',
'CLASS1',
'ASSOC1',
'ASSOC2',
'Flags']
[21]:
# HTML display
# table.show_in_browser(jsviewer=True)
# table.show_in_notebook(jsviewer=True)
Accessing the table¶
[22]:
# The header keywords are stored as a dict
# table.meta
table.meta["TSMIN"]
[22]:
25
[23]:
# First row
table[0]
[23]:
Row index=0
| Source_Name | RAJ2000 | DEJ2000 | GLON | GLAT | Conf_68_SemiMajor | Conf_68_SemiMinor | Conf_68_PosAng | Conf_95_SemiMajor | Conf_95_SemiMinor | Conf_95_PosAng | ROI_num | Signif_Avg | Pivot_Energy | Flux_Density | Unc_Flux_Density | Flux1000 | Unc_Flux1000 | Energy_Flux100 | Unc_Energy_Flux100 | Signif_Curve | SpectrumType | Spectral_Index | Unc_Spectral_Index | beta | Unc_beta | Cutoff | Unc_Cutoff | Exp_Index | Unc_Exp_Index | PowerLaw_Index | Flux30_100 | Unc_Flux30_100 | nuFnu30_100 | Sqrt_TS30_100 | Flux100_300 | Unc_Flux100_300 [2] | nuFnu100_300 | Sqrt_TS100_300 | Flux300_1000 | Unc_Flux300_1000 [2] | nuFnu300_1000 | Sqrt_TS300_1000 | Flux1000_3000 | Unc_Flux1000_3000 [2] | nuFnu1000_3000 | Sqrt_TS1000_3000 | Flux3000_10000 | Unc_Flux3000_10000 [2] | nuFnu3000_10000 | Sqrt_TS3000_10000 | Flux10000_100000 | Unc_Flux10000_100000 [2] | nuFnu10000_100000 | Sqrt_TS10000_100000 | Variability_Index | Signif_Peak | Flux_Peak | Unc_Flux_Peak | Time_Peak | Peak_Interval | Flux_History [48] | Unc_Flux_History [48,2] | Extended_Source_Name | 0FGL_Name | 1FGL_Name | 2FGL_Name | 1FHL_Name | ASSOC_GAM1 | ASSOC_GAM2 | ASSOC_GAM3 | TEVCAT_FLAG | ASSOC_TEV | CLASS1 | ASSOC1 | ASSOC2 | Flags |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| deg | deg | deg | deg | deg | deg | deg | deg | deg | deg | MeV | ph / (cm2 MeV s) | ph / (cm2 MeV s) | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | erg / (cm2 s) | MeV | MeV | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | ph / (cm2 s) | ph / (cm2 s) | erg / (cm2 s) | ph / (cm2 s) | ph / (cm2 s) | s | s | ph / (cm2 s) | ph / (cm2 s) | ||||||||||||||||||||||||||||||||||
| bytes18 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | int16 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | bytes16 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float32 | float64 | float32 | float32 | float32 | bytes18 | bytes17 | bytes18 | bytes18 | bytes18 | bytes15 | bytes14 | bytes15 | bytes1 | bytes21 | bytes5 | bytes26 | bytes26 | int16 |
| 3FGL J0000.1+6545 | 0.0377 | 65.7517 | 117.6939 | 3.4030 | 0.0628 | 0.0481 | 41.030 | 0.1019 | 0.0780 | 41.030 | 185 | 6.813 | 1159.08 | 1.0053e-12 | 1.4911e-13 | 1.0160e-09 | 1.5759212e-10 | 1.3570e-11 | 2.0732e-12 | 3.40 | PowerLaw | 2.4110 | 0.0823 | -inf | -inf | -inf | -inf | -inf | -inf | 2.4110 | nan | nan | nan | nan | 1.8083163e-08 | -8.39548e-09 .. 8.236045e-09 | 4.1752856e-12 | 2.169 | 6.9407333e-09 | -1.3528098e-09 .. 1.3669199e-09 | 4.543647e-12 | 5.270 | 1.2366493e-09 | -2.2512775e-10 .. 2.3434676e-10 | 2.8553431e-12 | 6.022 | 5.78052e-11 | -4.04457e-11 .. 4.8533958e-11 | 3.7841312e-13 | 1.509 | 2.8396294e-11 | -1.484277e-11 .. 1.91547e-11 | 4.3172041e-13 | 2.421 | 40.753517 | -inf | -inf | -inf | -inf | -inf | 6.361567e-08 .. 0.0 | -2.777704e-08 .. 2.0913118e-08 | 2FGL J2359.6+6543c | N | 4 |
[24]:
# Spectral index of the 5 first entries
table[:5]["Spectral_Index"]
[24]:
<Column name='Spectral_Index' dtype='float32' format='{:8.4f}' length=5>
| 2.4110 |
| 1.8667 |
| 2.7328 |
| 2.1493 |
| 2.3051 |
[25]:
# Which source has the lowest spectral index?
row = table[np.argmin(table["Spectral_Index"])]
print(
"Hardest source: ",
row["Source_Name"],
row["CLASS1"],
row["Spectral_Index"],
)
# Which source has the largest spectral index?
row = table[np.argmax(table["Spectral_Index"])]
print(
"Softest source: ",
row["Source_Name"],
row["CLASS1"],
row["Spectral_Index"],
)
Hardest source: 3FGL J2051.3-0828 PSR 0.50000006
Softest source: 3FGL J0534.5+2201s PWN 5.715894
Quantities and SkyCoords from a Table¶
[26]:
fluxes = table["nuFnu1000_3000"].quantity
fluxes
[26]:
$$[2.8553431 \times 10^{-12},~3.5327822 \times 10^{-13},~6.7680339 \times 10^{-13},~\dots,~7.3150578 \times 10^{-13},~1.0483972 \times 10^{-12},~4.7821024 \times 10^{-13}] \; \mathrm{\frac{erg}{s\,cm^{2}}}$$
[27]:
coord = SkyCoord(table["GLON"], table["GLAT"], frame="galactic")
coord.fk5
[27]:
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
[(3.87169046e-02, 65.75181035), (6.14199915e-02, -37.64829374),
(2.54454515e-01, 63.2441095 ), ..., (3.59736654e+02, 39.44771292),
(3.59829367e+02, -30.64638844), (3.59881529e+02, -20.88248953)]>
Selections in a Table¶
Here we select Sources according to their class and do some whole sky chart
[28]:
# Get coordinates of FSRQs
fsrq = np.where(
np.logical_or(table["CLASS1"] == "fsrq ", table["CLASS1"] == "FSQR ")
)
[29]:
# This is here for plotting purpose...
# glon = glon.wrap_at(180*u.degree)
# Open figure
fig = plt.figure(figsize=(14, 8))
ax = fig.add_subplot(111, projection="aitoff")
ax.scatter(
coord[fsrq].l.wrap_at(180 * u.degree).radian,
coord[fsrq].b.radian,
color="k",
label="FSRQ",
)
ax.grid(True)
ax.legend()
# ax.invert_xaxis() -> This does not work for projections...
[29]:
<matplotlib.legend.Legend at 0x112a83390>
[30]:
# Now do it for a series of classes
fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(111, projection="aitoff")
source_classes = ["", "psr", "spp", "fsrq", "bll", "bin"]
for source_class in source_classes:
# We select elements with correct class in upper or lower characters
index = np.array(
[_.strip().lower() == source_class for _ in table["CLASS1"]]
)
label = source_class if source_class else "unid"
ax.scatter(
coord[index].l.wrap_at(180 * u.degree).radian,
coord[index].b.radian,
label=label,
)
ax.grid(True)
ax.legend();
Creating tables¶
A Table is basically a dict mapping column names to column values, where a column value is a Numpy array (or Quantity object, which is a Numpy array sub-class). This implies that adding columns to a table after creation is nice and easy, but adding a row is hard and slow, basically all data has to be copied and all objects that make up a Table have to be re-created.
Here’s one way to create a Table from scratch: put the data into a list of dicts, and then call the Table constructor with the rows option.
[31]:
rows = [dict(a=42, b="spam"), dict(a=43, b="ham")]
my_table = Table(rows=rows)
my_table
[31]:
Table length=2
| a | b |
|---|---|
| int64 | str4 |
| 42 | spam |
| 43 | ham |
Writing tables¶
Writing tables to files is easy, you can just give the filename and format you want. If you run a script repeatedly you might want to add overwrite=True.
[32]:
# Examples how to write a table in different formats
# Uncomment if you really want to do it
# my_table.write('/tmp/table.fits', format='fits')
# my_table.write('/tmp/table.fits.gz', format='fits')
# my_table.write('/tmp/table.ecsv', format='ascii.ecsv')
FITS (and some other formats, e.g. HDF5) support writing multiple tables to a single file. The table.write API doesn’t support that directly yet. Here’s how you can currently write multiple tables to a FITS file: you have to convert the astropy.table.Table objects to astropy.io.fits.BinTable objects, and then store them in a astropy.io.fits.HDUList objects and call HDUList.writeto.
[33]:
my_table2 = Table(data=dict(a=[1, 2, 3]))
hdu_list = fits.HDUList(
[
fits.PrimaryHDU(), # need an empty primary HDU
fits.table_to_hdu(my_table),
fits.table_to_hdu(my_table2),
]
)
# hdu_list.writeto('tables.fits')
Tables and pandas¶
pandas is one of the most-used packages in the scientific Python stack. Numpy provides the ndarray object and functions that operate on ndarray objects. Pandas provides the Dataframe and Series objects, which roughly correspond to the Astropy Table and Column objects. While both pandas.Dataframe and astropy.table.Table can often be used to work with tabular data, each has features that the other doesn’t. When Astropy was started, it
was decided to not base it on pandas.Dataframe, but to introduce Table, mainly because pandas.Dataframe doesn’t support multi-dimensional columns, but FITS does and astronomers use sometimes.
But pandas.Dataframe has a ton of features that Table doesn’t, and is highly optimised, so if you find something to be hard with Table, you can convert it to a Dataframe and do your work there. As explained in the interfacing with the pandas package page in the Astropy docs, it is easy to go back and forth between Table and Dataframe:
table = Table.from_pandas(dataframe)
dataframe = table.to_pandas()
Let’s try it out with the Fermi-LAT catalog.
One little trick is needed when converting to a dataframe: we need to drop the multi-dimensional columns that the 3FGL catalog uses for a few columns (flux up/down errors, and lightcurves):
[34]:
scalar_colnames = tuple(
name for name in table.colnames if len(table[name].shape) <= 1
)
data_frame = table[scalar_colnames].to_pandas()
[35]:
# If you want to have a quick-look at the dataframe:
# data_frame
# data_frame.info()
# data_frame.describe()
[36]:
# Just do demonstrate one of the useful DataFrame methods,
# this is how you can count the number of sources in each class:
data_frame["CLASS1"].value_counts()
[36]:
b' ' 1010
b'bll ' 642
b'bcu ' 568
b'fsrq ' 446
b'PSR ' 143
b'spp ' 49
b'FSRQ ' 38
b'psr ' 24
b'BLL ' 18
b'glc ' 15
b'SNR ' 12
b'rdg ' 12
b'snr ' 11
b'PWN ' 10
b'BCU ' 5
b'sbg ' 4
b'ssrq ' 3
b'HMB ' 3
b'RDG ' 3
b'nlsy1' 3
b'agn ' 3
b'NLSY1' 2
b'GAL ' 2
b'pwn ' 2
b'NOV ' 1
b'SFR ' 1
b'css ' 1
b'BIN ' 1
b'sey ' 1
b'gal ' 1
Name: CLASS1, dtype: int64
If you’d like to learn more about pandas, have a look here or here.
Exercises¶
- When searched for the hardest and softest sources in 3FGL we did not look at the type of spectrum (PL, ECPL etc), find the hardest and softest PL sources instead.
- Replot the full sky chart of sources in ra-dec instead of galactic coordinates
- Find the 3FGL sources visible from Paris now
[37]: