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

• Try online

• You may download all the notebooks as a tar file.

• Source files: flux_profiles.ipynb | flux_profiles.py

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

%matplotlib inline
import matplotlib.pyplot as plt
from astropy import units as u
from astropy.coordinates import SkyCoord
import numpy as np

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

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:

counts_image = dataset.counts.sum_over_axes()
counts_image.smooth("0.1 deg").plot(stretch="sqrt");

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

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:

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:

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,

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

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:

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:

%%time
profile = flux_profile_estimator.run(datasets=dataset)

CPU times: user 48 s, sys: 721 ms, total: 48.7 s
Wall time: 1min

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():

ax = profile.plot(sed_type="dnde")
ax.set_yscale("linear")

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:

profile.sqrt_ts_threshold_ul = 2

ax = profile.plot(sed_type="eflux")
ax.set_yscale("linear")

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:

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

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

plt.ylabel("Counts ")

Text(0, 0.5, 'Counts ')

Serialisation and I/O

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

profile.write(
    filename="flux_profile_fermi.fits",
    format="profile",
    overwrite=True,
    sed_type="dnde",
)

profile_new = FluxPoints.read(
    filename="flux_profile_fermi.fits", format="profile"
)

No reference model set for FluxMaps. Assuming point source with E^-2 spectrum.

ax = profile_new.plot()
ax.set_yscale("linear")

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

table = profile.to_table(format="profile", formatted=True)
table

Table length=51

x_min	x_max	x_ref	e_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]
deg	deg	deg	keV	keV	keV	1 / (cm2 s TeV)	keV / (cm2 s TeV)	keV2 / (cm2 s TeV)
float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float32	float64	bool	float64	bool
-0.1960784313725492	0.1960784313725492	0.0	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	nan	nan	nan	nan	nan	nan	0.0	0.000	False	0.0	False
0.1960784313725492	0.5882352941176467	0.392156862745098	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	nan	nan	nan	nan	nan	nan	0.0	0.000	False	0.0	False
0.5882352941176467	0.9803921568627443	0.7843137254901955	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.172	0.124	0.434	2.059	1.435	162.23048986115967	17.448376	-818.848	True	163.0	True
0.9803921568627443	1.3725490196078436	1.176470588235294	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	1.894	0.207	2.322	120.591	10.981	450.20181889217974	192.60909	-3013.446	False	448.0	True
1.3725490196078436	1.764705882352942	1.5686274509803928	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	3.146	0.240	3.639	280.648	16.753	601.6822045013222	320.15994	-4348.966	False	599.0	True
1.764705882352942	2.15686274509804	1.960784313725491	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.736	0.183	1.116	18.991	4.358	357.6650304886396	74.95968	-2218.623	False	354.0	True
2.15686274509804	2.549019607843138	2.3529411764705888	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.656	0.186	1.042	14.214	3.770	363.83338136609376	66.89659	-2452.717	False	367.0	True
2.549019607843138	2.941176470588236	2.745098039215687	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.450	0.179	0.821	7.011	2.648	339.4818273314676	45.943058	-2174.005	False	339.0	True
2.941176470588236	3.3333333333333344	3.137254901960785	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	2.421	0.225	2.885	174.146	13.196	528.0599574489612	247.276	-3887.531	False	531.0	True
3.3333333333333344	3.725490196078432	3.529411764705883	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	2.190	0.209	2.621	169.698	13.027	457.3561643539827	223.84384	-3183.973	False	458.0	True
3.725490196078432	4.117647058823529	3.9215686274509807	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	2.868	0.232	3.345	243.821	15.615	564.4775740614491	293.2505	-4201.048	False	566.0	True
4.117647058823529	4.509803921568627	4.313725490196078	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	1.599	0.205	2.023	81.901	9.050	451.7627320041883	163.63295	-3040.782	False	448.0	True
4.509803921568627	4.901960784313727	4.7058823529411775	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.550	0.182	0.928	10.243	3.201	356.28517250965956	56.336945	-2313.852	False	356.0	True
4.901960784313727	5.294117647058824	5.098039215686276	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.210	0.171	0.566	1.586	1.259	315.49805616415784	21.543198	-1979.035	True	315.0	True
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
13.921568627451002	14.313725490196079	14.11764705882354	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.377	0.187	0.763	4.407	2.099	385.1566764851361	39.272778	-2526.278	False	384.0	True
14.313725490196079	14.705882352941181	14.50980392156863	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	-0.244	0.168	0.103	2.022	-1.422	313.3736443463295	-25.44988	-1936.261	True	313.0	True
14.705882352941181	15.098039215686308	14.901960784313744	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.307	0.179	0.679	3.133	1.770	366.322547445814	32.01158	-2178.012	True	357.0	True
15.098039215686308	15.490196078431385	15.294117647058847	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	-0.083	0.167	0.265	0.240	-0.490	313.4101919680402	-8.618814	-1914.901	True	311.0	True
15.490196078431385	15.882352941176464	15.686274509803924	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.228	0.172	0.584	1.863	1.365	330.89113194006126	23.80037	-2016.119	True	327.0	True
15.882352941176464	16.274509803921593	16.078431372549026	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.427	0.170	0.780	6.979	2.642	324.69019708074285	44.554733	-1957.576	False	321.0	True
16.274509803921593	16.666666666666693	16.470588235294144	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.760	0.195	1.163	17.553	4.190	428.005206597181	79.45124	-2767.621	False	421.0	True
16.666666666666693	17.05882352941177	16.862745098039234	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.675	0.197	1.082	13.277	3.644	437.13415063137586	70.588745	-2830.574	False	431.0	True
17.05882352941177	17.4509803921569	17.254901960784338	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.256	0.183	0.634	2.057	1.434	381.24750962187824	26.727816	-2356.744	True	374.0	True
17.4509803921569	17.843137254902004	17.647058823529452	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.174	0.182	0.551	0.951	0.975	366.45334911783	18.242867	-2463.364	True	370.0	True
17.843137254902004	18.23529411764708	18.03921568627454	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.593	0.192	0.990	10.705	3.272	408.05504964497743	62.12828	-2803.387	False	410.0	True
18.23529411764708	18.627450980392155	18.431372549019617	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.306	0.173	0.666	3.342	1.828	334.14141424669504	32.074398	-2186.572	True	336.0	True
18.627450980392155	19.01960784313726	18.823529411764707	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	0.010	0.124	0.269	0.005	0.071	171.8115649139593	0.9959075	-877.032	True	172.0	True
19.01960784313726	19.411764705882383	19.21568627450982	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	nan	nan	nan	nan	nan	nan	0.0	0.000	False	0.0	False
19.411764705882383	19.803921568627494	19.60784313725494	70710678.119	10000000.000	500000000.000	4.428e-10	3.043e-01	9.166e+06	nan	nan	nan	nan	nan	nan	0.0	0.000	False	0.0	False

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

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:

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)

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

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),
)

profile = flux_profile_estimator.run(datasets=dataset)

We can directly plot the result:

ax = profile.plot(axis_name="projected-distance", sed_type="flux")

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:

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

And now plot the points together with the likelihood profiles:

ax = profile_high.plot(sed_type="eflux", color="tab:orange")
profile_high.plot_ts_profiles(ax=ax, sed_type="eflux")
ax.set_yscale("linear")