.. 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.19 machine : x86_64 system : Linux Gammapy package: version : 1.3.dev222+g8af9dbfc9 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.0 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.1 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.8.4 pandas : not installed healpy : 1.16.6 iminuit : 2.25.2 sherpa : 4.16.0 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.2 ray : 2.12.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev .. 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"], ) .. 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=0x7f4142719be0) null_value : 0 parallel_backend : None reoptimize : False selection_optional : ['ul'] source : 0 spectrum : 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-209 .. 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") 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 210-216 Serialisation and I/O ~~~~~~~~~~~~~~~~~~~~~ The profile can be serialised using `~gammapy.estimators.FluxPoints.write`, given a specific format: .. GENERATED FROM PYTHON SOURCE LINES 216-231 .. 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 232-235 The profile can be serialised to a `~astropy.table.Table` object using: .. GENERATED FROM PYTHON SOURCE LINES 235-240 .. 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 241-244 No we can also estimate a radial profile starting from the Galactic center: .. GENERATED FROM PYTHON SOURCE LINES 244-252 .. 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 253-255 Again we first illustrate the regions: .. GENERATED FROM PYTHON SOURCE LINES 255-264 .. 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 265-267 This time we define two energy bins and include the fit statistic profile in the computation: .. GENERATED FROM PYTHON SOURCE LINES 267-274 .. 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"], ) .. GENERATED FROM PYTHON SOURCE LINES 275-276 The configuration of the fit statistic profile is done throught the norm parameter: .. GENERATED FROM PYTHON SOURCE LINES 276-278 .. code-block:: Python flux_profile_estimator.norm.scan_values = np.linspace(-1, 5, 11) .. GENERATED FROM PYTHON SOURCE LINES 279-280 Now we can run the estimator, .. GENERATED FROM PYTHON SOURCE LINES 280-284 .. code-block:: Python profile = flux_profile_estimator.run(datasets=dataset) .. GENERATED FROM PYTHON SOURCE LINES 285-287 and plot the result: .. GENERATED FROM PYTHON SOURCE LINES 287-292 .. 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 293-296 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 296-301 .. code-block:: Python profile_high = profile.slice_by_idx({"energy": slice(1, 2)}) plt.show() .. GENERATED FROM PYTHON SOURCE LINES 302-304 And now plot the points together with the likelihood profiles: .. GENERATED FROM PYTHON SOURCE LINES 304-312 .. 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 20.114 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/main?urlpath=lab/tree/notebooks/dev/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 ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_