.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/details/datasets.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_details_datasets.py: Datasets - Reduced data, IRFs, models ===================================== Learn how to work with datasets Introduction ------------ `gammapy.datasets` are a crucial part of the gammapy API. `~gammapy.datasets.Dataset` objects constitute ``DL4`` data - binned counts, IRFs, models and the associated likelihoods. `~gammapy.datasets.Datasets` from the end product of the data reduction stage, see :doc:`/tutorials/details/makers` tutorial and are passed on to the `~gammapy.modeling.Fit` or estimator classes for modelling and fitting purposes. To find the different types of `~gammapy.datasets.Dataset` objects that are supported see :ref:`datasets-types`: Setup ----- .. GENERATED FROM PYTHON SOURCE LINES 23-44 .. code-block:: Python import astropy.units as u from astropy.coordinates import SkyCoord from regions import CircleSkyRegion import matplotlib.pyplot as plt from IPython.display import display # %matplotlib inline from gammapy.data import GTI from gammapy.datasets import ( Datasets, FluxPointsDataset, MapDataset, SpectrumDatasetOnOff, ) from gammapy.estimators import FluxPoints from gammapy.maps import MapAxis, WcsGeom from gammapy.modeling.models import FoVBackgroundModel, PowerLawSpectralModel, SkyModel from gammapy.utils.scripts import make_path .. GENERATED FROM PYTHON SOURCE LINES 45-62 `~gammapy.datasets.MapDataset` ------------------------------ The counts, exposure, background, masks, and IRF maps are bundled together in a data structure named `~gammapy.datasets.MapDataset`. While the ``counts``, and ``background`` maps are binned in reconstructed energy and must have the same geometry, the IRF maps can have a different spatial (coarsely binned and larger) geometry and spectral range (binned in true energies). It is usually recommended that the true energy bin should be larger and more finely sampled and the reco energy bin. Creating an empty dataset ~~~~~~~~~~~~~~~~~~~~~~~~~ An empty `~gammapy.datasets.MapDataset` can be directly instantiated from any `~gammapy.maps.WcsGeom` object: .. GENERATED FROM PYTHON SOURCE LINES 62-76 .. code-block:: Python energy_axis = MapAxis.from_energy_bounds(1, 10, nbin=11, name="energy", unit="TeV") geom = WcsGeom.create( skydir=(83.63, 22.01), axes=[energy_axis], width=5 * u.deg, binsz=0.05 * u.deg, frame="icrs", ) dataset_empty = MapDataset.create(geom=geom, name="my-dataset") .. GENERATED FROM PYTHON SOURCE LINES 77-81 It is good practice to define a name for the dataset, such that you can identify it later by name. However, if you define a name it **must** be unique. Now we can already print the dataset: .. GENERATED FROM PYTHON SOURCE LINES 81-85 .. code-block:: Python print(dataset_empty) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : my-dataset Total counts : 0 Total background counts : 0.00 Total excess counts : 0.00 Predicted counts : 0.00 Predicted background counts : 0.00 Predicted excess counts : nan Exposure min : 0.00e+00 m2 s Exposure max : 0.00e+00 m2 s Number of total bins : 110000 Number of fit bins : 0 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 86-92 The printout shows the key summary information of the dataset, such as total counts, fit statistics, model information etc. `~gammapy.datasets.MapDataset.from_geom` has additional keywords, that can be used to define the binning of the IRF related maps: .. GENERATED FROM PYTHON SOURCE LINES 92-113 .. code-block:: Python # choose a different true energy binning for the exposure, psf and edisp energy_axis_true = MapAxis.from_energy_bounds( 0.1, 100, nbin=11, name="energy_true", unit="TeV", per_decade=True ) # choose a different rad axis binning for the psf rad_axis = MapAxis.from_bounds(0, 5, nbin=50, unit="deg", name="rad") gti = GTI.create(0 * u.s, 1000 * u.s) dataset_empty = MapDataset.create( geom=geom, energy_axis_true=energy_axis_true, rad_axis=rad_axis, binsz_irf=0.1, gti=gti, name="dataset-empty", ) .. GENERATED FROM PYTHON SOURCE LINES 114-116 To see the geometry of each map, we can use: .. GENERATED FROM PYTHON SOURCE LINES 116-120 .. code-block:: Python print(dataset_empty.geoms) .. rst-class:: sphx-glr-script-out .. code-block:: none {'geom': , 'geom_exposure': , 'geom_psf': , 'geom_edisp': } .. GENERATED FROM PYTHON SOURCE LINES 121-124 Another way to create a `~gammapy.datasets.MapDataset` is to just read an existing one from a FITS file: .. GENERATED FROM PYTHON SOURCE LINES 124-132 .. code-block:: Python dataset_cta = MapDataset.read( "$GAMMAPY_DATA/cta-1dc-gc/cta-1dc-gc.fits.gz", name="dataset-cta" ) print(dataset_cta) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : dataset-cta Total counts : 104317 Total background counts : 91507.70 Total excess counts : 12809.30 Predicted counts : 91507.69 Predicted background counts : 91507.70 Predicted excess counts : nan Exposure min : 6.28e+07 m2 s Exposure max : 1.90e+10 m2 s Number of total bins : 768000 Number of fit bins : 691680 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 133-136 Accessing contents of a dataset ------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 139-142 To further explore the contents of a `~gammapy.datasets.Dataset`, you can use e.g. `~gammapy.datasets.MapDataset.info_dict()` .. GENERATED FROM PYTHON SOURCE LINES 144-146 For a quick info, use .. GENERATED FROM PYTHON SOURCE LINES 146-149 .. code-block:: Python print(dataset_cta.info_dict()) .. rst-class:: sphx-glr-script-out .. code-block:: none {'name': 'dataset-cta', 'counts': 104317, 'excess': 12809.3046875, 'sqrt_ts': np.float64(41.41009347393684), 'background': 91507.6953125, 'npred': 91507.68628538586, 'npred_background': 91507.6953125, 'npred_signal': nan, 'exposure_min': , 'exposure_max': , 'livetime': , 'ontime': , 'counts_rate': , 'background_rate': , 'excess_rate': , 'n_bins': 768000, 'n_fit_bins': 691680, 'stat_type': 'cash', 'stat_sum': nan} .. GENERATED FROM PYTHON SOURCE LINES 150-152 For a quick view, use .. GENERATED FROM PYTHON SOURCE LINES 152-156 .. code-block:: Python dataset_cta.peek() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_001.png :alt: Exposure (at FoV center), Containment radius (at FoV center), Energy Dispersion (at FoV center), Exposure map at 1.17 TeV, Counts map, Model npred :srcset: /tutorials/details/images/sphx_glr_datasets_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 157-159 And access individual maps like: .. GENERATED FROM PYTHON SOURCE LINES 159-165 .. code-block:: Python plt.figure() counts_image = dataset_cta.counts.sum_over_axes() counts_image.smooth("0.1 deg").plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_002.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 166-169 Of course you can also access IRF related maps, e.g. the psf as `~gammapy.irf.PSFMap`: .. GENERATED FROM PYTHON SOURCE LINES 169-173 .. code-block:: Python print(dataset_cta.psf) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'rad', 'energy_true'] shape : (np.int64(16), np.int64(12), 66, 14) ndim : 4 unit : 1 / sr dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 174-176 And use any method on the `~gammapy.irf.PSFMap` object: .. GENERATED FROM PYTHON SOURCE LINES 176-179 .. code-block:: Python radius = dataset_cta.psf.containment_radius(energy_true=1 * u.TeV, fraction=0.95) print(radius) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.08675] deg .. GENERATED FROM PYTHON SOURCE LINES 180-186 .. code-block:: Python edisp_kernel = dataset_cta.edisp.get_edisp_kernel() edisp_kernel.plot_matrix() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_003.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 187-191 The `~gammapy.datasets.MapDataset` typically also contains the information on the residual hadronic background, stored in `~gammapy.datasets.MapDataset.background` as a map: .. GENERATED FROM PYTHON SOURCE LINES 191-195 .. code-block:: Python print(dataset_cta.background) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (np.int64(320), np.int64(240), 10) ndim : 3 unit : dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 196-199 As a next step we define a minimal model on the dataset using the `~gammapy.datasets.MapDataset.models` setter: .. GENERATED FROM PYTHON SOURCE LINES 199-207 .. code-block:: Python model = SkyModel.create("pl", "point", name="gc") model.spatial_model.position = SkyCoord("0d", "0d", frame="galactic") model_bkg = FoVBackgroundModel(dataset_name="dataset-cta") dataset_cta.models = [model, model_bkg] .. GENERATED FROM PYTHON SOURCE LINES 208-211 Assigning models to datasets is covered in more detail in :doc:`/tutorials/details/model_management`. Printing the dataset will now show the model components: .. GENERATED FROM PYTHON SOURCE LINES 211-215 .. code-block:: Python print(dataset_cta) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : dataset-cta Total counts : 104317 Total background counts : 91507.70 Total excess counts : 12809.31 Predicted counts : 91719.65 Predicted background counts : 91507.69 Predicted excess counts : 211.96 Exposure min : 6.28e+07 m2 s Exposure max : 1.90e+10 m2 s Number of total bins : 768000 Number of fit bins : 691680 Fit statistic type : cash Fit statistic value (-2 log(L)) : 424650.15 Number of models : 2 Number of parameters : 8 Number of free parameters : 5 Component 0: SkyModel Name : gc Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: index : 2.000 +/- 0.00 amplitude : 1.00e-12 +/- 0.0e+00 1 / (TeV s cm2) reference (frozen): 1.000 TeV lon_0 : 4.650 +/- 0.00 rad lat_0 : -0.505 +/- 0.00 rad Component 1: FoVBackgroundModel Name : dataset-cta-bkg Datasets names : ['dataset-cta'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.000 +/- 0.00 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 216-219 Now we can use the `~gammapy.datasets.MapDataset.npred` method to get a map of the total predicted counts of the model: .. GENERATED FROM PYTHON SOURCE LINES 219-225 .. code-block:: Python npred = dataset_cta.npred() npred.sum_over_axes().plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_004.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 226-229 To get the predicted counts from an individual model component we can use: .. GENERATED FROM PYTHON SOURCE LINES 229-235 .. code-block:: Python npred_source = dataset_cta.npred_signal(model_names=["gc"]) npred_source.sum_over_axes().plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_005.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_005.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 236-241 `~gammapy.datasets.MapDataset.background` contains the background map computed from the IRF. Internally it will be combined with a `~gammapy.modeling.models.FoVBackgroundModel`, to allow for adjusting the background model during a fit. To get the model corrected background, one can use `~gammapy.datasets.MapDataset.npred_background`. .. GENERATED FROM PYTHON SOURCE LINES 241-247 .. code-block:: Python npred_background = dataset_cta.npred_background() npred_background.sum_over_axes().plot() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_006.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 248-269 Using masks ~~~~~~~~~~~ There are two masks that can be set on a `~gammapy.datasets.MapDataset`, `~gammapy.datasets.MapDataset.mask_safe` and `~gammapy.datasets.MapDataset.mask_fit`. - The `~gammapy.datasets.MapDataset.mask_safe` is computed during the data reduction process according to the specified selection cuts, and should not be changed by the user. - During modelling and fitting, the user might want to additionally ignore some parts of a reduced dataset, e.g. to restrict the fit to a specific energy range or to ignore parts of the region of interest. This should be done by applying the `~gammapy.datasets.MapDataset.mask_fit`. To see details of applying masks, please refer to :ref:`masks-for-fitting`. Both the `~gammapy.datasets.MapDataset.mask_fit` and `~gammapy.datasets.MapDataset.mask_safe` must have the same `~gammapy.maps.Map.geom` as the `~gammapy.datasets.MapDataset.counts` and `~gammapy.datasets.MapDataset.background` maps. For example to see the safe data range: .. GENERATED FROM PYTHON SOURCE LINES 269-274 .. code-block:: Python dataset_cta.mask_safe.plot_grid(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_007.png :alt: Energy 100 GeV - 158 GeV, Energy 158 GeV - 251 GeV, Energy 251 GeV - 398 GeV, Energy 398 GeV - 631 GeV, Energy 631 GeV - 1.00 TeV, Energy 1.00 TeV - 1.58 TeV, Energy 1.58 TeV - 2.51 TeV, Energy 2.51 TeV - 3.98 TeV, Energy 3.98 TeV - 6.31 TeV, Energy 6.31 TeV - 10.0 TeV :srcset: /tutorials/details/images/sphx_glr_datasets_007.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 275-278 In addition it is possible to define a custom `~gammapy.datasets.MapDataset.mask_fit`. Here we apply a mask fit in energy and space: .. GENERATED FROM PYTHON SOURCE LINES 278-291 .. code-block:: Python region = CircleSkyRegion(SkyCoord("0d", "0d", frame="galactic"), 1.5 * u.deg) geom = dataset_cta.counts.geom mask_space = geom.region_mask([region]) mask_energy = geom.energy_mask(0.3 * u.TeV, 8 * u.TeV) dataset_cta.mask_fit = mask_space & mask_energy dataset_cta.mask_fit.plot_grid(vmin=0, vmax=1, add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_008.png :alt: Energy 100 GeV - 158 GeV, Energy 158 GeV - 251 GeV, Energy 251 GeV - 398 GeV, Energy 398 GeV - 631 GeV, Energy 631 GeV - 1.00 TeV, Energy 1.00 TeV - 1.58 TeV, Energy 1.58 TeV - 2.51 TeV, Energy 2.51 TeV - 3.98 TeV, Energy 3.98 TeV - 6.31 TeV, Energy 6.31 TeV - 10.0 TeV :srcset: /tutorials/details/images/sphx_glr_datasets_008.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 292-301 To access the energy range defined by the mask you can use: - `~gammapy.datasets.MapDataset.energy_range_safe` : energy range defined by the `~gammapy.datasets.MapDataset.mask_safe` - `~gammapy.datasets.MapDataset.energy_range_fit` : energy range defined by the `~gammapy.datasets.MapDataset.mask_fit` - `~gammapy.datasets.MapDataset.energy_range` : the final energy range used in likelihood computation These methods return two maps, with the `min` and `max` energy values at each spatial pixel .. GENERATED FROM PYTHON SOURCE LINES 301-309 .. code-block:: Python e_min, e_max = dataset_cta.energy_range # To see the low energy threshold at each point e_min.plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_009.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_009.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 310-311 To see the high energy threshold at each point .. GENERATED FROM PYTHON SOURCE LINES 311-316 .. code-block:: Python e_max.plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_010.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_010.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 317-321 Just as for `~gammapy.maps.Map` objects it is possible to cutout a whole `~gammapy.datasets.MapDataset`, which will perform the cutout for all maps in parallel. Optionally one can provide a new name to the resulting dataset: .. GENERATED FROM PYTHON SOURCE LINES 321-332 .. code-block:: Python cutout = dataset_cta.cutout( position=SkyCoord("0d", "0d", frame="galactic"), width=2 * u.deg, name="cta-cutout", ) cutout.counts.sum_over_axes().plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_011.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_011.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 333-335 It is also possible to slice a `~gammapy.datasets.MapDataset` in energy: .. GENERATED FROM PYTHON SOURCE LINES 335-343 .. code-block:: Python sliced = dataset_cta.slice_by_energy( energy_min=1 * u.TeV, energy_max=5 * u.TeV, name="slice-energy" ) sliced.counts.plot_grid(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_012.png :alt: Energy 1.00 TeV - 1.58 TeV, Energy 1.58 TeV - 2.51 TeV, Energy 2.51 TeV - 3.98 TeV :srcset: /tutorials/details/images/sphx_glr_datasets_012.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 344-347 The same operation will be applied to all other maps contained in the datasets such as `~gammapy.datasets.MapDataset.mask_fit`: .. GENERATED FROM PYTHON SOURCE LINES 347-352 .. code-block:: Python sliced.mask_fit.plot_grid() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_013.png :alt: Energy 1.00 TeV - 1.58 TeV, Energy 1.58 TeV - 2.51 TeV, Energy 2.51 TeV - 3.98 TeV :srcset: /tutorials/details/images/sphx_glr_datasets_013.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 353-360 Resampling datasets ~~~~~~~~~~~~~~~~~~~ It can often be useful to coarsely rebin an initially computed datasets by a specified factor. This can be done in either spatial or energy axes: .. GENERATED FROM PYTHON SOURCE LINES 360-367 .. code-block:: Python plt.figure() downsampled = dataset_cta.downsample(factor=8) downsampled.counts.sum_over_axes().plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_014.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_014.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 368-370 And the same down-sampling process is possible along the energy axis: .. GENERATED FROM PYTHON SOURCE LINES 370-378 .. code-block:: Python downsampled_energy = dataset_cta.downsample( factor=5, axis_name="energy", name="downsampled-energy" ) downsampled_energy.counts.plot_grid(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_015.png :alt: Energy 100 GeV - 1.00 TeV, Energy 1.00 TeV - 10.0 TeV :srcset: /tutorials/details/images/sphx_glr_datasets_015.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 379-382 In the printout one can see that the actual number of counts is preserved during the down-sampling: .. GENERATED FROM PYTHON SOURCE LINES 382-386 .. code-block:: Python print(downsampled_energy, dataset_cta) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : downsampled-energy Total counts : 104317 Total background counts : 91507.69 Total excess counts : 12809.31 Predicted counts : 91507.69 Predicted background counts : 91507.69 Predicted excess counts : nan Exposure min : 6.28e+07 m2 s Exposure max : 1.90e+10 m2 s Number of total bins : 153600 Number of fit bins : 22608 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 MapDataset ---------- Name : dataset-cta Total counts : 104317 Total background counts : 91507.70 Total excess counts : 12809.31 Predicted counts : 91719.65 Predicted background counts : 91507.69 Predicted excess counts : 211.96 Exposure min : 6.28e+07 m2 s Exposure max : 1.90e+10 m2 s Number of total bins : 768000 Number of fit bins : 67824 Fit statistic type : cash Fit statistic value (-2 log(L)) : 44319.08 Number of models : 2 Number of parameters : 8 Number of free parameters : 5 Component 0: SkyModel Name : gc Datasets names : None Spectral model type : PowerLawSpectralModel Spatial model type : PointSpatialModel Temporal model type : Parameters: index : 2.000 +/- 0.00 amplitude : 1.00e-12 +/- 0.0e+00 1 / (TeV s cm2) reference (frozen): 1.000 TeV lon_0 : 4.650 +/- 0.00 rad lat_0 : -0.505 +/- 0.00 rad Component 1: FoVBackgroundModel Name : dataset-cta-bkg Datasets names : ['dataset-cta'] Spectral model type : PowerLawNormSpectralModel Parameters: tilt (frozen): 0.000 norm : 1.000 +/- 0.00 reference (frozen): 1.000 TeV .. GENERATED FROM PYTHON SOURCE LINES 387-390 We can also resample the finer binned datasets to an arbitrary coarser energy binning using: .. GENERATED FROM PYTHON SOURCE LINES 390-397 .. code-block:: Python energy_axis_new = MapAxis.from_energy_edges([0.1, 0.3, 1, 3, 10] * u.TeV) resampled = dataset_cta.resample_energy_axis(energy_axis=energy_axis_new) resampled.counts.plot_grid(ncols=2, add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_016.png :alt: Energy 100 GeV - 251 GeV, Energy 251 GeV - 1.00 TeV, Energy 1.00 TeV - 2.51 TeV, Energy 2.51 TeV - 10.0 TeV :srcset: /tutorials/details/images/sphx_glr_datasets_016.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 398-401 To squash the whole dataset into a single energy bin there is the `~gammapy.datasets.MapDataset.to_image()` convenience method: .. GENERATED FROM PYTHON SOURCE LINES 401-407 .. code-block:: Python dataset_image = dataset_cta.to_image() dataset_image.counts.plot(add_cbar=True) plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_017.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_017.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 408-418 `~gammapy.datasets.SpectrumDataset` ----------------------------------- `~gammapy.datasets.SpectrumDataset` inherits from a `~gammapy.datasets.MapDataset`, and is specially adapted for 1D spectral analysis, and uses a `~gammapy.maps.RegionGeom` instead of a `~gammapy.maps.WcsGeom`. A `~gammapy.datasets.MapDataset` can be converted to a `~gammapy.datasets.SpectrumDataset`, by summing the ``counts`` and ``background`` inside the ``on_region``, which can then be used for classical spectral analysis. Containment correction is feasible only for circular regions. .. GENERATED FROM PYTHON SOURCE LINES 418-430 .. code-block:: Python region = CircleSkyRegion(SkyCoord(0, 0, unit="deg", frame="galactic"), 0.5 * u.deg) spectrum_dataset = dataset_cta.to_spectrum_dataset( region, containment_correction=True, name="spectrum-dataset" ) # For a quick look spectrum_dataset.peek() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_018.png :alt: Counts, Exposure, Energy Dispersion :srcset: /tutorials/details/images/sphx_glr_datasets_018.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 431-436 A `~gammapy.datasets.MapDataset` can also be integrated over the ``on_region`` to create a `~gammapy.datasets.MapDataset` with a `~gammapy.maps.RegionGeom`. Complex regions can be handled and since the full IRFs are used, containment correction is not required. .. GENERATED FROM PYTHON SOURCE LINES 436-441 .. code-block:: Python reg_dataset = dataset_cta.to_region_map_dataset(region, name="region-map-dataset") print(reg_dataset) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : region-map-dataset Total counts : 5644 Total background counts : 3323.66 Total excess counts : 2320.34 Predicted counts : 3323.66 Predicted background counts : 3323.66 Predicted excess counts : nan Exposure min : 4.66e+08 m2 s Exposure max : 1.89e+10 m2 s Number of total bins : 10 Number of fit bins : 6 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 442-448 `~gammapy.datasets.FluxPointsDataset` ------------------------------------- `~gammapy.datasets.FluxPointsDataset` is a `~gammapy.datasets.Dataset` container for precomputed flux points, which can be then used in fitting. .. GENERATED FROM PYTHON SOURCE LINES 448-459 .. code-block:: Python flux_points = FluxPoints.read( "$GAMMAPY_DATA/tests/spectrum/flux_points/diff_flux_points.fits" ) model = SkyModel(spectral_model=PowerLawSpectralModel(index=2.3)) fp_dataset = FluxPointsDataset(data=flux_points, models=model) fp_dataset.plot_spectrum() plt.show() .. image-sg:: /tutorials/details/images/sphx_glr_datasets_019.png :alt: datasets :srcset: /tutorials/details/images/sphx_glr_datasets_019.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 460-465 The masks on `~gammapy.datasets.FluxPointsDataset` are `~numpy.ndarray` objects and the data is a `~gammapy.estimators.FluxPoints` object. The `~gammapy.datasets.FluxPointsDataset.mask_safe`, by default, masks the upper limit points. .. GENERATED FROM PYTHON SOURCE LINES 465-468 .. code-block:: Python print(fp_dataset.mask_safe) # Note: the mask here is simply a numpy array .. rst-class:: sphx-glr-script-out .. code-block:: none [[[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[ True]] [[False]] [[ True]] [[False]] [[False]] [[ True]] [[False]] [[False]]] .. GENERATED FROM PYTHON SOURCE LINES 469-472 .. code-block:: Python print(fp_dataset.data) # is a `FluxPoints` object .. rst-class:: sphx-glr-script-out .. code-block:: none FluxPoints ---------- geom : RegionGeom axes : ['lon', 'lat', 'energy'] shape : (1, 1, 24) quantities : ['norm', 'norm_err', 'norm_ul'] ref. model : pl n_sigma : 1 n_sigma_ul : 2.0 sqrt_ts_threshold_ul : 2 sed type init : dnde .. GENERATED FROM PYTHON SOURCE LINES 473-477 .. code-block:: Python print(fp_dataset.data_shape()) # number of data points .. rst-class:: sphx-glr-script-out .. code-block:: none (24,) .. GENERATED FROM PYTHON SOURCE LINES 478-481 For an example of fitting `~gammapy.estimators.FluxPoints`, see :doc:`/tutorials/analysis-1d/sed_fitting`, and for using source catalogs see :doc:`/tutorials/details/catalog`. .. GENERATED FROM PYTHON SOURCE LINES 485-505 Datasets -------- `~gammapy.datasets.Datasets` are a collection of `~gammapy.datasets.Dataset` objects. They can be of the same type, or of different types, eg: mix of `~gammapy.datasets.FluxPointsDataset`, `~gammapy.datasets.MapDataset` and `~gammapy.datasets.SpectrumDataset`. For modelling and fitting of a list of `~gammapy.datasets.Dataset` objects, you can either: - (A) Do a joint fitting of all the datasets together **OR** - (B) Stack the datasets together, and then fit them. `~gammapy.datasets.Datasets` is a convenient tool to handle joint fitting of simultaneous datasets. As an example, please see :doc:`/tutorials/analysis-3d/analysis_mwl`. To see how stacking is performed, please see :ref:`stack`. To create a `~gammapy.datasets.Datasets` object, pass a list of `~gammapy.datasets.Dataset` on init, eg .. GENERATED FROM PYTHON SOURCE LINES 505-511 .. code-block:: Python datasets = Datasets([dataset_empty, dataset_cta]) print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : MapDataset Name : dataset-empty Instrument : Models : Dataset 1: Type : MapDataset Name : dataset-cta Instrument : Models : ['gc', 'dataset-cta-bkg'] .. GENERATED FROM PYTHON SOURCE LINES 512-516 If all the datasets have the same type we can also print an info table, collecting all the information from the individual calls to `~gammapy.datasets.Dataset.info_dict()`: .. GENERATED FROM PYTHON SOURCE LINES 516-522 .. code-block:: Python display(datasets.info_table()) # quick info of all datasets print(datasets.names) # unique name of each dataset .. rst-class:: sphx-glr-script-out .. code-block:: none name counts excess ... stat_type stat_sum ... ------------- ------ ------------------ ... --------- ----------------- dataset-empty 0 0.0 ... cash nan dataset-cta 104317 12809.313714614138 ... cash 44319.07871368968 ['dataset-empty', 'dataset-cta'] .. GENERATED FROM PYTHON SOURCE LINES 523-525 We can access individual datasets in `~gammapy.datasets.Datasets` object by name: .. GENERATED FROM PYTHON SOURCE LINES 525-529 .. code-block:: Python print(datasets["dataset-empty"]) # extracts the first dataset .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : dataset-empty Total counts : 0 Total background counts : 0.00 Total excess counts : 0.00 Predicted counts : 0.00 Predicted background counts : 0.00 Predicted excess counts : nan Exposure min : 0.00e+00 m2 s Exposure max : 0.00e+00 m2 s Number of total bins : 110000 Number of fit bins : 0 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 530-532 Or by index: .. GENERATED FROM PYTHON SOURCE LINES 532-536 .. code-block:: Python print(datasets[0]) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDataset ---------- Name : dataset-empty Total counts : 0 Total background counts : 0.00 Total excess counts : 0.00 Predicted counts : 0.00 Predicted background counts : 0.00 Predicted excess counts : nan Exposure min : 0.00e+00 m2 s Exposure max : 0.00e+00 m2 s Number of total bins : 110000 Number of fit bins : 0 Fit statistic type : cash Fit statistic value (-2 log(L)) : nan Number of models : 0 Number of parameters : 0 Number of free parameters : 0 .. GENERATED FROM PYTHON SOURCE LINES 537-539 Other list type operations work as well such as: .. GENERATED FROM PYTHON SOURCE LINES 539-545 .. code-block:: Python # Use python list convention to remove/add datasets, eg: datasets.remove("dataset-empty") print(datasets.names) .. rst-class:: sphx-glr-script-out .. code-block:: none ['dataset-cta'] .. GENERATED FROM PYTHON SOURCE LINES 546-548 Or .. GENERATED FROM PYTHON SOURCE LINES 548-553 .. code-block:: Python datasets.append(spectrum_dataset) print(datasets.names) .. rst-class:: sphx-glr-script-out .. code-block:: none ['dataset-cta', 'spectrum-dataset'] .. GENERATED FROM PYTHON SOURCE LINES 554-557 Let’s create a list of spectrum datasets to illustrate some more functionality: .. GENERATED FROM PYTHON SOURCE LINES 557-569 .. code-block:: Python datasets = Datasets() path = make_path("$GAMMAPY_DATA/joint-crab/spectra/hess") for filename in path.glob("pha_*.fits"): dataset = SpectrumDatasetOnOff.read(filename) datasets.append(dataset) print(datasets) .. rst-class:: sphx-glr-script-out .. code-block:: none Datasets -------- Dataset 0: Type : SpectrumDatasetOnOff Name : 23523 Instrument : Models : Dataset 1: Type : SpectrumDatasetOnOff Name : 23559 Instrument : Models : Dataset 2: Type : SpectrumDatasetOnOff Name : 23592 Instrument : Models : Dataset 3: Type : SpectrumDatasetOnOff Name : 23526 Instrument : Models : .. GENERATED FROM PYTHON SOURCE LINES 570-572 Now we can stack all datasets using `~gammapy.datasets.Datasets.stack_reduce()`: .. GENERATED FROM PYTHON SOURCE LINES 572-577 .. code-block:: Python stacked = datasets.stack_reduce(name="stacked") print(stacked) .. rst-class:: sphx-glr-script-out .. code-block:: none SpectrumDatasetOnOff -------------------- Name : stacked Total counts : 459 Total background counts : 27.50 Total excess counts : 431.50 Predicted counts : 45.27 Predicted background counts : 45.27 Predicted excess counts : nan Exposure min : 2.13e+06 m2 s Exposure max : 2.63e+09 m2 s Number of total bins : 80 Number of fit bins : 43 Fit statistic type : wstat Fit statistic value (-2 log(L)) : 1530.36 Number of models : 0 Number of parameters : 0 Number of free parameters : 0 Total counts_off : 645 Acceptance : 168 Acceptance off : 4497 .. GENERATED FROM PYTHON SOURCE LINES 578-580 Or slice all datasets by a given energy range: .. GENERATED FROM PYTHON SOURCE LINES 580-583 .. code-block:: Python datasets_sliced = datasets.slice_by_energy(energy_min="1 TeV", energy_max="10 TeV") print(datasets_sliced.energy_ranges) .. rst-class:: sphx-glr-script-out .. code-block:: none (, ) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.059 seconds) .. _sphx_glr_download_tutorials_details_datasets.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/details/datasets.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: datasets.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: datasets.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: datasets.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_