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

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:

Setup

In [1]:
# to make plots appear in the notebook
%matplotlib inline
import matplotlib.pyplot as plt
In [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 gammapy.extern.pathlib import Path
import numpy as np
import astropy
import os

print("numpy:", np.__version__)
print("astropy:", astropy.__version__)
numpy: 1.13.0
astropy: 3.0.4
In [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

In [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)
In [5]:
# And convert its value to an equivalent unit
flux.to("erg cm-2 s-1")
Out[5]:
$$1.3058974 \times 10^{-11} \; \mathrm{\frac{erg}{s\,cm^{2}}}$$
In [6]:
flux.to("W/m2")
Out[6]:
$$1.3058974 \times 10^{-14} \; \mathrm{\frac{W}{m^{2}}}$$

More generally a Quantity is a numpy array with a unit.

In [7]:
E = np.logspace(1, 4, 10) * u.GeV
E.to("TeV")
Out[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.

In [8]:
x_eff = 30 * u.mbarn
density = 1 * u.cm ** -3

interaction_time = (density * x_eff * cst.c) ** -1

interaction_time.to("Myr")
Out[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

In [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. / 3. * 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
In [10]:
# Now plot it
E_elec = np.logspace(-1., 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);
../_images/notebooks_astropy_introduction_15_0.png

A frequent issue is homogeneity. One can use decorators to ensure it.

In [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. / 3. * 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.

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

In [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?

In [14]:
now = Time.now()
print(now)
print(now.mjd)
2018-09-24 10:38:14.353572
58385.44322168486
In [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=2018-09-24 10:38:14.353572, location=(4200910.643257838, 172456.78503911156, 4780088.658775934) m, pressure=0.0 hPa, temperature=0.0 deg_C, relative_humidity=0.0, obswl=1.0 micron): (az, alt) in deg
    ( 278.69436218,  22.05936223)>

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.

In [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]
In [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

In [18]:
table.info()
<Table length=3034>
        name          dtype   shape        unit       n_bad
-------------------- ------- ------- ---------------- -----
         Source_Name bytes18                              0
             RAJ2000 float32                      deg     0
             DEJ2000 float32                      deg     0
                GLON float32                      deg     0
                GLAT float32                      deg     0
   Conf_68_SemiMajor float32                      deg    28
   Conf_68_SemiMinor float32                      deg    28
      Conf_68_PosAng float32                      deg    28
   Conf_95_SemiMajor float32                      deg    28
   Conf_95_SemiMinor float32                      deg    28
      Conf_95_PosAng float32                      deg    28
             ROI_num   int16                              0
          Signif_Avg float32                              5
        Pivot_Energy float32                      MeV     0
        Flux_Density float32         ph / (cm2 MeV s)     0
    Unc_Flux_Density float32         ph / (cm2 MeV s)     5
            Flux1000 float32             ph / (cm2 s)     0
        Unc_Flux1000 float32             ph / (cm2 s)     5
      Energy_Flux100 float32            erg / (cm2 s)     0
  Unc_Energy_Flux100 float32            erg / (cm2 s)     5
        Signif_Curve float32                              0
        SpectrumType bytes16                              0
      Spectral_Index float32                              0
  Unc_Spectral_Index float32                              9
                beta float32                           2639
            Unc_beta float32                           2680
              Cutoff float32                      MeV  2918
          Unc_Cutoff float32                      MeV  2918
           Exp_Index float32                           2918
       Unc_Exp_Index float32                           3028
      PowerLaw_Index float32                              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                           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                              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                              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                              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                              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                              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
In [19]:
# Statistics on each column
table.info("stats")
<Table length=3034>
        name              mean          std          min         max     n_bad
-------------------- ------------- ------------- ----------- ----------- -----
         Source_Name            --            --          --          --     0
             RAJ2000       185.638       100.585      0.0377     359.881     0
             DEJ2000      -2.23514       41.5371    -87.6185     88.7739     0
                GLON       186.273       107.572    0.126161     359.969     0
                GLAT       1.22728       36.4998    -87.9516     86.3679     0
   Conf_68_SemiMajor          -inf           nan        -inf    0.658726    28
   Conf_68_SemiMinor          -inf           nan        -inf    0.320945    28
      Conf_68_PosAng     -0.219671       51.9263      -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.219671       51.9263      -89.94       89.92    28
             ROI_num 424.777191826 238.846231996           1         840     0
          Signif_Avg          -inf           nan        -inf     1048.96     5
        Pivot_Energy       2104.46        2578.2     100.798     37202.0     0
        Flux_Density   1.29389e-11   2.17231e-10 4.21515e-16 1.14476e-08     0
    Unc_Flux_Density          -inf           nan        -inf 6.86476e-10     5
            Flux1000    2.7507e-09   2.76475e-08 4.90582e-12 1.29768e-06     0
        Unc_Flux1000          -inf           nan        -inf 2.89396e-09     5
      Energy_Flux100   2.56011e-11   1.88859e-10 1.65227e-12 8.93008e-09     0
  Unc_Energy_Flux100          -inf           nan        -inf 2.52953e-11     5
        Signif_Curve       2.43161       3.87731         0.0      84.973     0
        SpectrumType            --            --          --          --     0
      Spectral_Index       2.21451      0.360539         0.5     5.71589     0
  Unc_Spectral_Index          -inf           nan        -inf     1.06113     9
                beta          -inf           nan        -inf         1.0  2639
            Unc_beta          -inf           nan        -inf    0.337608  2680
              Cutoff          -inf           nan        -inf     13326.4  2918
          Unc_Cutoff          -inf           nan        -inf     3750.56  2918
           Exp_Index          -inf           nan        -inf         1.0  2918
       Unc_Exp_Index          -inf           nan        -inf   0.0442999  3028
      PowerLaw_Index       2.24982      0.306354     1.09753     5.71589     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.37795e-08   1.20342e-07 1.00222e-14 5.37625e-06     0
     Unc_Flux100_300          -inf           nan        -inf 9.56701e-08   836
        nuFnu100_300   5.71922e-12   3.02685e-11 2.49123e-18 1.37177e-09     0
      Sqrt_TS100_300       3.81654       8.91583         0.0     215.416     5
        Flux300_1000   8.61936e-09   6.73293e-08 3.97788e-15 3.21672e-06     0
    Unc_Flux300_1000          -inf           nan        -inf 1.89716e-08   327
       nuFnu300_1000   5.90484e-12   4.79845e-11 2.91179e-18 2.29243e-09     0
     Sqrt_TS300_1000        7.4718       19.8981         0.0     568.663     5
       Flux1000_3000   2.19502e-09   2.27795e-08 1.50877e-15 1.07019e-06     0
   Unc_Flux1000_3000          -inf           nan        -inf 6.09948e-09   111
      nuFnu1000_3000   5.14613e-12   5.38784e-11 3.76523e-18 2.52534e-09     0
    Sqrt_TS1000_3000       8.60778       22.2479         0.0     711.632     0
      Flux3000_10000   4.83629e-10   4.64785e-09 9.33897e-16 2.16177e-07     0
  Unc_Flux3000_10000          -inf           nan        -inf 2.98606e-09   137
     nuFnu3000_10000   3.07335e-12   2.83225e-11  5.6769e-18 1.32102e-09     0
   Sqrt_TS3000_10000       6.98671       13.9194         0.0     465.899     3
    Flux10000_100000   7.39921e-11   3.22538e-10 1.03878e-17 1.13815e-08     0
Unc_Flux10000_100000          -inf           nan        -inf 7.88092e-10   670
   nuFnu10000_100000   1.16414e-12   4.63575e-12 9.38528e-20 1.16018e-10     0
 Sqrt_TS10000_100000       4.01938       5.47004         0.0     111.192     2
   Variability_Index       154.126       1230.71         0.0     60733.9     0
         Signif_Peak          -inf           nan        -inf     255.071  2387
           Flux_Peak          -inf           nan        -inf 2.12374e-05  2387
       Unc_Flux_Peak          -inf           nan        -inf 1.47859e-07  2387
           Time_Peak          -inf           nan        -inf 364203776.0  2387
       Peak_Interval          -inf           nan        -inf    2.63e+06  2387
        Flux_History   3.54401e-08   2.26052e-07         0.0 2.12374e-05     0
    Unc_Flux_History          -inf           nan        -inf 2.30386e-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.7735662492 265.711605004           0        2565     0
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
/Users/jer/anaconda/envs/gammapy-dev/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:1308: RuntimeWarning: invalid value encountered in subtract
  np.subtract(arr, avg, out=arr, casting='unsafe')
In [20]:
### list of column names
table.colnames
Out[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']
In [21]:
# HTML display
# table.show_in_browser(jsviewer=True)
# table.show_in_notebook(jsviewer=True)

Accessing the table

In [22]:
# The header keywords are stored as a dict
# table.meta
table.meta["TSMIN"]
Out[22]:
25
In [23]:
# First row
table[0]
Out[23]:
Row index=0
Source_NameRAJ2000DEJ2000GLONGLATConf_68_SemiMajorConf_68_SemiMinorConf_68_PosAngConf_95_SemiMajorConf_95_SemiMinorConf_95_PosAngROI_numSignif_AvgPivot_EnergyFlux_DensityUnc_Flux_DensityFlux1000Unc_Flux1000Energy_Flux100Unc_Energy_Flux100Signif_CurveSpectrumTypeSpectral_IndexUnc_Spectral_IndexbetaUnc_betaCutoffUnc_CutoffExp_IndexUnc_Exp_IndexPowerLaw_IndexFlux30_100Unc_Flux30_100nuFnu30_100Sqrt_TS30_100Flux100_300Unc_Flux100_300 [2]nuFnu100_300Sqrt_TS100_300Flux300_1000Unc_Flux300_1000 [2]nuFnu300_1000Sqrt_TS300_1000Flux1000_3000Unc_Flux1000_3000 [2]nuFnu1000_3000Sqrt_TS1000_3000Flux3000_10000Unc_Flux3000_10000 [2]nuFnu3000_10000Sqrt_TS3000_10000Flux10000_100000Unc_Flux10000_100000 [2]nuFnu10000_100000Sqrt_TS10000_100000Variability_IndexSignif_PeakFlux_PeakUnc_Flux_PeakTime_PeakPeak_IntervalFlux_History [48]Unc_Flux_History [48,2]Extended_Source_Name0FGL_Name1FGL_Name2FGL_Name1FHL_NameASSOC_GAM1ASSOC_GAM2ASSOC_GAM3TEVCAT_FLAGASSOC_TEVCLASS1ASSOC1ASSOC2Flags
degdegdegdegdegdegdegdegdegdegMeVph / (cm2 MeV s)ph / (cm2 MeV s)ph / (cm2 s)ph / (cm2 s)erg / (cm2 s)erg / (cm2 s)MeVMeVph / (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)ssph / (cm2 s)ph / (cm2 s)
bytes18float32float32float32float32float32float32float32float32float32float32int16float32float32float32float32float32float32float32float32float32bytes16float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float32float64float32float32float32bytes18bytes17bytes18bytes18bytes18bytes15bytes14bytes15bytes1bytes21bytes5bytes26bytes26int16
3FGL J0000.1+65450.037765.7517117.6943.402960.06284450.048104741.030.10190.07841.031856.813131159.081.00534e-121.49106e-131.01596e-091.57592e-101.35695e-112.07317e-123.39662PowerLaw2.410960.0823208-inf-inf-inf-inf-inf-inf2.41096nannannannan1.80832e-08-8.39548e-09 .. 8.23604e-094.17529e-122.168756.94073e-09-1.35281e-09 .. 1.36692e-094.54365e-125.269711.23665e-09-2.25128e-10 .. 2.34347e-102.85534e-126.022245.78052e-11-4.04457e-11 .. 4.8534e-113.78413e-131.509192.83963e-11-1.48428e-11 .. 1.91547e-114.3172e-132.4211740.7535-inf-inf-inf-inf-inf6.36157e-08 .. 0.0-2.7777e-08 .. 2.09131e-082FGL J2359.6+6543cN4
In [24]:
# Spectral index of the 5 first entries
table[:5]["Spectral_Index"]
Out[24]:
<Column name='Spectral_Index' dtype='float32' length=5>
2.41096
1.86673
2.73279
2.14935
2.30511
In [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.5
Softest source:  3FGL J0534.5+2201s PWN   5.71589

Quantities and SkyCoords from a Table

In [26]:
fluxes = table["nuFnu1000_3000"].quantity
fluxes
Out[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}}}$$
In [27]:
coord = SkyCoord(table["GLON"], table["GLAT"], frame="galactic")
coord.fk5
Out[27]:
<SkyCoord (FK5: equinox=J2000.000): (ra, dec) in deg
    [(  3.87171473e-02,  65.75181023), (  6.14202362e-02, -37.64829309),
     (  2.54455001e-01,  63.2441096 ), ...,
     (  3.59736655e+02,  39.44771268), (  3.59829367e+02, -30.64638846),
     (  3.59881529e+02, -20.88248885)]>

Selections in a Table

Here we select Sources according to their class and do some whole sky chart

In [28]:
# Get coordinates of FSRQs
fsrq = np.where(
    np.logical_or(table["CLASS1"] == "fsrq ", table["CLASS1"] == "FSQR ")
)
In [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...
Out[29]:
<matplotlib.legend.Legend at 0x116ae7e10>
../_images/notebooks_astropy_introduction_43_1.png
In [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();
../_images/notebooks_astropy_introduction_44_0.png

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.

In [31]:
rows = [dict(a=42, b="spam"), dict(a=43, b="ham")]
my_table = Table(rows=rows)
my_table
Out[31]:
Table length=2
ab
int64str4
42spam
43ham

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.

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

In [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):

In [34]:
scalar_colnames = tuple(
    name for name in table.colnames if len(table[name].shape) <= 1
)
data_frame = table[scalar_colnames].to_pandas()
In [35]:
# If you want to have a quick-look at the dataframe:
# data_frame
# data_frame.info()
# data_frame.describe()
In [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()
Out[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'agn  '       3
b'RDG  '       3
b'HMB  '       3
b'nlsy1'       3
b'ssrq '       3
b'NLSY1'       2
b'pwn  '       2
b'GAL  '       2
b'css  '       1
b'NOV  '       1
b'gal  '       1
b'sey  '       1
b'SFR  '       1
b'BIN  '       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