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


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 0x13f9adf40>
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 6.87 s, sys: 121 ms, total: 6.99 s
Wall time: 6.76 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_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.2304898611595517.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.6822045016994320.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.3561643542095223.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.76273200537855163.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
5.2941176470588245.6862745098039235.490196078431373570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06-0.1470.1590.1840.824-0.908271.28752712430287-15.065349-1643.463True272.0True
5.6862745098039236.078431372549025.88235294117647170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3020.1620.6403.7661.941282.3287388310083731.023876-1692.050True282.0True
6.078431372549026.4705882352941186.27450980392156970710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.0930.1720.4510.3010.549317.73307633550269.594112-2034.334True319.0True
6.4705882352941186.8627450980392166.66666666666666770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06-0.5000.170-0.1477.804-2.794311.3167047326352-51.37142-1931.193True311.0True
............................................................
12.35294117647061912.74509803921569612.54901960784315770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4410.1970.8485.3812.320424.972106268764745.757168-2877.351False425.0True
12.74509803921569613.13725490196077312.94117647058823670710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4920.1940.8937.0182.649416.758847090610551.088707-2705.136False412.0True
13.13725490196077313.52941176470587513.33333333333332570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.9590.2051.38225.8795.087455.2250730519383599.72924-3216.032False458.0True
13.52941176470587513.92156862745100213.7254901960784470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6110.1840.99312.4653.531376.6673329372988563.567978-2403.516False373.0True
13.92156862745100214.3137254901960814.11764705882354170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3770.1870.7634.4072.099385.156676485136139.272778-2526.278False384.0True
14.3137254901960814.70588235294118514.50980392156863370710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06-0.2440.1680.1032.022-1.422313.3736443463295-25.44988-1936.261True313.0True
14.70588235294118515.09803921568631114.90196078431374770710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3070.1790.6793.1331.770366.32254744581432.01158-2178.012True357.0True
15.09803921568631115.49019607843138915.29411764705885170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+06-0.0830.1670.2650.240-0.490313.4101919680869-8.618814-1914.901True311.0True
15.49019607843138915.88235294117646615.68627450980392670710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.2280.1720.5841.8631.365330.8911319400612623.80037-2016.119True327.0True
15.88235294117646616.27450980392159616.0784313725490370710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.4270.1700.7806.9792.642324.6901970807428544.554733-1957.576False321.0True
16.27450980392159616.66666666666669616.47058823529414470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.7600.1951.16317.5534.190428.00520659718179.45124-2767.621False421.0True
16.66666666666669617.05882352941177516.86274509803923470710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.6750.1971.08213.2773.644437.1341506313758670.588745-2830.574False431.0True
17.05882352941177517.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.23529411764708318.03921568627454570710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.5930.1920.99010.7053.272408.0550496449774362.12828-2803.387False410.0True
18.23529411764708318.62745098039215818.4313725490196270710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.3060.1730.6663.3421.828334.1414142466950432.074398-2186.572True336.0True
18.62745098039215819.0196078431372618.8235294117647170710678.11910000000.000500000000.0004.428e-103.043e-019.166e+060.0100.1240.2690.0050.071171.81156491395930.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");

/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/scipy/optimize/_numdiff.py:579: RuntimeWarning: invalid value encountered in true_divide
J_transposed[i] = df / dx


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

[ ]: