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");
/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,
../../../_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 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")
../../../_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_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"),
    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");
/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)
../../../_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
[ ]: