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

Flux Profile Estimation#

This tutorial shows how to estimate flux profiles.

Prerequisites#

Knowledge of 3D data reduction and datasets used in Gammapy, see for instance the first analysis tutorial.

Context#

A useful tool to study and compare the saptial distribution of flux in images and data cubes is the measurement of flxu profiles. Flux profiles can show spatial correlations of gamma-ray data with e.g. gas maps or other type of gamma-ray data. Most commonly flux profiles are measured along some preferred coordinate axis, either radially distance from a source of interest, along longitude and latitude coordinate axes or along the path defined by two spatial coordinates.

Proposed Approach#

Flux profile estimation essentially works by estimating flux points for a set of predefined spatially connected regions. For radial flux profiles the shape of the regions are annuli with a common center, for linear profiles it’s typically a rectangular shape.

We will work on a pre-computed MapDataset of Fermi-LAT data, use Region to define the structure of the bins of the flux profile and run the actually profile extraction using the FluxProfileEstimator

Setup and Imports#

[1]:
%matplotlib inline
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.coordinates import SkyCoord
import numpy as np
[2]:
from gammapy.datasets import MapDataset
from gammapy.estimators import FluxProfileEstimator, FluxPoints
from gammapy.maps import RegionGeom
from gammapy.utils.regions import (
    make_concentric_annulus_sky_regions,
    make_orthogonal_rectangle_sky_regions,
)
from gammapy.modeling.models import PowerLawSpectralModel

Read and Introduce Data#

[3]:
dataset = MapDataset.read(
    "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz", name="fermi-dataset"
)

This is what the counts image we will work with looks like:

[4]:
counts_image = dataset.counts.sum_over_axes()
counts_image.smooth("0.1 deg").plot(stretch="sqrt");
../../../_images/tutorials_analysis_3D_flux_profiles_8_0.png

There are 400x200 pixels in the dataset and 11 energy bins between 10 GeV and 2 TeV:

[5]:
print(dataset.counts)
WcsNDMap

        geom  : WcsGeom
        axes  : ['lon', 'lat', 'energy']
        shape : (400, 200, 11)
        ndim  : 3
        unit  :
        dtype : >i4

Profile Estimation#

Configuration#

We start by defining a list of spatially connected regions along the galactic longitude axis. For this there is a helper function make_orthogonal_rectangle_sky_regions. The individual region bins for the profile have a height of 3 deg and in total there are 31 bins. The starts from lon = 10 deg tand goes to lon = 350 deg. In addition we have to specify the wcs to take into account possible projections effects on the region definition:

[6]:
regions = make_orthogonal_rectangle_sky_regions(
    start_pos=SkyCoord("10d", "0d", frame="galactic"),
    end_pos=SkyCoord("350d", "0d", frame="galactic"),
    wcs=counts_image.geom.wcs,
    height="3 deg",
    nbin=51,
)

We can use the RegionGeom object to illustrate the regions on top of the counts image:

