.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-3d/flux_profiles.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_analysis-3d_flux_profiles.py: 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 `~gammapy.datasets.MapDataset` of Fermi-LAT data, use `~regions.SkyRegion` to define the structure of the bins of the flux profile and run the flux profile extraction using the `~gammapy.estimators.FluxProfileEstimator` .. GENERATED FROM PYTHON SOURCE LINES 37-45 .. code-block:: Python import numpy as np from astropy import units as u from astropy.coordinates import SkyCoord # %matplotlib inline import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 46-49 Setup ----- .. GENERATED FROM PYTHON SOURCE LINES 49-60 .. code-block:: Python 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 from gammapy.utils.regions import ( make_concentric_annulus_sky_regions, make_orthogonal_rectangle_sky_regions, ) .. GENERATED FROM PYTHON SOURCE LINES 61-64 Read and Introduce Data ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 64-70 .. code-block:: Python dataset = MapDataset.read( "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz", name="fermi-dataset" ) .. GENERATED FROM PYTHON SOURCE LINES 71-73 This is what the counts image we will work with looks like: .. GENERATED FROM PYTHON SOURCE LINES 73-78 .. code-block:: Python counts_image = dataset.counts.sum_over_axes() counts_image.smooth("0.1 deg").plot(stretch="sqrt") plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_001.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 79-82 There are 400x200 pixels in the dataset and 11 energy bins between 10 GeV and 2 TeV: .. GENERATED FROM PYTHON SOURCE LINES 82-86 .. code-block:: Python print(dataset.counts) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (np.int64(400), np.int64(200), 11) ndim : 3 unit : dtype : >i4 .. GENERATED FROM PYTHON SOURCE LINES 87-101 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 `~gammapy.utils.regions.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: .. GENERATED FROM PYTHON SOURCE LINES 101-111 .. code-block:: Python regions = make_orthogonal_rectangle_sky_regions( start_pos=SkyCoord("9d", "0d", frame="galactic"), end_pos=SkyCoord("351d", "0d", frame="galactic"), wcs=counts_image.geom.wcs, height="3 deg", nbin=49, ) .. GENERATED FROM PYTHON SOURCE LINES 112-115 We can use the `~gammapy.maps.RegionGeom` object to illustrate the regions on top of the counts image: .. GENERATED FROM PYTHON SOURCE LINES 115-122 .. code-block:: Python geom = RegionGeom.create(region=regions) ax = counts_image.smooth("0.1 deg").plot(stretch="sqrt") geom.plot_region(ax=ax, color="w") plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_002.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.11/site-packages/regions/shapes/rectangle.py:205: UserWarning: Setting the 'color' property will override the edgecolor or facecolor properties. return Rectangle(xy=xy, width=width, height=height, .. GENERATED FROM PYTHON SOURCE LINES 123-127 Next we create the `~gammapy.estimators.FluxProfileEstimator`. For the estimation of the flux profile we assume a spectral model with a power-law shape and an index of 2.3 .. GENERATED FROM PYTHON SOURCE LINES 127-136 .. code-block:: Python flux_profile_estimator = FluxProfileEstimator( regions=regions, spectral_model=PowerLawSpectralModel(index=2.3), energy_edges=[10, 2000] * u.GeV, selection_optional=["ul"], ) .. GENERATED FROM PYTHON SOURCE LINES 137-139 We can see the full configuration by printing the estimator object: .. GENERATED FROM PYTHON SOURCE LINES 139-143 .. code-block:: Python print(flux_profile_estimator) .. rst-class:: sphx-glr-script-out .. code-block:: none FluxProfileEstimator -------------------- energy_edges : [ 10. 2000.] GeV fit : n_jobs : None n_sigma : 1 n_sigma_sensitivity : 5 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=0x7efc0a154250) null_value : 0 parallel_backend : None reoptimize : False selection_optional : ['ul'] source : 0 spectral_model : PowerLawSpectralModel sum_over_energy_groups : False .. GENERATED FROM PYTHON SOURCE LINES 144-149 Run Estimation ~~~~~~~~~~~~~~ Now we can run the profile estimation and explore the results: .. GENERATED FROM PYTHON SOURCE LINES 149-155 .. code-block:: Python profile = flux_profile_estimator.run(datasets=dataset) print(profile) .. rst-class:: sphx-glr-script-out .. code-block:: none FluxPoints ---------- geom : RegionGeom axes : ['lon', 'lat', 'energy', 'projected-distance'] shape : (1, 1, 1, 49) quantities : ['norm', 'norm_err', 'norm_ul', 'ts', 'npred', 'npred_excess', 'stat', 'stat_null', 'counts', 'success'] ref. model : pl n_sigma : 1 n_sigma_ul : 2 sqrt_ts_threshold_ul : 2 sed type init : likelihood .. GENERATED FROM PYTHON SOURCE LINES 156-165 We can see the flux profile is represented by a `~gammapy.estimators.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 `~gammapy.estimators.FluxPoints.plot`: .. GENERATED FROM PYTHON SOURCE LINES 165-170 .. code-block:: Python ax = profile.plot(sed_type="dnde") ax.set_yscale("linear") plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_003.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 171-175 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: .. GENERATED FROM PYTHON SOURCE LINES 175-184 .. code-block:: Python profile.sqrt_ts_threshold_ul = 2 plt.figure() ax = profile.plot(sed_type="eflux") ax.set_yscale("linear") plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_004.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 185-189 We can also plot any other quantity of interest, that is defined on the `~gammapy.estimators.FluxPoints` result object. E.g. the predicted total counts, background counts and excess counts: .. GENERATED FROM PYTHON SOURCE LINES 189-202 .. code-block:: Python 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") ax.legend() plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_005.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 203-209 Serialisation and I/O ~~~~~~~~~~~~~~~~~~~~~ The profile can be serialised using `~gammapy.estimators.FluxPoints.write`, given a specific format: .. GENERATED FROM PYTHON SOURCE LINES 209-224 .. code-block:: Python 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() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_006.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 225-228 The profile can be serialised to a `~astropy.table.Table` object using: .. GENERATED FROM PYTHON SOURCE LINES 228-233 .. code-block:: Python table = profile.to_table(format="profile", formatted=True) display(table) .. rst-class:: sphx-glr-script-out .. code-block:: none x_min x_max x_ref ... counts success deg deg deg ... -------------------- ------------------- ------------------- ... ------ ------- -0.18367346938775508 0.18367346938775508 0.0 ... 328.0 True 0.18367346938775508 0.5510204081632653 0.36734693877551017 ... 607.0 True 0.5510204081632653 0.9183673469387759 0.7346938775510206 ... 405.0 True 0.9183673469387759 1.285714285714286 1.102040816326531 ... 322.0 True 1.285714285714286 1.6530612244897958 1.469387755102041 ... 344.0 True 1.6530612244897958 2.0204081632653064 1.836734693877551 ... 311.0 True 2.0204081632653064 2.387755102040817 2.204081632653062 ... 484.0 True 2.387755102040817 2.755102040816327 2.571428571428572 ... 537.0 True 2.755102040816327 3.1224489795918373 2.9387755102040822 ... 487.0 True ... ... ... ... ... ... 14.142857142857132 14.510204081632645 14.326530612244888 ... 267.0 True 14.510204081632645 14.877551020408161 14.693877551020403 ... 291.0 True 14.877551020408161 15.244897959183646 15.061224489795904 ... 357.0 True 15.244897959183646 15.612244897959162 15.428571428571404 ... 374.0 True 15.612244897959162 15.979591836734677 15.79591836734692 ... 386.0 True 15.979591836734677 16.34693877551019 16.163265306122433 ... 368.0 True 16.34693877551019 16.714285714285705 16.530612244897945 ... 321.0 True 16.714285714285705 17.081632653061217 16.89795918367346 ... 361.0 True 17.081632653061217 17.44897959183673 17.265306122448973 ... 385.0 True 17.44897959183673 17.81632653061222 17.632653061224474 ... 319.0 True Length = 49 rows .. GENERATED FROM PYTHON SOURCE LINES 234-237 No we can also estimate a radial profile starting from the Galactic center: .. GENERATED FROM PYTHON SOURCE LINES 237-245 .. code-block:: Python regions = make_concentric_annulus_sky_regions( center=SkyCoord("0d", "0d", frame="galactic"), radius_max="1.5 deg", nbin=11, ) .. GENERATED FROM PYTHON SOURCE LINES 246-248 Again we first illustrate the regions: .. GENERATED FROM PYTHON SOURCE LINES 248-257 .. code-block:: Python 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() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_007.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_007.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.11/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) .. GENERATED FROM PYTHON SOURCE LINES 258-260 This time we define two energy bins and include the fit statistic profile in the computation: .. GENERATED FROM PYTHON SOURCE LINES 260-267 .. code-block:: Python flux_profile_estimator = FluxProfileEstimator( regions=regions, spectral_model=PowerLawSpectralModel(index=2.3), energy_edges=[10, 100, 2000] * u.GeV, selection_optional=["ul", "scan"], ) .. GENERATED FROM PYTHON SOURCE LINES 268-269 The configuration of the fit statistic profile is done throught the norm parameter: .. GENERATED FROM PYTHON SOURCE LINES 269-271 .. code-block:: Python flux_profile_estimator.norm.scan_values = np.linspace(-1, 5, 11) .. GENERATED FROM PYTHON SOURCE LINES 272-273 Now we can run the estimator, .. GENERATED FROM PYTHON SOURCE LINES 273-277 .. code-block:: Python profile = flux_profile_estimator.run(datasets=dataset) .. GENERATED FROM PYTHON SOURCE LINES 278-280 and plot the result: .. GENERATED FROM PYTHON SOURCE LINES 280-285 .. code-block:: Python profile.plot(axis_name="projected-distance", sed_type="flux") plt.show() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_008.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 286-289 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: .. GENERATED FROM PYTHON SOURCE LINES 289-294 .. code-block:: Python profile_high = profile.slice_by_idx({"energy": slice(1, 2)}) plt.show() .. GENERATED FROM PYTHON SOURCE LINES 295-297 And now plot the points together with the likelihood profiles: .. GENERATED FROM PYTHON SOURCE LINES 297-304 .. code-block:: Python 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() .. image-sg:: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_009.png :alt: flux profiles :srcset: /tutorials/analysis-3d/images/sphx_glr_flux_profiles_009.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 22.077 seconds) .. _sphx_glr_download_tutorials_analysis-3d_flux_profiles.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v2.0?urlpath=lab/tree/notebooks/2.0/tutorials/analysis-3d/flux_profiles.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: flux_profiles.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: flux_profiles.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: flux_profiles.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_