Note
Go to the end to download the full example code or to run this example in your browser via Binder
Flux Profile Estimation#
Learn how to estimate flux profiles on a Fermi-LAT dataset.
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 spatial distribution of flux in images and data cubes is the measurement of flux 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
SkyRegion to define the structure of the bins of the flux profile and
run the flux profile extraction using the FluxProfileEstimator
import numpy as np
from astropy import units as u
from astropy.coordinates import SkyCoord
# %matplotlib inline
import matplotlib.pyplot as plt
Setup#
from IPython.display import display
from gammapy.datasets import MapDataset
from gammapy.estimators import FluxPoints, FluxProfileEstimator
from gammapy.maps import RegionGeom
from gammapy.modeling.models import PowerLawSpectralModel
Check setup#
from gammapy.utils.check import check_tutorials_setup
from gammapy.utils.regions import (
make_concentric_annulus_sky_regions,
make_orthogonal_rectangle_sky_regions,
)
check_tutorials_setup()
System:
python_executable : /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/bin/python
python_version : 3.9.18
machine : x86_64
system : Darwin
Gammapy package:
version : 1.1.dev2518+gc7608209d.d20231208
path : /Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy
Other packages:
numpy : 1.26.1
scipy : 1.11.3
astropy : 5.2.2
regions : 0.7
click : 8.1.7
yaml : 6.0.1
IPython : 8.16.1
jupyterlab : not installed
matplotlib : 3.8.0
pandas : not installed
healpy : 1.16.6
iminuit : 2.24.0
sherpa : 4.16.0
naima : 0.10.0
emcee : 3.1.4
corner : 2.2.2
ray : 2.7.1
Gammapy environment variables:
GAMMAPY_DATA : /Users/mregeard/Workspace/dev/code/gammapy/gammapy-data
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")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:83: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
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.
Its starts from lon = 10 deg and 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")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.9/site-packages/regions/shapes/rectangle.py:208: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties.
return Rectangle(xy=xy, width=width, height=height,
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:127: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
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 0x14d976940>
n_jobs : None
n_sigma : 1
n_sigma_ul : 2
norm : Parameter(name='norm', value=1.0, factor=1.0, scale=1.0, unit=Unit(dimensionless), min=nan, max=nan, frozen=False, prior=None, id=0x14d976130)
null_value : 0
parallel_backend : None
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:
profile = flux_profile_estimator.run(datasets=dataset)
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 plot:
ax = profile.plot(sed_type="dnde")
ax.set_yscale("linear")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:176: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
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
plt.figure()
ax = profile.plot(sed_type="eflux")
ax.set_yscale("linear")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:190: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
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"]
fig, ax = plt.subplots()
for quantity in quantities:
profile[quantity].plot(ax=ax, label=quantity.title())
ax.set_ylabel("Counts")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:207: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
Serialisation and I/O#
The profile can be serialised using 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")
ax = profile_new.plot()
ax.set_yscale("linear")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:229: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
The profile can be serialised to a Table object
using:
table = profile.to_table(format="profile", formatted=True)
display(table)
x_min x_max x_ref e_ref e_min ... stat is_ul counts success
deg deg deg keV keV ...
------------------- ------------------ ------------------ ------------ ------------ ... --------- ----- ------ -------
-0.1960784313725492 0.1960784313725492 0.0 70710678.119 10000000.000 ... 0.000 False 0.0 False
0.1960784313725492 0.5882352941176467 0.392156862745098 70710678.119 10000000.000 ... 0.000 False 0.0 False
0.5882352941176467 0.9803921568627443 0.7843137254901955 70710678.119 10000000.000 ... -818.848 True 163.0 True
0.9803921568627443 1.3725490196078436 1.176470588235294 70710678.119 10000000.000 ... -3013.446 False 448.0 True
1.3725490196078436 1.764705882352942 1.5686274509803928 70710678.119 10000000.000 ... -4348.966 False 599.0 True
1.764705882352942 2.15686274509804 1.960784313725491 70710678.119 10000000.000 ... -2218.623 False 354.0 True
2.15686274509804 2.549019607843138 2.3529411764705888 70710678.119 10000000.000 ... -2452.717 False 367.0 True
2.549019607843138 2.941176470588236 2.745098039215687 70710678.119 10000000.000 ... -2174.005 False 339.0 True
2.941176470588236 3.3333333333333344 3.137254901960785 70710678.119 10000000.000 ... -3887.531 False 531.0 True
3.3333333333333344 3.725490196078432 3.529411764705883 70710678.119 10000000.000 ... -3183.973 False 458.0 True
3.725490196078432 4.117647058823529 3.9215686274509807 70710678.119 10000000.000 ... -4201.048 False 566.0 True
4.117647058823529 4.509803921568627 4.313725490196078 70710678.119 10000000.000 ... -3040.782 False 448.0 True
4.509803921568627 4.901960784313727 4.7058823529411775 70710678.119 10000000.000 ... -2313.852 False 356.0 True
4.901960784313727 5.294117647058824 5.098039215686276 70710678.119 10000000.000 ... -1979.035 True 315.0 True
5.294117647058824 5.686274509803923 5.4901960784313735 70710678.119 10000000.000 ... -1643.463 True 272.0 True
5.686274509803923 6.07843137254902 5.882352941176471 70710678.119 10000000.000 ... -1692.050 True 282.0 True
6.07843137254902 6.470588235294118 6.274509803921569 70710678.119 10000000.000 ... -2034.334 True 319.0 True
6.470588235294118 6.862745098039216 6.666666666666667 70710678.119 10000000.000 ... -1931.193 True 311.0 True
6.862745098039216 7.254901960784315 7.058823529411765 70710678.119 10000000.000 ... -1861.277 True 303.0 True
7.254901960784315 7.647058823529413 7.450980392156864 70710678.119 10000000.000 ... -1577.587 True 263.0 True
7.647058823529413 8.039215686274511 7.843137254901962 70710678.119 10000000.000 ... -1754.528 True 288.0 True
8.039215686274511 8.43137254901961 8.235294117647062 70710678.119 10000000.000 ... -2098.681 True 328.0 True
... ... ... ... ... ... ... ... ... ...
10.784313725490232 11.176470588235311 10.980392156862772 70710678.119 10000000.000 ... -3224.239 False 473.0 True
11.176470588235311 11.568627450980387 11.372549019607849 70710678.119 10000000.000 ... -3256.816 False 471.0 True
11.568627450980387 11.960784313725515 11.764705882352951 70710678.119 10000000.000 ... -3129.927 False 453.0 True
11.960784313725515 12.352941176470619 12.156862745098067 70710678.119 10000000.000 ... -2782.182 True 408.0 True
12.352941176470619 12.745098039215696 12.549019607843157 70710678.119 10000000.000 ... -2877.351 False 425.0 True
12.745098039215696 13.137254901960773 12.941176470588236 70710678.119 10000000.000 ... -2705.136 False 412.0 True
13.137254901960773 13.529411764705875 13.333333333333325 70710678.119 10000000.000 ... -3216.032 False 458.0 True
13.529411764705875 13.921568627451002 13.72549019607844 70710678.119 10000000.000 ... -2403.516 False 373.0 True
13.921568627451002 14.313725490196079 14.11764705882354 70710678.119 10000000.000 ... -2526.278 False 384.0 True
14.313725490196079 14.705882352941181 14.50980392156863 70710678.119 10000000.000 ... -1936.261 True 313.0 True
14.705882352941181 15.098039215686308 14.901960784313744 70710678.119 10000000.000 ... -2178.012 True 357.0 True
15.098039215686308 15.490196078431385 15.294117647058847 70710678.119 10000000.000 ... -1914.901 True 311.0 True
15.490196078431385 15.882352941176464 15.686274509803924 70710678.119 10000000.000 ... -2016.119 True 327.0 True
15.882352941176464 16.274509803921593 16.078431372549026 70710678.119 10000000.000 ... -1957.576 False 321.0 True
16.274509803921593 16.666666666666693 16.470588235294144 70710678.119 10000000.000 ... -2767.621 False 421.0 True
16.666666666666693 17.05882352941177 16.862745098039234 70710678.119 10000000.000 ... -2830.574 False 431.0 True
17.05882352941177 17.4509803921569 17.254901960784338 70710678.119 10000000.000 ... -2356.744 True 374.0 True
17.4509803921569 17.843137254902004 17.647058823529452 70710678.119 10000000.000 ... -2463.364 True 370.0 True
17.843137254902004 18.23529411764708 18.03921568627454 70710678.119 10000000.000 ... -2803.387 False 410.0 True
18.23529411764708 18.627450980392155 18.431372549019617 70710678.119 10000000.000 ... -2186.572 True 336.0 True
18.627450980392155 19.01960784313726 18.823529411764707 70710678.119 10000000.000 ... -877.032 True 172.0 True
19.01960784313726 19.411764705882383 19.21568627450982 70710678.119 10000000.000 ... 0.000 False 0.0 False
19.411764705882383 19.803921568627494 19.60784313725494 70710678.119 10000000.000 ... 0.000 False 0.0 False
Length = 51 rows
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")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/.tox/build_docs/lib/python3.9/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)
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:262: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
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"],
)
The configuration of the fit statistic profile is done throught the norm parameter:
flux_profile_estimator.norm.scan_values = np.linspace(-1, 5, 11)
Now we can run the estimator,
profile = flux_profile_estimator.run(datasets=dataset)
and plot the result:
profile.plot(axis_name="projected-distance", sed_type="flux")
plt.show()
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:290: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
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)})
plt.show()
And now plot the points together with the likelihood profiles:
fig, ax = plt.subplots()
profile_high.plot(ax=ax, sed_type="eflux", color="tab:orange")
profile_high.plot_ts_profiles(ax=ax, sed_type="eflux")
ax.set_yscale("linear")
plt.show()
# sphinx_gallery_thumbnail_number = 2
/Users/mregeard/Workspace/dev/code/gammapy/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:310: UserWarning: FigureCanvasAgg is non-interactive, and thus cannot be shown
plt.show()
Total running time of the script: ( 0 minutes 10.087 seconds)