[7]:
geom = RegionGeom.create(region=regions)
ax = counts_image.smooth("0.1 deg").plot(stretch="sqrt")
geom.plot_region(ax=ax, color="w");
/Users/terrier/Code/anaconda3/envs/gammapy-dev/lib/python3.8/site-packages/regions/shapes/rectangle.py:206: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties.
  return Rectangle(xy=xy, width=width, height=height,
../../../_images/tutorials_analysis_3D_flux_profiles_14_1.png

Next we create the FluxProfileEstimator. For the estimation of the flux profile we assume a spectral model with a power-law shape and an index of 2.3

[8]:
flux_profile_estimator = FluxProfileEstimator(
    regions=regions,
    spectrum=PowerLawSpectralModel(index=2.3),
    energy_edges=[10, 2000] * u.GeV,
    selection_optional=["ul"],
)

We can see the full configuration by printing the estimator object:

[9]:
print(flux_profile_estimator)
FluxProfileEstimator
--------------------

  energy_edges           : [  10. 2000.] GeV
  fit                    : <gammapy.modeling.fit.Fit object at 0x146e7c070>
  n_sigma                : 1
  n_sigma_ul             : 2
  norm_max               : 5
  norm_min               : 0.2
  norm_n_values          : 11
  norm_values            : None
  null_value             : 0
  reoptimize             : False
  selection_optional     : ['ul']
  source                 : 0
  spectrum               : PowerLawSpectralModel
  sum_over_energy_groups : False

Run Estimation#

Now we can run the profile estimation and explore the results:

[10]:
%%time
profile = flux_profile_estimator.run(datasets=dataset)
CPU times: user 48 s, sys: 721 ms, total: 48.7 s
Wall time: 1min
[11]:
print(profile)
FluxPoints
----------

  geom                   : RegionGeom
  axes                   : ['lon', 'lat', 'energy', 'projected-distance']
  shape                  : (1, 1, 1, 51)
  quantities             : ['norm', 'norm_err', 'norm_ul', 'ts', 'npred', 'npred_excess', 'stat', 'counts', 'success']
  ref. model             : pl
  n_sigma                : 1
  n_sigma_ul             : 2
  sqrt_ts_threshold_ul   : 2
  sed type init          : likelihood

We can see the flux profile is represented by a FluxPoints object with a projected-distance axis, which defines the main axis the flux profile is measured along. The lon and lat axes can be ignored.

Plotting Results#

Let us directly plot the result using FluxPoints.plot():

[12]:
ax = profile.plot(sed_type="dnde")
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_23_0.png

Based on the spectral model we specified above we can also plot in any other sed type, e.g. energy flux and define a different threshold when to plot upper limits:

[13]:
profile.sqrt_ts_threshold_ul = 2

ax = profile.plot(sed_type="eflux")
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_25_0.png

We can also plot any other quantity of interest, that is defined on the FluxPoints result object. E.g. the predicted total counts, background counts and excess counts:

[14]:
quantities = ["npred", "npred_excess", "npred_background"]

for quantity in quantities:
    profile[quantity].plot(label=quantity.title())

plt.ylabel("Counts ")
[14]:
Text(0, 0.5, 'Counts ')
../../../_images/tutorials_analysis_3D_flux_profiles_27_1.png

Serialisation and I/O#

The profile can be serialised using FluxPoints.write(), given a specific format:

[15]:
profile.write(
    filename="flux_profile_fermi.fits",
    format="profile",
    overwrite=True,
    sed_type="dnde",
)
[16]:
profile_new = FluxPoints.read(
    filename="flux_profile_fermi.fits", format="profile"
)
No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.
[17]:
ax = profile_new.plot()
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_31_0.png

The profile can be serialised to a ~astropy.table.Table object using:

[18]:
table = profile.to_table(format="profile", formatted=True)
table
[18]:
Table length=51
x_minx_maxx_refe_ref [1]e_min [1]e_max [1]ref_dnde [1]ref_flux [1]ref_eflux [1]norm [1]norm_err [1]norm_ul [1]ts [1]sqrt_ts [1]npred [1,1]npred_excess [1,1]stat [1]is_ul [1]counts [1,1]success [1]
degdegdegkeVkeVkeV1 / (cm2 s TeV)keV / (cm2 s TeV)keV2 / (cm2 s TeV)
float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float32float64boolfloat64bool
-0.19607843137254920.19607843137254920.070710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
0.19607843137254920.58823529411764670.39215686274509870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
0.58823529411764670.98039215686274430.784313725490195570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.1720.1240.4342.0591.435162.2304898611596717.448376-818.848True163.0True
0.98039215686274431.37254901960784361.17647058823529470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+061.8940.2072.322120.59110.981450.20181889217974192.60909-3013.446False448.0True
1.37254901960784361.7647058823529421.568627450980392870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+063.1460.2403.639280.64816.753601.6822045013222320.15994-4348.966False599.0True
1.7647058823529422.156862745098041.96078431372549170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7360.1831.11618.9914.358357.665030488639674.95968-2218.623False354.0True
2.156862745098042.5490196078431382.352941176470588870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6560.1861.04214.2143.770363.8333813660937666.89659-2452.717False367.0True
2.5490196078431382.9411764705882362.74509803921568770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4500.1790.8217.0112.648339.481827331467645.943058-2174.005False339.0True
2.9411764705882363.33333333333333443.13725490196078570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+062.4210.2252.885174.14613.196528.0599574489612247.276-3887.531False531.0True
3.33333333333333443.7254901960784323.52941176470588370710678.11910000000.000500000000.0004.428e-103.043e-019.166e+062.1900.2092.621169.69813.027457.3561643539827223.84384-3183.973False458.0True
3.7254901960784324.1176470588235293.921568627450980770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+062.8680.2323.345243.82115.615564.4775740614491293.2505-4201.048False566.0True
4.1176470588235294.5098039215686274.31372549019607870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+061.5990.2052.02381.9019.050451.7627320041883163.63295-3040.782False448.0True
4.5098039215686274.9019607843137274.705882352941177570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.5500.1820.92810.2433.201356.2851725096595656.336945-2313.852False356.0True
4.9019607843137275.2941176470588245.09803921568627670710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.2100.1710.5661.5861.259315.4980561641578421.543198-1979.035True315.0True
............................................................
13.92156862745100214.31372549019607914.1176470588235470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3770.1870.7634.4072.099385.156676485136139.272778-2526.278False384.0True
14.31372549019607914.70588235294118114.5098039215686370710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06-0.2440.1680.1032.022-1.422313.3736443463295-25.44988-1936.261True313.0True
14.70588235294118115.09803921568630814.90196078431374470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3070.1790.6793.1331.770366.32254744581432.01158-2178.012True357.0True
15.09803921568630815.49019607843138515.29411764705884770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06-0.0830.1670.2650.240-0.490313.4101919680402-8.618814-1914.901True311.0True
15.49019607843138515.88235294117646415.68627450980392470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.2280.1720.5841.8631.365330.8911319400612623.80037-2016.119True327.0True
15.88235294117646416.27450980392159316.07843137254902670710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4270.1700.7806.9792.642324.6901970807428544.554733-1957.576False321.0True
16.27450980392159316.66666666666669316.47058823529414470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7600.1951.16317.5534.190428.00520659718179.45124-2767.621False421.0True
16.66666666666669317.0588235294117716.86274509803923470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6750.1971.08213.2773.644437.1341506313758670.588745-2830.574False431.0True
17.0588235294117717.450980392156917.25490196078433870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.2560.1830.6342.0571.434381.2475096218782426.727816-2356.744True374.0True
17.450980392156917.84313725490200417.64705882352945270710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.1740.1820.5510.9510.975366.4533491178318.242867-2463.364True370.0True
17.84313725490200418.2352941176470818.0392156862745470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.5930.1920.99010.7053.272408.0550496449774362.12828-2803.387False410.0True
18.2352941176470818.62745098039215518.43137254901961770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3060.1730.6663.3421.828334.1414142466950432.074398-2186.572True336.0True
18.62745098039215519.0196078431372618.82352941176470770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.0100.1240.2690.0050.071171.81156491395930.9959075-877.032True172.0True
19.0196078431372619.41176470588238319.2156862745098270710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
19.41176470588238319.80392156862749419.6078431372549470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False

No we can also estimate a radial profile starting from the Galactic center:

[19]:
regions = make_concentric_annulus_sky_regions(
    center=SkyCoord("0d", "0d", frame="galactic"),
    radius_max="1.5 deg",
    nbin=11,
)

Again we first illustrate the regions:

[20]:
geom = RegionGeom.create(region=regions)
gc_image = counts_image.cutout(
    position=SkyCoord("0d", "0d", frame="galactic"), width=3 * u.deg
)
ax = gc_image.smooth("0.1 deg").plot(stretch="sqrt")
geom.plot_region(ax=ax, color="w");
/Users/terrier/Code/anaconda3/envs/gammapy-dev/lib/python3.8/site-packages/regions/core/compound.py:160: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties.
  patch = mpatches.PathPatch(path, **mpl_kwargs)
../../../_images/tutorials_analysis_3D_flux_profiles_37_1.png

This time we define two energy bins and include the fit statistic profile in the computation:

[21]:
flux_profile_estimator = FluxProfileEstimator(
    regions=regions,
    spectrum=PowerLawSpectralModel(index=2.3),
    energy_edges=[10, 100, 2000] * u.GeV,
    selection_optional=["ul", "scan"],
    norm_values=np.linspace(-1, 5, 11),
)
[22]:
profile = flux_profile_estimator.run(datasets=dataset)

We can directly plot the result:

[23]:
ax = profile.plot(axis_name="projected-distance", sed_type="flux")
../../../_images/tutorials_analysis_3D_flux_profiles_42_0.png

However because of the powerlaw spectrum the flux at high energies is much lower. To extract the profile at high energies only we can use:

[24]:
profile_high = profile.slice_by_idx({"energy": slice(1, 2)})

And now plot the points together with the likelihood profiles:

[25]:
ax = profile_high.plot(sed_type="eflux", color="tab:orange")
profile_high.plot_ts_profiles(ax=ax, sed_type="eflux")
ax.set_yscale("linear")
../../../_images/tutorials_analysis_3D_flux_profiles_46_0.png
[ ]: