.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/api/maps.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_api_maps.py: Maps ==== A thorough tutorial to work with WCS maps. .. figure:: ../../_static/gammapy_maps.png :alt: Gammapy Maps Illustration Gammapy Maps Illustration Introduction ------------ The `~gammapy.maps` submodule contains classes for representing pixilised data on the sky with an arbitrary number of non-spatial dimensions such as energy, time, event class or any possible user-defined dimension (illustrated in the image above). The main `Map` data structure features a uniform API for `WCS `__ as well as `HEALPix `__ based images. The API also generalizes simple image based operations such as smoothing, interpolation and reprojection to the arbitrary extra dimensions and makes working with (2 + N)-dimensional hypercubes as easy as working with a simple 2D image. Further information is also provided on the `~gammapy.maps` docs page. In the following introduction we will learn all the basics of working with WCS based maps. HEALPix based maps will be covered in a future tutorial. Make sure you have worked through the :doc:`Gammapy overview `, because a solid knowledge about working with `SkyCoord` and `Quantity` objects as well as `Numpy `__ is required for this tutorial. This notebook is rather lengthy, but getting to know the `Map` data structure in detail is essential for working with Gammapy and will allow you to fulfill complex analysis tasks with very few and simple code in future! .. GENERATED FROM PYTHON SOURCE LINES 43-46 Setup ----- .. GENERATED FROM PYTHON SOURCE LINES 46-70 .. code-block:: Python import os # %matplotlib inline import numpy as np from astropy import units as u from astropy.convolution import convolve from astropy.coordinates import SkyCoord from astropy.io import fits from astropy.table import Table from astropy.time import Time import matplotlib.pyplot as plt from IPython.display import display from gammapy.data import EventList from gammapy.maps import ( LabelMapAxis, Map, MapAxes, MapAxis, TimeMapAxis, WcsGeom, WcsNDMap, ) .. GENERATED FROM PYTHON SOURCE LINES 71-73 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 73-78 .. code-block:: Python from gammapy.utils.check import check_tutorials_setup 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.dev241+g0271bebfc 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.20.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev .. GENERATED FROM PYTHON SOURCE LINES 79-88 Creating WCS Maps ----------------- Using Factory Methods ~~~~~~~~~~~~~~~~~~~~~ Maps are most easily created using the `~gammapy.maps.Map.create` factory method: .. GENERATED FROM PYTHON SOURCE LINES 88-92 .. code-block:: Python m_allsky = Map.create() .. GENERATED FROM PYTHON SOURCE LINES 93-98 Calling `~gammapy.maps.Map.create` without any further arguments creates by default an allsky WCS map using a CAR projection, ICRS coordinates and a pixel size of 1 deg. This can be easily checked by printing the `~gammapy.maps.Map.geom` attribute of the map: .. GENERATED FROM PYTHON SOURCE LINES 98-102 .. code-block:: Python print(m_allsky.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat'] shape : (3600, 1800) ndim : 2 frame : icrs projection : CAR center : 0.0 deg, 0.0 deg width : 360.0 deg x 180.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 103-112 The `~gammapy.maps.Map.geom` attribute is a `~gammapy.maps.Geom` object, that defines the basic geometry of the map, such as size of the pixels, width and height of the image, coordinate system etc., but we will learn more about this object later. Besides the ``.geom`` attribute the map has also a ``.data`` attribute, which is just a plain ``~numpy.ndarray`` and stores the data associated with this map: .. GENERATED FROM PYTHON SOURCE LINES 112-116 .. code-block:: Python print(m_allsky.data) .. rst-class:: sphx-glr-script-out .. code-block:: none [[0. 0. 0. ... 0. 0. 0.] [0. 0. 0. ... 0. 0. 0.] [0. 0. 0. ... 0. 0. 0.] ... [0. 0. 0. ... 0. 0. 0.] [0. 0. 0. ... 0. 0. 0.] [0. 0. 0. ... 0. 0. 0.]] .. GENERATED FROM PYTHON SOURCE LINES 117-122 By default maps are filled with zeros. The ``map_type`` argument can be used to control the pixelization scheme (WCS or HPX). .. GENERATED FROM PYTHON SOURCE LINES 122-132 .. code-block:: Python position = SkyCoord(0.0, 5.0, frame="galactic", unit="deg") # Create a WCS Map m_wcs = Map.create(binsz=0.1, map_type="wcs", skydir=position, width=10.0) # Create a HPX Map m_hpx = Map.create(binsz=0.1, map_type="hpx", skydir=position, width=10.0) .. GENERATED FROM PYTHON SOURCE LINES 133-136 Here is an example that creates a WCS map centered on the Galactic center and now uses Galactic coordinates: .. GENERATED FROM PYTHON SOURCE LINES 136-144 .. code-block:: Python skydir = SkyCoord(0, 0, frame="galactic", unit="deg") m_gc = Map.create( binsz=0.02, width=(10, 5), skydir=skydir, frame="galactic", proj="TAN" ) print(m_gc.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 frame : galactic projection : TAN center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 145-150 In addition we have defined a TAN projection, a pixel size of ``0.02`` deg and a width of the map of ``10 deg x 5 deg``. The `width` argument also takes scalar value instead of a tuple, which is interpreted as both the width and height of the map, so that a quadratic map is created. .. GENERATED FROM PYTHON SOURCE LINES 153-164 Creating from a Map Geometry ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ As we have seen in the first examples, the `~gammapy.maps.Map` object couples the data (stored as a `~numpy.ndarray`) with a `~gammapy.maps.Geom` object. The `~gammapy.maps.~Geom` object can be seen as a generalization of an `astropy.wcs.WCS` object, providing the information on how the data maps to physical coordinate systems. In some cases e.g. when creating many maps with the same WCS geometry it can be advantageous to first create the map geometry independent of the map object it-self: .. GENERATED FROM PYTHON SOURCE LINES 164-168 .. code-block:: Python wcs_geom = WcsGeom.create(binsz=0.02, width=(10, 5), skydir=(0, 0), frame="galactic") .. GENERATED FROM PYTHON SOURCE LINES 169-172 And then create the map objects from the ``wcs_geom`` geometry specification: .. GENERATED FROM PYTHON SOURCE LINES 172-179 .. code-block:: Python maps = {} for name in ["counts", "background"]: maps[name] = Map.from_geom(wcs_geom) .. GENERATED FROM PYTHON SOURCE LINES 180-183 The `~gammapy.maps.Geom` object also has a few helpful methods. E.g. we can check whether a given position on the sky is contained in the map geometry: .. GENERATED FROM PYTHON SOURCE LINES 183-189 .. code-block:: Python # define the position of the Galactic center and anti-center positions = SkyCoord([0, 180], [0, 0], frame="galactic", unit="deg") wcs_geom.contains(positions) .. rst-class:: sphx-glr-script-out .. code-block:: none array([ True, False]) .. GENERATED FROM PYTHON SOURCE LINES 190-192 Or get the image center of the map: .. GENERATED FROM PYTHON SOURCE LINES 192-196 .. code-block:: Python print(wcs_geom.center_skydir) .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 197-199 Or we can also retrieve the solid angle per pixel of the map: .. GENERATED FROM PYTHON SOURCE LINES 199-203 .. code-block:: Python print(wcs_geom.solid_angle()) .. rst-class:: sphx-glr-script-out .. code-block:: none [[1.21731921e-07 1.21731921e-07 1.21731921e-07 ... 1.21731921e-07 1.21731921e-07 1.21731921e-07] [1.21733761e-07 1.21733761e-07 1.21733761e-07 ... 1.21733761e-07 1.21733761e-07 1.21733761e-07] [1.21735587e-07 1.21735587e-07 1.21735587e-07 ... 1.21735587e-07 1.21735587e-07 1.21735587e-07] ... [1.21735587e-07 1.21735587e-07 1.21735587e-07 ... 1.21735587e-07 1.21735587e-07 1.21735587e-07] [1.21733761e-07 1.21733761e-07 1.21733761e-07 ... 1.21733761e-07 1.21733761e-07 1.21733761e-07] [1.21731921e-07 1.21731921e-07 1.21731921e-07 ... 1.21731921e-07 1.21731921e-07 1.21731921e-07]] sr .. GENERATED FROM PYTHON SOURCE LINES 204-212 Adding Non-Spatial Axes ~~~~~~~~~~~~~~~~~~~~~~~ In many analysis scenarios we would like to add extra dimension to the maps to study e.g. energy or time dependency of the data. Those non-spatial dimensions are handled with the `~gammapy.maps.MapAxis` object. Let us first define an energy axis, with 4 bins: .. GENERATED FROM PYTHON SOURCE LINES 212-219 .. code-block:: Python energy_axis = MapAxis.from_bounds( 1, 100, nbin=4, unit="TeV", name="energy", interp="log" ) print(energy_axis) .. rst-class:: sphx-glr-script-out .. code-block:: none MapAxis name : energy unit : 'TeV' nbins : 4 node type : edges edges min : 1.0e+00 TeV edges max : 1.0e+02 TeV interp : log .. GENERATED FROM PYTHON SOURCE LINES 220-225 Where ``interp='log'`` specifies that a logarithmic spacing is used between the bins, equivalent to ``np.logspace(0, 2, 4)``. This `~gammapy.maps.MapAxis` object we can now pass to `~gammapy.maps.Map.create()` using the ``axes=`` argument: .. GENERATED FROM PYTHON SOURCE LINES 225-230 .. code-block:: Python m_cube = Map.create(binsz=0.02, width=(10, 5), frame="galactic", axes=[energy_axis]) print(m_cube.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat', 'energy'] shape : (500, 250, 4) ndim : 3 frame : galactic projection : CAR center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 231-239 Now we see that besides ``lon`` and ``lat`` the map has an additional axes named ``energy`` with 4 bins. The total dimension of the map is now ``ndim=3``. We can also add further axes by passing a list of `~gammapy.maps.MapAxis` objects. To demonstrate this we create a time axis with linearly spaced bins and pass both axes to `Map.create()`: .. GENERATED FROM PYTHON SOURCE LINES 239-248 .. code-block:: Python time_axis = MapAxis.from_bounds(0, 24, nbin=24, unit="hour", name="time", interp="lin") m_4d = Map.create( binsz=0.02, width=(10, 5), frame="galactic", axes=[energy_axis, time_axis] ) print(m_4d.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat', 'energy', 'time'] shape : (500, 250, 4, 24) ndim : 4 frame : galactic projection : CAR center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 249-259 The `~gammapy.maps.MapAxis` object internally stores the coordinates or “position values” associated with every map axis bin or “node”. We distinguish between two node types: ``"edges"`` and ``"center"``. The node type ``"edges"``\ (which is also the default) specifies that the data associated with this axis is integrated between the edges of the bin (e.g. counts data). The node type ``"center"`` specifies that the data is given at the center of the bin (e.g. exposure or differential fluxes). The edges of the bins can be checked with `~gammapy.maps.MapAxis.edges` attribute: .. GENERATED FROM PYTHON SOURCE LINES 259-263 .. code-block:: Python print(energy_axis.edges) .. rst-class:: sphx-glr-script-out .. code-block:: none [ 1. 3.16227766 10. 31.6227766 100. ] TeV .. GENERATED FROM PYTHON SOURCE LINES 264-267 The numbers are given in the units we specified above, which can be checked again with: .. GENERATED FROM PYTHON SOURCE LINES 267-271 .. code-block:: Python print(energy_axis.unit) .. rst-class:: sphx-glr-script-out .. code-block:: none TeV .. GENERATED FROM PYTHON SOURCE LINES 272-275 The centers of the axis bins can be checked with the `~gammapy.maps.MapAxis.center` attribute: .. GENERATED FROM PYTHON SOURCE LINES 275-278 .. code-block:: Python print(energy_axis.center) .. rst-class:: sphx-glr-script-out .. code-block:: none [ 1.77827941 5.62341325 17.7827941 56.23413252] TeV .. GENERATED FROM PYTHON SOURCE LINES 279-285 Adding Non-contiguous axes ~~~~~~~~~~~~~~~~~~~~~~~~~~ Non-spatial map axes can also be handled through two other objects known as the `~gammapy.maps.TimeMapAxis` and the `~gammapy.maps.LabelMapAxis`. .. GENERATED FROM PYTHON SOURCE LINES 288-294 TimeMapAxis ^^^^^^^^^^^ The `~gammapy.maps.TimeMapAxis` object provides an axis for non-adjacent time intervals. .. GENERATED FROM PYTHON SOURCE LINES 294-303 .. code-block:: Python time_map_axis = TimeMapAxis( edges_min=[1, 5, 10, 15] * u.day, edges_max=[2, 7, 13, 18] * u.day, reference_time=Time("2020-03-19"), ) print(time_map_axis) .. rst-class:: sphx-glr-script-out .. code-block:: none TimeMapAxis ----------- name : time nbins : 4 reference time : 2020-03-19 00:00:00.000 scale : utc time min. : 2020-03-20 00:00:00.000 time max. : 2020-04-06 00:00:00.000 total time : 216.0 h .. GENERATED FROM PYTHON SOURCE LINES 304-307 This ``time_map_axis`` can then be utilised in a similar way to the previous implementation to create a `~gammapy.maps.Map`. .. GENERATED FROM PYTHON SOURCE LINES 307-313 .. code-block:: Python map_4d = Map.create( binsz=0.02, width=(10, 5), frame="galactic", axes=[energy_axis, time_map_axis] ) print(map_4d.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat', 'energy', 'time'] shape : (500, 250, 4, 4) ndim : 4 frame : galactic projection : CAR center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 314-318 It is possible to utilise the `~gammapy.maps.TimeMapAxis.slice` attrribute to create new a `~gammapy.maps.TimeMapAxis`. Here we are slicing between the first and third axis to extract the subsection of the axis between indice 0 and 2. .. GENERATED FROM PYTHON SOURCE LINES 318-321 .. code-block:: Python print(time_map_axis.slice([0, 2])) .. rst-class:: sphx-glr-script-out .. code-block:: none TimeMapAxis ----------- name : time nbins : 2 reference time : 2020-03-19 00:00:00.000 scale : utc time min. : 2020-03-20 00:00:00.000 time max. : 2020-04-01 00:00:00.000 total time : 96.0 h .. GENERATED FROM PYTHON SOURCE LINES 322-325 It is also possible to `~gammapy.maps.TimeMapAxis.squash` the axis, which squashes the existing axis into one bin. This creates a new axis between the extreme edges of the initial axis. .. GENERATED FROM PYTHON SOURCE LINES 325-329 .. code-block:: Python print(time_map_axis.squash()) .. rst-class:: sphx-glr-script-out .. code-block:: none TimeMapAxis ----------- name : time nbins : 1 reference time : 2020-03-19 00:00:00.000 scale : utc time min. : 2020-03-20 00:00:00.000 time max. : 2020-04-06 00:00:00.000 total time : 408.0 h .. GENERATED FROM PYTHON SOURCE LINES 330-332 The `~gammapy.maps.TimeMapAxis.is_contiguous` method returns a boolean which indicates whether the `~gammapy.maps.TimeMapAxis` is contiguous or not. .. GENERATED FROM PYTHON SOURCE LINES 332-335 .. code-block:: Python print(time_map_axis.is_contiguous) .. rst-class:: sphx-glr-script-out .. code-block:: none False .. GENERATED FROM PYTHON SOURCE LINES 336-339 As we have a non-contiguous axis we can print the array of bin edges for both the minimum axis edges (`~gammapy.maps.TimeMapAxis.edges_min`) and the maximum axis edges (`~gammapy.maps.TimeMapAxis.edges_max`). .. GENERATED FROM PYTHON SOURCE LINES 339-344 .. code-block:: Python print(time_map_axis.edges_min) print(time_map_axis.edges_max) .. rst-class:: sphx-glr-script-out .. code-block:: none [ 1. 5. 10. 15.] d [ 2. 7. 13. 18.] d .. GENERATED FROM PYTHON SOURCE LINES 345-348 Next, we use the `~gammapy.maps.TimeMapAxis.to_contiguous` functionality to create a contiguous axis and expose `~gammapy.maps.TimeMapAxis.edges`. This method returns a `~astropy.units.Quantity` with respect to the reference time. .. GENERATED FROM PYTHON SOURCE LINES 348-356 .. code-block:: Python time_map_axis_contiguous = time_map_axis.to_contiguous() print(time_map_axis_contiguous.is_contiguous) print(time_map_axis_contiguous.edges) .. rst-class:: sphx-glr-script-out .. code-block:: none True [ 1. 2. 5. 7. 10. 13. 15. 18.] d .. GENERATED FROM PYTHON SOURCE LINES 357-358 The `~gammapy.maps.TimeMapAxis.time_edges` will return the `~astropy.time.Time` object directly .. GENERATED FROM PYTHON SOURCE LINES 358-362 .. code-block:: Python print(time_map_axis_contiguous.time_edges) .. rst-class:: sphx-glr-script-out .. code-block:: none ['2020-03-20 00:00:00.000' '2020-03-21 00:00:00.000' '2020-03-24 00:00:00.000' '2020-03-26 00:00:00.000' '2020-03-29 00:00:00.000' '2020-04-01 00:00:00.000' '2020-04-03 00:00:00.000' '2020-04-06 00:00:00.000'] .. GENERATED FROM PYTHON SOURCE LINES 363-370 `~gammapy.maps.TimeMapAxis` also has both functionalities of `~gammapy.maps.TimeMapAxis.coord_to_pix` and `~gammapy.maps.TimeMapAxis.coord_to_idx`. The `~gammapy.maps.TimeMapAxis.coord_to_idx` attribute will give the index of the ``time`` specified, similarly for `~gammapy.maps.TimeMapAxis.coord_to_pix` which returns the pixel. A linear interpolation is assumed. Start by choosing a time which we know is within the `~gammapy.maps.TimeMapAxis` and see the results. .. GENERATED FROM PYTHON SOURCE LINES 370-379 .. code-block:: Python time = Time(time_map_axis.time_max.mjd[0], format="mjd") print(time_map_axis.coord_to_pix(time)) print(time_map_axis.coord_to_idx(time)) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.5] 0 .. GENERATED FROM PYTHON SOURCE LINES 380-381 This functionality can also be used with an array of `~astropy.time.Time` values. .. GENERATED FROM PYTHON SOURCE LINES 381-388 .. code-block:: Python times = Time(time_map_axis.time_max.mjd, format="mjd") print(time_map_axis.coord_to_pix(times)) print(time_map_axis.coord_to_idx(times)) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.5 1.5 2.5 3.5] [0 1 2 3] .. GENERATED FROM PYTHON SOURCE LINES 389-395 Note here we take a `~astropy.time.Time` which is outside the edges. A linear interpolation is assumed for both methods, therefore for a time outside the ``time_map_axis`` there is no extrapolation and -1 is returned. Note: due to this, the `~gammapy.maps.TimeMapAxis.coord_to_pix` method will return ``nan`` and the `~gammapy.maps.TimeMapAxis.coord_to_idx` method returns -1. .. GENERATED FROM PYTHON SOURCE LINES 395-401 .. code-block:: Python print(time_map_axis.coord_to_pix(Time(time.mjd + 1, format="mjd"))) print(time_map_axis.coord_to_idx(Time(time.mjd + 1, format="mjd"))) .. rst-class:: sphx-glr-script-out .. code-block:: none [nan] -1 .. GENERATED FROM PYTHON SOURCE LINES 402-408 LabelMapAxis ^^^^^^^^^^^^ The `~gammapy.maps.LabelMapAxis` object allows for handling of labels for map axes. It provides an axis for non-numeric entries. .. GENERATED FROM PYTHON SOURCE LINES 408-415 .. code-block:: Python label_axis = LabelMapAxis( labels=["dataset-1", "dataset-2", "dataset-3"], name="dataset" ) print(label_axis) .. rst-class:: sphx-glr-script-out .. code-block:: none LabelMapAxis ------------ name : dataset nbins : 3 node type : label labels : ['dataset-1', 'dataset-2', 'dataset-3'] .. GENERATED FROM PYTHON SOURCE LINES 416-417 The labels can be checked using the `~gammapy.maps.LabelMapAxis.center` attribute: .. GENERATED FROM PYTHON SOURCE LINES 417-420 .. code-block:: Python print(label_axis.center) .. rst-class:: sphx-glr-script-out .. code-block:: none ['dataset-1' 'dataset-2' 'dataset-3'] .. GENERATED FROM PYTHON SOURCE LINES 421-422 To obtain the position of the label, one can utilise the `~gammapy.maps.LabelMapAxis.coord_to_pix` attribute .. GENERATED FROM PYTHON SOURCE LINES 422-425 .. code-block:: Python print(label_axis.coord_to_pix(["dataset-3"])) .. rst-class:: sphx-glr-script-out .. code-block:: none [2.] .. GENERATED FROM PYTHON SOURCE LINES 426-431 To adapt and create new axes the following attributes can be utilised: `~gammapy.maps.LabelMapAxis.concatenate`, `~gammapy.maps.LabelMapAxis.slice` and `~gammapy.maps.LabelMapAxis.squash`. Combining two different `~gammapy.maps.LabelMapAxis` is done in the following way: .. GENERATED FROM PYTHON SOURCE LINES 431-436 .. code-block:: Python label_axis2 = LabelMapAxis(labels=["dataset-a", "dataset-b"], name="dataset") print(label_axis.concatenate(label_axis2)) .. rst-class:: sphx-glr-script-out .. code-block:: none LabelMapAxis ------------ name : dataset nbins : 5 node type : label labels : ['dataset-1', 'dataset-2', 'dataset-3', 'dataset-a', 'dataset-b'] .. GENERATED FROM PYTHON SOURCE LINES 437-439 A new `~gammapy.maps.LabelMapAxis` can be created by slicing an already existing one. Here we are slicing between the second and third bins to extract the subsection. .. GENERATED FROM PYTHON SOURCE LINES 439-441 .. code-block:: Python print(label_axis.slice([1, 2])) .. rst-class:: sphx-glr-script-out .. code-block:: none LabelMapAxis ------------ name : dataset nbins : 2 node type : label labels : ['dataset-2', 'dataset-3'] .. GENERATED FROM PYTHON SOURCE LINES 442-443 A new axis object can be created by squashing the axis into a single bin. .. GENERATED FROM PYTHON SOURCE LINES 443-447 .. code-block:: Python print(label_axis.squash()) .. rst-class:: sphx-glr-script-out .. code-block:: none LabelMapAxis ------------ name : dataset nbins : 1 node type : label labels : ['dataset-1...dataset-3'] .. GENERATED FROM PYTHON SOURCE LINES 448-452 Mixing the three previous axis types (`~gammapy.maps.MapAxis`, `~gammapy.maps.TimeMapAxis` and `~gammapy.maps.LabelMapAxis`) would be done like so: .. GENERATED FROM PYTHON SOURCE LINES 452-459 .. code-block:: Python axes = MapAxes(axes=[energy_axis, time_map_axis, label_axis]) hdu = axes.to_table_hdu(format="gadf") table = Table.read(hdu) display(table) .. rst-class:: sphx-glr-script-out .. code-block:: none CHANNEL ENERGY E_MIN ... TIME_MIN TIME_MAX DATASET TeV TeV ... d d ------- ------------------ ------------------ ... -------- -------- --------- 0 1.778279410038923 1.0 ... 1.0 2.0 dataset-1 1 1.778279410038923 1.0 ... 1.0 2.0 dataset-2 2 1.778279410038923 1.0 ... 1.0 2.0 dataset-3 3 5.623413251903492 3.1622776601683795 ... 1.0 2.0 dataset-1 4 5.623413251903492 3.1622776601683795 ... 1.0 2.0 dataset-2 5 5.623413251903492 3.1622776601683795 ... 1.0 2.0 dataset-3 6 17.78279410038923 10.000000000000002 ... 1.0 2.0 dataset-1 7 17.78279410038923 10.000000000000002 ... 1.0 2.0 dataset-2 8 17.78279410038923 10.000000000000002 ... 1.0 2.0 dataset-3 9 56.234132519034915 31.62277660168379 ... 1.0 2.0 dataset-1 ... ... ... ... ... ... ... 38 1.778279410038923 1.0 ... 15.0 18.0 dataset-3 39 5.623413251903492 3.1622776601683795 ... 15.0 18.0 dataset-1 40 5.623413251903492 3.1622776601683795 ... 15.0 18.0 dataset-2 41 5.623413251903492 3.1622776601683795 ... 15.0 18.0 dataset-3 42 17.78279410038923 10.000000000000002 ... 15.0 18.0 dataset-1 43 17.78279410038923 10.000000000000002 ... 15.0 18.0 dataset-2 44 17.78279410038923 10.000000000000002 ... 15.0 18.0 dataset-3 45 56.234132519034915 31.62277660168379 ... 15.0 18.0 dataset-1 46 56.234132519034915 31.62277660168379 ... 15.0 18.0 dataset-2 47 56.234132519034915 31.62277660168379 ... 15.0 18.0 dataset-3 Length = 48 rows .. GENERATED FROM PYTHON SOURCE LINES 460-491 Reading and Writing ------------------- Gammapy `~gammapy.maps.Map` objects are serialized using the Flexible Image Transport Format (FITS). Depending on the pixelization scheme (HEALPix or WCS) and presence of non-spatial dimensions the actual convention to write the FITS file is different. By default Gammapy uses a generic convention named ``"gadf"``, which will support WCS and HEALPix formats as well as an arbitrary number of non-spatial axes. The convention is documented in detail on the `Gamma Astro Data Formats `__ page. Other conventions required by specific software (e.g. the Fermi Science Tools) are supported as well. At the moment those are the following - ``"fgst-ccube"``: Fermi counts cube format. - ``"fgst-ltcube"``: Fermi livetime cube format. - ``"fgst-bexpcube"``: Fermi exposure cube format - ``"fgst-template"``: Fermi Galactic diffuse and source template format. - ``"fgst-srcmap"`` and ``"fgst-srcmap-sparse"``: Fermi source map and sparse source map format. The conventions listed above only support an additional energy axis. Reading Maps ~~~~~~~~~~~~ Reading FITS files is mainly exposed via the `~gammapy.maps.Map.read()` method. Let us take a look at a first example: .. GENERATED FROM PYTHON SOURCE LINES 491-497 .. code-block:: Python filename = "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz" m_3fhl_gc = Map.read(filename) print(m_3fhl_gc) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (400, 200) ndim : 2 unit : dtype : >i8 .. GENERATED FROM PYTHON SOURCE LINES 498-506 If ``map_type`` argument is not given when calling read a map object will be instantiated with the pixelization of the input HDU. By default ``Map.read()`` will try to find the first valid data hdu in the filename and read the data from there. If multiple HDUs are present in the FITS file, the desired one can be chosen with the additional `hdu=` argument: .. GENERATED FROM PYTHON SOURCE LINES 506-511 .. code-block:: Python m_3fhl_gc = Map.read(filename, hdu="PRIMARY") print(m_3fhl_gc) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (400, 200) ndim : 2 unit : dtype : >i8 .. GENERATED FROM PYTHON SOURCE LINES 512-517 In rare cases e.g. when the FITS file is not valid or meta data is missing from the header it can be necessary to modify the header of a certain HDU before creating the `Map` object. In this case we can use `astropy.io.fits` directly to read the FITS file: .. GENERATED FROM PYTHON SOURCE LINES 517-523 .. code-block:: Python filename = os.environ["GAMMAPY_DATA"] + "/fermi-3fhl-gc/fermi-3fhl-gc-exposure.fits.gz" hdulist = fits.open(filename) print(hdulist.info()) .. rst-class:: sphx-glr-script-out .. code-block:: none Filename: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/dev/fermi-3fhl-gc/fermi-3fhl-gc-exposure.fits.gz No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 23 (400, 200) float32 None .. GENERATED FROM PYTHON SOURCE LINES 524-527 And then modify the header keyword and use `Map.from_hdulist()` to create the `Map` object after: .. GENERATED FROM PYTHON SOURCE LINES 527-532 .. code-block:: Python hdulist["PRIMARY"].header["BUNIT"] = "cm2 s" print(Map.from_hdulist(hdulist=hdulist)) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (400, 200) ndim : 2 unit : cm2 s dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 533-539 Writing Maps ~~~~~~~~~~~~ Writing FITS files on disk via the `Map.write()` method. Here is a first example: .. GENERATED FROM PYTHON SOURCE LINES 539-543 .. code-block:: Python m_cube.write("example_cube.fits", overwrite=True) .. GENERATED FROM PYTHON SOURCE LINES 544-548 By default Gammapy does not overwrite files. In this example we set `overwrite=True` in case the cell gets executed multiple times. Now we can read back the cube from disk using `Map.read()`: .. GENERATED FROM PYTHON SOURCE LINES 548-553 .. code-block:: Python m_cube = Map.read("example_cube.fits") print(m_cube) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (500, 250, 4) ndim : 3 unit : dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 554-557 We can also choose a different FITS convention to write the example cube in a format compatible to the Fermi Galactic diffuse background model: .. GENERATED FROM PYTHON SOURCE LINES 557-561 .. code-block:: Python m_cube.write("example_cube_fgst.fits", format="fgst-template", overwrite=True) .. GENERATED FROM PYTHON SOURCE LINES 562-565 To understand a little bit better the generic `gadf` convention we use `Map.to_hdulist()` to generate a list of FITS HDUs first: .. GENERATED FROM PYTHON SOURCE LINES 565-570 .. code-block:: Python hdulist = m_4d.to_hdulist(format="gadf") print(hdulist.info()) .. rst-class:: sphx-glr-script-out .. code-block:: none Filename: (No file associated with this HDUList) No. Name Ver Type Cards Dimensions Format 0 PRIMARY 1 PrimaryHDU 30 (500, 250, 4, 24) float32 1 PRIMARY_BANDS 1 BinTableHDU 33 96R x 7C ['K', 'D', 'D', 'D', 'D', 'D', 'D'] None .. GENERATED FROM PYTHON SOURCE LINES 571-575 As we can see the `HDUList` object contains to HDUs. The first one named `PRIMARY` contains the data array with shape corresponding to our data and the WCS information stored in the header: .. GENERATED FROM PYTHON SOURCE LINES 575-579 .. code-block:: Python print(hdulist["PRIMARY"].header) .. rst-class:: sphx-glr-script-out .. code-block:: none SIMPLE = T / conforms to FITS standard BITPIX = -32 / array data type NAXIS = 4 / number of array dimensions NAXIS1 = 500 NAXIS2 = 250 NAXIS3 = 4 NAXIS4 = 24 EXTEND = T WCSAXES = 2 / Number of coordinate axes CRPIX1 = 250.5 / Pixel coordinate of reference point CRPIX2 = 125.5 / Pixel coordinate of reference point CDELT1 = -0.02 / [deg] Coordinate increment at reference point CDELT2 = 0.02 / [deg] Coordinate increment at reference point CUNIT1 = 'deg' / Units of coordinate increment and value CUNIT2 = 'deg' / Units of coordinate increment and value CTYPE1 = 'GLON-CAR' / galactic longitude, plate caree projection CTYPE2 = 'GLAT-CAR' / galactic latitude, plate caree projection CRVAL1 = 0.0 / [deg] Coordinate value at reference point CRVAL2 = 0.0 / [deg] Coordinate value at reference point LONPOLE = 0.0 / [deg] Native longitude of celestial pole LATPOLE = 90.0 / [deg] Native latitude of celestial pole MJDREF = 0.0 / [d] MJD of fiducial time AXCOLS1 = 'E_MIN,E_MAX' INTERP1 = 'log ' AXCOLS2 = 'TIME_MIN,TIME_MAX' INTERP2 = 'lin ' WCSSHAPE= '(500,250,4,24)' BANDSHDU= 'PRIMARY_BANDS' META = '{} ' BUNIT = '' END .. GENERATED FROM PYTHON SOURCE LINES 580-585 The second HDU is a `BinTableHDU` named `PRIMARY_BANDS` contains the information on the non-spatial axes such as name, order, unit, min, max and center values of the axis bins. We use an `astropy.table.Table` to show the information: .. GENERATED FROM PYTHON SOURCE LINES 585-589 .. code-block:: Python print(Table.read(hdulist["PRIMARY_BANDS"])) .. rst-class:: sphx-glr-script-out .. code-block:: none CHANNEL ENERGY E_MIN ... TIME TIME_MIN TIME_MAX TeV TeV ... h h h ------- ------------------ ------------------ ... ---- -------- -------- 0 1.778279410038923 1.0 ... 0.5 0.0 1.0 1 5.623413251903492 3.1622776601683795 ... 0.5 0.0 1.0 2 17.78279410038923 10.000000000000002 ... 0.5 0.0 1.0 3 56.234132519034915 31.62277660168379 ... 0.5 0.0 1.0 4 1.778279410038923 1.0 ... 1.5 1.0 2.0 5 5.623413251903492 3.1622776601683795 ... 1.5 1.0 2.0 6 17.78279410038923 10.000000000000002 ... 1.5 1.0 2.0 7 56.234132519034915 31.62277660168379 ... 1.5 1.0 2.0 8 1.778279410038923 1.0 ... 2.5 2.0 3.0 9 5.623413251903492 3.1622776601683795 ... 2.5 2.0 3.0 ... ... ... ... ... ... ... 86 17.78279410038923 10.000000000000002 ... 21.5 21.0 22.0 87 56.234132519034915 31.62277660168379 ... 21.5 21.0 22.0 88 1.778279410038923 1.0 ... 22.5 22.0 23.0 89 5.623413251903492 3.1622776601683795 ... 22.5 22.0 23.0 90 17.78279410038923 10.000000000000002 ... 22.5 22.0 23.0 91 56.234132519034915 31.62277660168379 ... 22.5 22.0 23.0 92 1.778279410038923 1.0 ... 23.5 23.0 24.0 93 5.623413251903492 3.1622776601683795 ... 23.5 23.0 24.0 94 17.78279410038923 10.000000000000002 ... 23.5 23.0 24.0 95 56.234132519034915 31.62277660168379 ... 23.5 23.0 24.0 Length = 96 rows .. GENERATED FROM PYTHON SOURCE LINES 590-594 Maps can be serialized to a sparse data format by calling write with `sparse=True`. This will write all non-zero pixels in the map to a data table appropriate to the pixelization scheme. .. GENERATED FROM PYTHON SOURCE LINES 594-600 .. code-block:: Python m = Map.create(binsz=0.1, map_type="wcs", width=10.0) m.write("file.fits", hdu="IMAGE", sparse=True, overwrite=True) m = Map.read("file.fits", hdu="IMAGE", map_type="wcs") .. GENERATED FROM PYTHON SOURCE LINES 601-619 Accessing Data -------------- How to get data values ~~~~~~~~~~~~~~~~~~~~~~ All map objects have a set of accessor methods, which can be used to access or update the contents of the map irrespective of its underlying representation. Those accessor methods accept as their first argument a coordinate `tuple` containing scalars, `list`, or `numpy.ndarray` with one tuple element for each dimension. Some methods additionally accept a `dict` or `MapCoord` argument, of which both allow to assign coordinates by axis name. Let us first begin with the `~gammapy.maps.Map.get_by_idx()` method, that accepts a tuple of indices. The order of the indices corresponds to the axis order of the map: .. GENERATED FROM PYTHON SOURCE LINES 619-623 .. code-block:: Python print(m_gc.get_by_idx((50, 30))) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.] .. GENERATED FROM PYTHON SOURCE LINES 624-628 **Important:** Gammapy uses a reversed index order in the map API with the longitude axes first. To achieve the same by directly indexing into the numpy array we have to call: .. GENERATED FROM PYTHON SOURCE LINES 628-632 .. code-block:: Python print(m_gc.data[([30], [50])]) .. rst-class:: sphx-glr-script-out .. code-block:: none [0.] .. GENERATED FROM PYTHON SOURCE LINES 633-636 To check the order of the axes you can always print the ``.geom``` attribute: .. GENERATED FROM PYTHON SOURCE LINES 636-640 .. code-block:: Python print(m_gc.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 frame : galactic projection : TAN center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 641-645 To access values directly by sky coordinates we can use the `~gammapy.maps.Map.get_by_coord()` method. This time we pass in a `dict`, specifying the axes names corresponding to the given coordinates: .. GENERATED FROM PYTHON SOURCE LINES 645-649 .. code-block:: Python print(m_gc.get_by_coord({"lon": [0, 180], "lat": [0, 0]})) .. rst-class:: sphx-glr-script-out .. code-block:: none [ 0. nan] .. GENERATED FROM PYTHON SOURCE LINES 650-659 The units of the coordinates are assumed to be in degrees in the coordinate system used by the map. If the coordinates do not correspond to the exact pixel center, the value of the nearest pixel center will be returned. For positions outside the map geometry `np.nan` is returned. The coordinate or idx arrays follow normal `Numpy broadcasting rules `__. So the following works as expected: .. GENERATED FROM PYTHON SOURCE LINES 659-664 .. code-block:: Python lons = np.linspace(-4, 4, 10) print(m_gc.get_by_coord({"lon": lons, "lat": 0})) .. rst-class:: sphx-glr-script-out .. code-block:: none [0. 0. 0. 0. 0. 0. 0. 0. 0. 0.] .. GENERATED FROM PYTHON SOURCE LINES 665-668 Or as an even more advanced example, we can provide `lats` as column vector and broadcasting to a 2D result array will be applied: .. GENERATED FROM PYTHON SOURCE LINES 668-674 .. code-block:: Python lons = np.linspace(-4, 4, 8) lats = np.linspace(-4, 4, 8).reshape(-1, 1) print(m_gc.get_by_coord({"lon": lons, "lat": lats})) .. rst-class:: sphx-glr-script-out .. code-block:: none [[nan nan nan nan nan nan nan nan] [nan nan nan nan nan nan nan nan] [ 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0.] [ 0. 0. 0. 0. 0. 0. 0. 0.] [nan nan nan nan nan nan nan nan] [nan nan nan nan nan nan nan nan]] .. GENERATED FROM PYTHON SOURCE LINES 675-687 Indexing and Slicing Sub-Maps ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ When you have worked with Numpy arrays in the past you are probably familiar with the concept of indexing and slicing into data arrays. To support slicing of non-spatial axes of `Map` objects, the `Map` object has a `~gammapy.maps.Map.slice_by_idx()` method, which allows to extract sub-maps from a larger map. The following example demonstrates how to get the map at the energy bin number 3: .. GENERATED FROM PYTHON SOURCE LINES 687-692 .. code-block:: Python m_sub = m_cube.slice_by_idx({"energy": 3}) print(m_sub) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 unit : dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 693-700 Note that the returned object is again a `~gammapy.maps.Map` with updated axes information. In this case, because we extracted only a single image, the energy axes is dropped from the map. To extract a sub-cube with a sliced energy axes we can use a normal ``slice()`` object: .. GENERATED FROM PYTHON SOURCE LINES 700-705 .. code-block:: Python m_sub = m_cube.slice_by_idx({"energy": slice(1, 3)}) print(m_sub) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy'] shape : (500, 250, 2) ndim : 3 unit : dtype : >f4 .. GENERATED FROM PYTHON SOURCE LINES 706-712 Note that the returned object is also a `~gammapy.maps.Map` object, but this time with updated energy axis specification. Slicing of multiple dimensions is supported by adding further entries to the dict passed to `~gammapy.maps.Map.slice_by_idx()` .. GENERATED FROM PYTHON SOURCE LINES 712-717 .. code-block:: Python m_sub = m_4d.slice_by_idx({"energy": slice(1, 3), "time": slice(4, 10)}) print(m_sub) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy', 'time'] shape : (500, 250, 2, 6) ndim : 4 unit : dtype : float32 .. GENERATED FROM PYTHON SOURCE LINES 718-722 For convenience there is also a `~gammapy.maps.Map.get_image_by_coord()` method which allows to access image planes at given non-spatial physical coordinates. This method also supports `~astropy.units.Quantity` objects: .. GENERATED FROM PYTHON SOURCE LINES 722-727 .. code-block:: Python image = m_4d.get_image_by_coord({"energy": 4 * u.TeV, "time": 5 * u.h}) print(image.geom) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsGeom axes : ['lon', 'lat'] shape : (500, 250) ndim : 2 frame : galactic projection : CAR center : 0.0 deg, 0.0 deg width : 10.0 deg x 5.0 deg wcs ref : 0.0 deg, 0.0 deg .. GENERATED FROM PYTHON SOURCE LINES 728-737 Iterating by image ~~~~~~~~~~~~~~~~~~ For maps with non-spatial dimensions the `~gammapy.maps.Map.iter_by_image_data` method can be used to loop over image slices. The image plane index `idx` is returned in data order, so that the data array can be indexed directly. Here is an example for an in-place convolution of an image using `~astropy.convolution.convolve` to interpolate NaN values: .. GENERATED FROM PYTHON SOURCE LINES 737-750 .. code-block:: Python axis1 = MapAxis([1, 10, 100], interp="log", name="energy") axis2 = MapAxis([1, 2, 3], interp="lin", name="time") m = Map.create(width=(5, 3), axes=[axis1, axis2], binsz=0.1) m.data[:, :, 15:18, 20:25] = np.nan for img, idx in m.iter_by_image_data(): kernel = np.ones((5, 5)) m.data[idx] = convolve(img, kernel) assert not np.isnan(m.data).any() .. GENERATED FROM PYTHON SOURCE LINES 751-760 Modifying Data -------------- How to set data values ~~~~~~~~~~~~~~~~~~~~~~ To modify and set map data values the `Map` object features as well a `~gammapy.maps.Map.set_by_idx()` method: .. GENERATED FROM PYTHON SOURCE LINES 760-764 .. code-block:: Python m_cube.set_by_idx(idx=(10, 20, 3), vals=42) .. GENERATED FROM PYTHON SOURCE LINES 765-767 here we check that data have been updated: .. GENERATED FROM PYTHON SOURCE LINES 767-771 .. code-block:: Python print(m_cube.get_by_idx((10, 20, 3))) .. rst-class:: sphx-glr-script-out .. code-block:: none [42.] .. GENERATED FROM PYTHON SOURCE LINES 772-775 Of course there is also a `~gammapy.maps.Map.set_by_coord()` method, which allows to set map data values in physical coordinates. .. GENERATED FROM PYTHON SOURCE LINES 775-779 .. code-block:: Python m_cube.set_by_coord({"lon": 0, "lat": 0, "energy": 2 * u.TeV}, vals=42) .. GENERATED FROM PYTHON SOURCE LINES 780-789 Again the `lon` and `lat` values are assumed to be given in degrees in the coordinate system used by the map. For the energy axis, the unit is the one specified on the axis (use ``m_cube.geom.axes[0].unit`` to check if needed…). All ``.xxx_by_coord()`` methods accept `~astropy.coordinates.SkyCoord` objects as well. In this case we have to use the ``"skycoord"`` keyword instead of ``"lon"`` and ``"lat"``: .. GENERATED FROM PYTHON SOURCE LINES 789-794 .. code-block:: Python skycoords = SkyCoord([1.2, 3.4], [-0.5, 1.1], frame="galactic", unit="deg") m_cube.set_by_coord({"skycoord": skycoords, "energy": 2 * u.TeV}, vals=42) .. GENERATED FROM PYTHON SOURCE LINES 795-800 Filling maps from event lists ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ This example shows how to fill a counts cube from an event list: .. GENERATED FROM PYTHON SOURCE LINES 800-814 .. code-block:: Python energy_axis = MapAxis.from_bounds( 10.0, 2e3, 12, interp="log", name="energy", unit="GeV" ) counts_3d = WcsNDMap.create( binsz=0.1, width=10.0, skydir=(0, 0), frame="galactic", axes=[energy_axis] ) events = EventList.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-events.fits.gz") counts_3d.fill_by_coord({"skycoord": events.radec, "energy": events.energy}) counts_3d.write("ccube.fits", format="fgst-ccube", overwrite=True) .. GENERATED FROM PYTHON SOURCE LINES 815-817 Alternatively you can use the `~gammapy.maps.Map.fill_events` method: .. GENERATED FROM PYTHON SOURCE LINES 817-825 .. code-block:: Python counts_3d = WcsNDMap.create( binsz=0.1, width=10.0, skydir=(0, 0), frame="galactic", axes=[energy_axis] ) counts_3d.fill_events(events) .. GENERATED FROM PYTHON SOURCE LINES 826-830 If you have a given map already, and want to make a counts image with the same geometry (not using the pixel data from the original map), you can also use the `~gammapy.maps.Map.fill_events` method. .. GENERATED FROM PYTHON SOURCE LINES 830-837 .. code-block:: Python events = EventList.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-events.fits.gz") reference_map = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz") counts = Map.from_geom(reference_map.geom) counts.fill_events(events) .. GENERATED FROM PYTHON SOURCE LINES 838-842 It works for IACT and Fermi-LAT events, for WCS or HEALPix map geometries, and also for extra axes. Especially energy axes are automatically handled correctly. .. GENERATED FROM PYTHON SOURCE LINES 845-866 Filling maps from interpolation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Maps support interpolation via the `~~gammapy.maps.Map.interp_by_coord` and `~~gammapy.maps.Map.interp_by_pix` methods. Currently, the following interpolation methods are supported: - ``"nearest"`` : Return value of nearest pixel (no interpolation). - ``"linear"`` : Interpolation with first order polynomial. This is the only interpolation method that is supported for all map types. - `quadratic` : Interpolation with second order polynomial. - `cubic` : Interpolation with third order polynomial. Note that ``"quadratic"`` and ``"cubic"`` interpolation are currently only supported for WCS-based maps with regular geometry (e.g. 2D or ND with the same geometry in every image plane). ``"linear"`` and higher order interpolation by pixel coordinates is only supported for WCS-based maps. In the following example we create a new map and fill it by interpolating another map: .. GENERATED FROM PYTHON SOURCE LINES 866-887 .. code-block:: Python # read map filename = "$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz" m_iem_gc = Map.read(filename) # create new geometry skydir = SkyCoord(266.4, -28.9, frame="icrs", unit="deg") wcs_geom_cel = WcsGeom.create(skydir=skydir, binsz=0.1, frame="icrs", width=(8, 4)) # create new empty map from geometry m_iem_10GeV = Map.from_geom(wcs_geom_cel) coords = m_iem_10GeV.geom.get_coord() # fill new map using interpolation m_iem_10GeV.data = m_iem_gc.interp_by_coord( {"skycoord": coords.skycoord, "energy_true": 10 * u.GeV}, method="linear", fill_value=np.nan, ) .. GENERATED FROM PYTHON SOURCE LINES 888-891 Interpolating onto a different geometry ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 894-898 For 3d geometries this operation can be performed directly using the `~gammapy.maps.Map.interp_to_geom()` method. This is very useful, ex: while using map arithmetic. .. GENERATED FROM PYTHON SOURCE LINES 898-915 .. code-block:: Python # create new geometry energy_axis = MapAxis.from_bounds( 10.0, 2e3, 6, interp="log", name="energy_true", unit="GeV" ) skydir = SkyCoord(266.4, -28.9, frame="icrs", unit="deg") wcs_geom_3d = WcsGeom.create( skydir=skydir, binsz=0.1, frame="icrs", width=(8, 4), axes=[energy_axis] ) # create the interpolated map m_iem_interp = m_iem_gc.interp_to_geom( wcs_geom_3d, preserve_counts=False, method="linear", fill_value=np.nan ) print(m_iem_interp) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat', 'energy_true'] shape : (80, 40, 6) ndim : 3 unit : 1 / (cm2 MeV s sr) dtype : float64 .. GENERATED FROM PYTHON SOURCE LINES 916-920 Note that ``preserve_counts=`` option should be true if the map is an integral quantity (e.g. counts) and false if the map is a differential quantity (e.g. intensity). .. GENERATED FROM PYTHON SOURCE LINES 923-932 Maps operations --------------- Basic operators ~~~~~~~~~~~~~~~ One can perform simple arithmetic on maps using the `+`, `-`, `*`, `/` operators, this works only for maps with the same geometry: .. GENERATED FROM PYTHON SOURCE LINES 932-938 .. code-block:: Python iem_plus_iem = m_iem_10GeV + m_iem_10GeV iem_minus_iem = m_iem_10GeV - m_iem_10GeV .. GENERATED FROM PYTHON SOURCE LINES 939-942 These operations can be applied between a Map and a scalar in that specific order: .. GENERATED FROM PYTHON SOURCE LINES 942-947 .. code-block:: Python iem_times_two = m_iem_10GeV * 2 # iem_times_two = 2 * m_iem_10GeV # this won't work .. GENERATED FROM PYTHON SOURCE LINES 948-951 The logic operators can also be applied on maps (the result is a map of boolean type): .. GENERATED FROM PYTHON SOURCE LINES 951-956 .. code-block:: Python is_null = iem_minus_iem == 0 print(is_null) .. rst-class:: sphx-glr-script-out .. code-block:: none WcsNDMap geom : WcsGeom axes : ['lon', 'lat'] shape : (80, 40) ndim : 2 unit : dtype : bool .. GENERATED FROM PYTHON SOURCE LINES 957-960 Here we check that the result is `True` for all the well-defined pixels (not `NaN`): .. GENERATED FROM PYTHON SOURCE LINES 960-964 .. code-block:: Python print(np.all(is_null.data[~np.isnan(iem_minus_iem)])) .. rst-class:: sphx-glr-script-out .. code-block:: none True .. GENERATED FROM PYTHON SOURCE LINES 965-972 Cutouts ~~~~~~~ The `WCSNDMap` objects features a `~gammapy.maps.Map.cutout()` method, which allows you to cut out a smaller part of a larger map. This can be useful, e.g. when working with all-sky diffuse maps. Here is an example: .. GENERATED FROM PYTHON SOURCE LINES 972-977 .. code-block:: Python position = SkyCoord(0, 0, frame="galactic", unit="deg") m_iem_cutout = m_iem_gc.cutout(position=position, width=(4 * u.deg, 2 * u.deg)) .. GENERATED FROM PYTHON SOURCE LINES 978-984 The returned object is again a `~gammapy.maps.Map` object with updated WCS information and data size. As one can see the cutout is automatically applied to all the non-spatial axes as well. The cutout width is given in the order of `(lon, lat)` and can be specified with units that will be handled correctly. .. GENERATED FROM PYTHON SOURCE LINES 987-1002 Visualizing and Plotting ------------------------ All map objects provide a `~gammapy.maps.Map.plot` method for generating a visualization of a map. This method returns figure, axes, and image objects that can be used to further tweak/customize the image. The `~gammapy.maps.Map.plot` method should be used with 2D maps, while 3D maps can be displayed with the `~gammapy.maps.Map.plot_interactive()` or `~gammapy.maps.Map.plot_grid()` methods. Image Plotting ~~~~~~~~~~~~~~ For debugging and inspecting the map data it is useful to plot or visualize the images planes contained in the map. .. GENERATED FROM PYTHON SOURCE LINES 1002-1007 .. code-block:: Python filename = "$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz" m_3fhl_gc = Map.read(filename) .. GENERATED FROM PYTHON SOURCE LINES 1008-1011 After reading the map we can now plot it on the screen by calling the ``.plot()`` method: .. GENERATED FROM PYTHON SOURCE LINES 1011-1016 .. code-block:: Python m_3fhl_gc.plot() plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_maps_001.png :alt: maps :srcset: /tutorials/api/images/sphx_glr_maps_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1017-1022 We can easily improve the plot by calling `~gammapy.maps.Map.smooth()` first and providing additional arguments to `~gammapy.maps.Map.plot()`. Most of them are passed further to `plt.imshow() `__: .. GENERATED FROM PYTHON SOURCE LINES 1022-1028 .. code-block:: Python smoothed = m_3fhl_gc.smooth(width=0.2 * u.deg, kernel="gauss") smoothed.plot(stretch="sqrt", add_cbar=True, vmax=4, cmap="inferno") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_maps_002.png :alt: maps :srcset: /tutorials/api/images/sphx_glr_maps_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1029-1034 We can use the `plt.rc_context() `__ context manager to further tweak the plot by adapting the figure and font size: .. GENERATED FROM PYTHON SOURCE LINES 1034-1042 .. code-block:: Python rc_params = {"figure.figsize": (12, 5.4), "font.size": 12} with plt.rc_context(rc=rc_params): smoothed = m_3fhl_gc.smooth(width=0.2 * u.deg, kernel="gauss") smoothed.plot(stretch="sqrt", add_cbar=True, vmax=4) plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_maps_003.png :alt: maps :srcset: /tutorials/api/images/sphx_glr_maps_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 1043-1052 Cube plotting ~~~~~~~~~~~~~ For maps with non-spatial dimensions the `~gammapy.maps.Map` object features an interactive plotting method, that works in jupyter notebooks only (Note: it requires the package `ipywidgets` to be installed). We first read a small example cutout from the Fermi Galactic diffuse model and display the data cube by calling `~gammapy.maps.Map.plot_interactive()`: .. GENERATED FROM PYTHON SOURCE LINES 1052-1062 .. code-block:: Python rc_params = { "figure.figsize": (12, 5.4), "font.size": 12, "axes.formatter.limits": (2, -2), } m_iem_gc.plot_interactive(add_cbar=True, stretch="sqrt", rc_params=rc_params) plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_maps_004.png :alt: maps :srcset: /tutorials/api/images/sphx_glr_maps_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none interactive(children=(SelectionSlider(continuous_update=False, description='Select energy_true:', layout=Layout(width='50%'), options=('58.5 MeV', '80.0 MeV', '109 MeV', '150 MeV', '205 MeV', '280 MeV', '383 MeV', '523 MeV', '716 MeV', '979 MeV', '1.34 GeV', '1.83 GeV', '2.50 GeV', '3.42 GeV', '4.68 GeV', '6.41 GeV', '8.76 GeV', '12.0 GeV', '16.4 GeV', '22.4 GeV', '30.6 GeV', '41.9 GeV', '57.3 GeV', '78.4 GeV', '107 GeV', '147 GeV', '201 GeV', '274 GeV', '375 GeV', '513 GeV'), style=SliderStyle(description_width='initial'), value='58.5 MeV'), RadioButtons(description='Select stretch:', index=1, options=('linear', 'sqrt', 'log'), style=DescriptionStyle(description_width='initial'), value='sqrt'), Output()), _dom_classes=('widget-interact',)) .. GENERATED FROM PYTHON SOURCE LINES 1063-1076 Now you can use the interactive slider to select an energy range and the corresponding image is displayed on the screen. You can also use the radio buttons to select your preferred image stretching. We have passed additional keywords using the `rc_params` argument to improve the figure and font size. Those keywords are directly passed to the `plt.rc_context() `__ context manager. Additionally all the slices of a 3D `~gammapy.maps.Map` can be displayed using the `~gammapy.maps.Map.plot_grid()` method. By default the colorbars bounds of the subplots are not the same, we can make them consistent using the `vmin` and `vmax` options: .. GENERATED FROM PYTHON SOURCE LINES 1076-1079 .. code-block:: Python counts_3d.plot_grid(ncols=4, figsize=(16, 12), vmin=0, vmax=100, stretch="log") plt.show() .. image-sg:: /tutorials/api/images/sphx_glr_maps_005.png :alt: Energy 10.0 GeV - 15.6 GeV, Energy 15.6 GeV - 24.2 GeV, Energy 24.2 GeV - 37.6 GeV, Energy 37.6 GeV - 58.5 GeV, Energy 58.5 GeV - 90.9 GeV, Energy 90.9 GeV - 141 GeV, Energy 141 GeV - 220 GeV, Energy 220 GeV - 342 GeV, Energy 342 GeV - 532 GeV, Energy 532 GeV - 827 GeV, Energy 827 GeV - 1.29 TeV, Energy 1.29 TeV - 2.00 TeV :srcset: /tutorials/api/images/sphx_glr_maps_005.png :class: sphx-glr-single-img .. _sphx_glr_download_tutorials_api_maps.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/api/maps.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: maps.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: maps.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_