.. 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-55 .. 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 .. GENERATED FROM PYTHON SOURCE LINES 56-58 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 58-67 .. code-block:: Python 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() .. rst-class:: sphx-glr-script-out .. code-block:: none System: python_executable : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/bin/python python_version : 3.9.20 machine : x86_64 system : Linux Gammapy package: version : 1.3 path : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.13.1 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.2 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.2 pandas : not installed healpy : 1.17.3 iminuit : 2.30.1 sherpa : 4.16.1 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.3 ray : 2.39.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3 .. GENERATED FROM PYTHON SOURCE LINES 68-71 Read and Introduce Data ----------------------- .. GENERATED FROM PYTHON SOURCE LINES 71-77 .. code-block:: Python dataset = MapDataset.read( "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc.fits.gz", name="fermi-dataset" ) .. GENERATED FROM PYTHON SOURCE LINES 78-80 This is what the counts image we will work with looks like: .. GENERATED FROM PYTHON SOURCE LINES 80-85 .. 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 86-89 There are 400x200 pixels in the dataset and 11 energy bins between 10 GeV and 2 TeV: .. GENERATED FROM PYTHON SOURCE LINES 89-93 .. code-block:: Python print(dataset.counts) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (400, 200, 11) ndim : 3 unit : dtype : >i4 .. GENERATED FROM PYTHON SOURCE LINES 94-108 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 108-118 .. code-block:: Python 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, ) .. GENERATED FROM PYTHON SOURCE LINES 119-122 We can use the `~gammapy.maps.RegionGeom` object to illustrate the regions on top of the counts image: .. GENERATED FROM PYTHON SOURCE LINES 122-129 .. 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.9/site-packages/regions/shapes/rectangle.py:207: 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 130-134 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 134-143 .. code-block:: Python flux_profile_estimator = FluxProfileEstimator( regions=regions, spectrum=PowerLawSpectralModel(index=2.3), energy_edges=[10, 2000] * u.GeV, selection_optional=["ul"], ) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:135: GammapyDeprecationWarning: "spectrum" was deprecated in version v1.3 and will be removed in a future version. Use argument "spectral_model" instead. flux_profile_estimator = FluxProfileEstimator( .. GENERATED FROM PYTHON SOURCE LINES 144-146 We can see the full configuration by printing the estimator object: .. GENERATED FROM PYTHON SOURCE LINES 146-150 .. 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_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=0x7f46b56a4940) 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 151-156 Run Estimation ~~~~~~~~~~~~~~ Now we can run the profile estimation and explore the results: .. GENERATED FROM PYTHON SOURCE LINES 158-163 .. 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, 51) 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 164-173 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 173-178 .. 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 179-183 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 183-192 .. 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 193-197 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 197-210 .. 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 211-217 Serialisation and I/O ~~~~~~~~~~~~~~~~~~~~~ The profile can be serialised using `~gammapy.estimators.FluxPoints.write`, given a specific format: .. GENERATED FROM PYTHON SOURCE LINES 217-232 .. 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 233-236 The profile can be serialised to a `~astropy.table.Table` object using: .. GENERATED FROM PYTHON SOURCE LINES 236-241 .. 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.1960784313725492 0.1960784313725492 0.0 ... 0.0 False 0.1960784313725492 0.5882352941176467 0.392156862745098 ... 0.0 False 0.5882352941176467 0.9803921568627443 0.7843137254901955 ... 163.0 True 0.9803921568627443 1.3725490196078436 1.176470588235294 ... 448.0 True 1.3725490196078436 1.764705882352942 1.5686274509803928 ... 599.0 True 1.764705882352942 2.15686274509804 1.960784313725491 ... 354.0 True 2.15686274509804 2.549019607843138 2.3529411764705888 ... 367.0 True 2.549019607843138 2.941176470588236 2.745098039215687 ... 339.0 True 2.941176470588236 3.3333333333333344 3.137254901960785 ... 531.0 True 3.3333333333333344 3.725490196078432 3.529411764705883 ... 458.0 True ... ... ... ... ... ... 15.882352941176466 16.274509803921596 16.07843137254903 ... 321.0 True 16.274509803921596 16.666666666666696 16.470588235294144 ... 421.0 True 16.666666666666696 17.058823529411775 16.862745098039234 ... 431.0 True 17.058823529411775 17.4509803921569 17.254901960784338 ... 374.0 True 17.4509803921569 17.843137254902004 17.647058823529452 ... 370.0 True 17.843137254902004 18.235294117647083 18.039215686274545 ... 410.0 True 18.235294117647083 18.627450980392158 18.43137254901962 ... 336.0 True 18.627450980392158 19.01960784313726 18.82352941176471 ... 172.0 True 19.01960784313726 19.41176470588239 19.215686274509824 ... 0.0 False 19.41176470588239 19.803921568627494 19.607843137254942 ... 0.0 False Length = 51 rows .. GENERATED FROM PYTHON SOURCE LINES 242-245 No we can also estimate a radial profile starting from the Galactic center: .. GENERATED FROM PYTHON SOURCE LINES 245-253 .. 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 254-256 Again we first illustrate the regions: .. GENERATED FROM PYTHON SOURCE LINES 256-265 .. 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.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) .. GENERATED FROM PYTHON SOURCE LINES 266-268 This time we define two energy bins and include the fit statistic profile in the computation: .. GENERATED FROM PYTHON SOURCE LINES 268-275 .. code-block:: Python flux_profile_estimator = FluxProfileEstimator( regions=regions, spectrum=PowerLawSpectralModel(index=2.3), energy_edges=[10, 100, 2000] * u.GeV, selection_optional=["ul", "scan"], ) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/examples/tutorials/analysis-3d/flux_profiles.py:269: GammapyDeprecationWarning: "spectrum" was deprecated in version v1.3 and will be removed in a future version. Use argument "spectral_model" instead. flux_profile_estimator = FluxProfileEstimator( .. GENERATED FROM PYTHON SOURCE LINES 276-277 The configuration of the fit statistic profile is done throught the norm parameter: .. GENERATED FROM PYTHON SOURCE LINES 277-279 .. code-block:: Python flux_profile_estimator.norm.scan_values = np.linspace(-1, 5, 11) .. GENERATED FROM PYTHON SOURCE LINES 280-281 Now we can run the estimator, .. GENERATED FROM PYTHON SOURCE LINES 281-285 .. code-block:: Python profile = flux_profile_estimator.run(datasets=dataset) .. GENERATED FROM PYTHON SOURCE LINES 286-288 and plot the result: .. GENERATED FROM PYTHON SOURCE LINES 288-293 .. 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 294-297 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 297-302 .. code-block:: Python profile_high = profile.slice_by_idx({"energy": slice(1, 2)}) plt.show() .. GENERATED FROM PYTHON SOURCE LINES 303-305 And now plot the points together with the likelihood profiles: .. GENERATED FROM PYTHON SOURCE LINES 305-313 .. 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() # sphinx_gallery_thumbnail_number = 2 .. 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.207 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/v1.3?urlpath=lab/tree/notebooks/1.3/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 `_