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


[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");


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");

/home/runner/micromamba-root/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

[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 0x7f39c55c5730>
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 17.9 s, sys: 18 ms, total: 18 s
Wall time: 18.1 s

[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")


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


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


### 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")


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_refe_mine_maxref_dnderef_fluxref_efluxnormnorm_errnorm_ultssqrt_tsnprednpred_excessstatis_ulcountssuccess
degdegdegkeVkeVkeV1 / (cm2 s TeV)keV / (cm2 s TeV)keV2 / (cm2 s TeV)
float64float64float64float64[1]float64[1]float64[1]float64[1]float64[1]float64[1]float64[1]float64[1]float64[1]float64[1]float64[1]float64[1,1]float32[1,1]float64[1]bool[1]float64[1,1]bool[1]
-0.196078431372549180.196078431372549180.070710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
0.196078431372549180.58823529411764660.392156862745097970710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
0.58823529411764660.98039215686274410.784313725490195370710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.1720.1240.4342.0591.435162.2304898612695517.448376-818.848True163.0True
0.98039215686274411.37254901960784341.176470588235293770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+061.8940.2072.322120.59110.981450.20181889164763192.60909-3013.446False448.0True
1.37254901960784341.7647058823529421.568627450980392870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+063.1460.2403.639280.64816.753601.6822045018936320.15994-4348.966False599.0True
1.7647058823529422.156862745098041.96078431372549170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7360.1831.11618.9914.358357.665030487863274.95968-2218.623False354.0True
2.156862745098042.5490196078431382.352941176470588870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6560.1861.04214.2143.770363.833381365701366.89659-2452.717False367.0True
2.5490196078431382.94117647058823552.745098039215686770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4500.1790.8217.0112.648339.4818273315226545.943058-2174.005False339.0True
2.94117647058823553.3333333333333343.137254901960784770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+062.4210.2252.885174.14613.196528.0599574498646247.276-3887.531False531.0True
............................................................
15.88235294117646416.27450980392159316.07843137254902670710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4270.1700.7806.9792.642324.690197080841944.554733-1957.576False321.0True
16.27450980392159316.66666666666669316.47058823529414470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7600.1951.16317.5534.190428.0052065967902379.45124-2767.621False421.0True
16.66666666666669317.0588235294117716.86274509803923470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6750.1971.08213.2773.644437.1341506314762570.588745-2830.574False431.0True
17.0588235294117717.450980392156917.25490196078433870710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.2560.1830.6342.0571.434381.2475096215588326.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.055049644950662.12828-2803.387False410.0True
18.2352941176470818.62745098039215518.43137254901961770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3060.1730.6663.3421.828334.14141424656932.074398-2186.572True336.0True
18.62745098039215519.0196078431372618.82352941176470770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.0100.1240.2690.0050.071171.811564914664730.9959075-877.032True172.0True
19.0196078431372619.4117647058823919.21568627450982470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06nannannannannannan0.00.000False0.0False
19.4117647058823919.80392156862749419.60784313725494270710678.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"),
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");

/home/runner/micromamba-root/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:

[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")


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

[ ]: