maps - Sky maps

Introduction

gammapy.maps contains classes for representing pixelized data structures with at least two spatial dimensions representing coordinates on a sphere (e.g. an image in celestial coordinates). These classes support an arbitrary number of non-spatial dimensions and can represent images (2D), cubes (3D), or hypercubes (4+D). Two pixelization schemes are supported:

• WCS : Projection onto a 2D cartesian grid following the conventions of the World Coordinate System (WCS). Pixels are square in projected coordinates and as such are not equal area in spherical coordinates.

• HEALPix : Hierarchical Equal Area Iso Latitude pixelation of the sphere. Pixels are equal area but have irregular shapes.

gammapy.maps is organized around two data structures: geometry classes inheriting from Geom and map classes inheriting from Map. A geometry defines the map boundaries, pixelization scheme, and provides methods for converting to/from map and pixel coordinates. A map owns a Geom instance as well as a data array containing map values. Where possible it is recommended to use the abstract Map interface for accessing or updating the contents of a map as this allows algorithms to be used interchangeably with different map representations. The following reviews methods of the abstract map interface. Documentation specific to WCS- and HEALPix-based maps is provided in HEALPix-based Maps.

Getting Started

All map objects have an abstract inteface provided through the methods of the Map. These methods can be used for accessing and manipulating the contents of a map without reference to the underlying data representation (e.g. whether a map uses WCS or HEALPix pixelization). For applications which do depend on the specific representation one can also work directly with the classes derived from Map. In the following we review some of the basic methods for working with map objects.

Constructing with Factory Methods

The Map class provides a create factory method to facilitate creating an empty map object from scratch. The map_type argument can be used to control the pixelization scheme (WCS or HPX).

from gammapy.maps import Map
from astropy.coordinates import SkyCoord

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)

Higher dimensional map objects (cubes and hypercubes) can be constructed by passing a list of MapAxis objects for non-spatial dimensions with the axes parameter:

from gammapy.maps import Map, MapAxis
from astropy.coordinates import SkyCoord

position = SkyCoord(0.0, 5.0, frame='galactic', unit='deg')
energy_axis = MapAxis.from_bounds(100., 1E5, 12, interp='log', name='energy', unit='GeV')

# Create a WCS Map
m_wcs = Map.create(binsz=0.1, map_type='wcs', skydir=position, width=10.0,
                      axes=[energy_axis])

# Create a HPX Map
m_hpx = Map.create(binsz=0.1, map_type='hpx', skydir=position, width=10.0,
                      axes=[energy_axis])

Multi-resolution maps (maps with a different pixel size or geometry in each image plane) can be constructed by passing a vector argument for any of the geometry parameters. This vector must have the same shape as the non-spatial dimensions of the map. The following example demonstrates creating an energy cube with a pixel size proportional to the Fermi-LAT PSF:

import numpy as np
from gammapy.maps import Map, MapAxis
from astropy.coordinates import SkyCoord

position = SkyCoord(0.0, 5.0, frame='galactic', unit='deg')
energy_axis = MapAxis.from_bounds(100., 1E5, 12, interp='log', name='energy', unit='GeV')

binsz = np.sqrt((3.0*(energy_axis.center/100.)**-0.8)**2 + 0.1**2)

# Create a WCS Map
m_wcs = Map.create(binsz=binsz, map_type='wcs', skydir=position, width=10.0,
                      axes=[energy_axis])

# Create a HPX Map
m_hpx = Map.create(binsz=binsz, map_type='hpx', skydir=position, width=10.0,
                      axes=[energy_axis])

Indexing and Slicing

All map objects feature a slice_by_idx() method, which can be used to slice and index non-spatial axes of the map to create arbitrary sub-maps. The method accepts a dict specifying the axes name and correspoding integer index or slice objects. When indexing an axis with an integer the corresponding axes is dropped from the returned sub-map. To keep the axes (with length 1) in the returned sub-map use a slice object of length one. This behaviour is equivalent to regular numpy array indexing. The following example demonstrates the use of slice_by_idx() on a map with a time and energy axes:

import numpy as np
from gammapy.maps import Map, MapAxis
from astropy.coordinates import SkyCoord

position = SkyCoord(0.0, 5.0, frame='galactic', unit='deg')
energy_axis = MapAxis.from_bounds(100., 1E5, 12, interp='log', unit='GeV', name='energy')
time_axis = MapAxis.from_bounds(0., 12, 12, interp='lin', unit='h', name='time')

# Create a WCS Map
m_wcs = Map.create(binsz=0.02, map_type='wcs', skydir=position, width=10.0,
                      axes=[energy_axis, time_axis])

# index first image plane of the energy axes and third from the time axis
m_wcs.slice_by_idx({'energy': 0, 'time': 2})

# index first image plane of the energy axes and keep time axis unchanged
m_wcs.slice_by_idx({'energy': 0})

# slice first three images of the energy axis at a fixed time
m_wcs.slice_by_idx({'energy': slice(0, 3), 'time': 0})

# slice first three images of the energy axis as well as time axis
m_wcs.slice_by_idx({'energy': slice(0, 3), 'time': slice(0, 3)})

Accessor Methods

All map objects have a set of accessor methods provided through the abstract Map class. These methods can be used to access or update the contents of the map irrespective of its underlying representation. Four types of accessor methods are provided:

• get : Return the value of the map at the pixel containing the given coordinate (get_by_idx, get_by_pix, get_by_coord).

• interp : Interpolate or extrapolate the value of the map at an arbitrary coordinate (see also Interpolation).

• set : Set the value of the map at the pixel containing the given coordinate (set_by_idx, set_by_pix, set_by_coord).

• fill : Increment the value of the map at the pixel containing the given coordinate with a unit weight or the value in the optional weights argument (fill_by_idx, fill_by_pix, fill_by_coord).

Accessor methods accept as their first argument a coordinate tuple containing scalars, lists, or numpy arrays with one tuple element for each dimension of the map. coord methods optionally support a dict or MapCoord argument.

When using tuple input the first two elements in the tuple should be longitude and latitude followed by one element for each non-spatial dimension. Map coordinates can be expressed in one of three coordinate systems:

• idx : Pixel indices. These are explicit (integer) pixel indices into the map.

• pix : Coordinates in pixel space. Pixel coordinates are continuous defined on the interval [0,N-1] where N is the number of pixels along a given map dimension with pixel centers at integer values. For methods that reference a discrete pixel, pixel coordinates wil be rounded to the nearest pixel index and passed to the corresponding idx method.

• coord : The true map coordinates including angles on the sky (longitude and latitude). This coordinate system supports three coordinate representations: tuple, dict, and MapCoord. The tuple representation should contain longitude and latitude in degrees followed by one coordinate array for each non-spatial dimension.

The coordinate system accepted by a given accessor method can be inferred from the suffix of the method name (e.g. get_by_idx). The following demonstrates how one can access the same pixels of a WCS map using each of the three coordinate systems:

from gammapy.maps import Map

m = Map.create(binsz=0.1, map_type='wcs', width=10.0)

vals = m.get_by_idx( ([49,50],[49,50]) )
vals = m.get_by_pix( ([49.0,50.0],[49.0,50.0]) )
vals = m.get_by_coord( ([-0.05,-0.05],[0.05,0.05]) )

Coordinate arguments obey normal numpy broadcasting rules. The coordinate tuple may contain any combination of scalars, lists or numpy arrays as long as they have compatible shapes. For instance a combination of scalar and vector arguments can be used to perform an operation along a slice of the map at a fixed value along that dimension. Multi-dimensional arguments can be use to broadcast a given operation across a grid of coordinate values.

from gammapy.maps import Map

m = Map.create(binsz=0.1, map_type='wcs', width=10.0)
coords = np.linspace(-4.0, 4.0, 9)

# Equivalent calls for accessing value at pixel (49,49)
vals = m.get_by_idx( (49,49) )
vals = m.get_by_idx( ([49],[49]) )
vals = m.get_by_idx( (np.array([49]), np.array([49])) )

# Retrieve map values along latitude at fixed longitude=0.0
vals = m.get_by_coord( (0.0, coords) )
# Retrieve map values on a 2D grid of latitude/longitude points
vals = m.get_by_coord( (coords[None,:], coords[:,None]) )
# Set map values along slice at longitude=0.0 to twice their existing value
m.set_by_coord((0.0, coords), 2.0*m.get_by_coord((0.0, coords)))

The set and fill methods can both be used to set pixel values. The following demonstrates how one can set pixel values:

from gammapy.maps import Map

m = Map.create(binsz=0.1, map_type='wcs', width=10.0)

m.set_by_coord(([-0.05, -0.05], [0.05, 0.05]), [0.5, 1.5])
m.fill_by_coord( ([-0.05, -0.05], [0.05, 0.05]), weights=[0.5, 1.5])

Interface with MapCoord and SkyCoord

The coord accessor methods accept dict, MapCoord, and SkyCoord arguments in addition to the standard tuple of ndarray argument. When using a tuple argument a SkyCoord can be used instead of longitude and latitude arrays. The coordinate frame of the SkyCoord will be transformed to match the coordinate system of the map.

import numpy as np
from astropy.coordinates import SkyCoord
from gammapy.maps import Map, MapCoord, MapAxis

lon = [0, 1]
lat = [1, 2]
energy = [100, 1000]
energy_axis = MapAxis.from_bounds(100, 1E5, 12, interp='log', name='energy')

skycoord = SkyCoord(lon, lat, unit='deg', frame='galactic')
m = Map.create(binsz=0.1, map_type='wcs', width=10.0,
              frame="galactic", axes=[energy_axis])

m.set_by_coord((skycoord, energy), [0.5, 1.5])
m.get_by_coord((skycoord, energy))

A MapCoord or dict argument can be used to interact with a map object without reference to the axis ordering of the map geometry:

coord = MapCoord.create(dict(lon=lon, lat=lat, energy=energy))
m.set_by_coord(coord, [0.5, 1.5])
m.get_by_coord(coord)
m.set_by_coord(dict(lon=lon, lat=lat, energy=energy), [0.5, 1.5])
m.get_by_coord(dict(lon=lon, lat=lat, energy=energy))

However when using the named axis interface the axis name string (e.g. as given by MapAxis.name) must match the name given in the method argument. The two spatial axes must always be named lon and lat.

Differential and integral maps

gammapy.maps supports both differential and integral maps, representing differential values at specific coordinates, or integral values within bins. This is achieved by specifying the node_type of a MapAxis. Quantities defined at bin centers should have a node_type of “center”, and quantities integrated in bins should have node_type of “edges”. Interpolation is defined only for differential quantities.

For the specific case of the energy axis, conventionally, true energies are have node_type “center” (usually used for IRFs and exposure) whereas the reconstructed energy axis has node_type “edges” (usually used for counts and background). Model evaluations are first computed on differential bins, and then multiplied by the bin volumes to finally return integrated maps, so the output predicted counts maps are integral with node_type “edges”.

MapCoord

MapCoord is an N-dimensional coordinate object that stores both spatial and non-spatial coordinates and is accepted by all coord methods. A MapCoord can be created with or without explicitly named axes with MapCoord.create. Axes of a MapCoord can be accessed by index, name, or attribute. A MapCoord without explicit axis names can be created by calling MapCoord.create with a tuple argument:

import numpy as np
from astropy.coordinates import SkyCoord
from gammapy.maps import MapCoord

lon = [0.0, 1.0]
lat = [1.0, 2.0]
energy = [100, 1000]
skycoord = SkyCoord(lon, lat, unit='deg', frame='galactic')

# Create a MapCoord from a tuple (no explicit axis names)
c = MapCoord.create((lon, lat, energy))
print(c[0], c['lon'], c.lon)
print(c[1], c['lat'], c.lat)
print(c[2], c['axis0'])

# Create a MapCoord from a tuple + SkyCoord (no explicit axis names)
c = MapCoord.create((skycoord, energy))
print(c[0], c['lon'], c.lon)
print(c[1], c['lat'], c.lat)
print(c[2], c['axis0'])

The first two elements of the tuple argument must contain longitude and latitude. Non-spatial axes are assigned a default name axis{I} where {I} is the index of the non-spatial dimension. MapCoord objects created without named axes must have the same axis ordering as the map geometry.

A MapCoord with named axes can be created by calling MapCoord.create with a dict:

c = MapCoord.create(dict(lon=lon, lat=lat, energy=energy))
print(c[0], c['lon'], c.lon)
print(c[1], c['lat'], c.lat)
print(c[2], c['energy'])

c = MapCoord.create({'energy': energy, 'lon': lon, 'lat': lat})
print(c[0], c['energy'])
print(c[1], c['lon'], c.lon)
print(c[2], c['lat'], c.lat)

c = MapCoord.create(dict(skycoord=skycoord, energy=energy))
print(c[0], c['lon'], c.lon)
print(c[1], c['lat'], c.lat)
print(c[2], c['energy'])

Spatial axes must be named lon and lat. MapCoord objects created with named axes do not need to have the same axis ordering as the map geometry. However the name of the axis must match the name of the corresponding map geometry axis.

Interpolation

Maps support interpolation via the interp_by_coord and 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.

from gammapy.maps import Map

m = Map.create(binsz=0.1, map_type='wcs', width=10.0)

m.interp_by_coord(([-0.05, -0.05], [0.05, 0.05]), interp='linear')
m.interp_by_coord(([-0.05, -0.05], [0.05, 0.05]), interp='cubic')

Iterating by image

For maps with non-spatial dimensions the iter_by_image 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:

import numpy as np
from astropy.convolution import convolve

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():
    kernel = np.ones((5, 5))
    m.data[idx] = convolve(img, kernel)

assert not np.isnan(m.data).any()

FITS I/O

Maps can be written to and read from a FITS file with the write and read methods:

from gammapy.maps import Map

m = Map.create(binsz=0.1, map_type='wcs', width=10.0)
m.write('file.fits', hdu='IMAGE')
m = Map.read('file.fits', hdu='IMAGE')

If map_type argument is not given when calling read a map object will be instantiated with the pixelization of the input HDU.

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.

from gammapy.maps import Map

m = Map.create(binsz=0.1, map_type='wcs', width=10.0)
m.write('file.fits', hdu='IMAGE', sparse=True)
m = Map.read('file.fits', hdu='IMAGE', map_type='wcs')

By default files will be written to the gamma-astro-data-format specification for sky maps (see here). The GADF format offers a number of enhancements over existing map formats such as support for writing multi-resolution maps, sparse maps, and cubes with different geometries to the same file. For backward compatibility with software using other formats, the conv keyword option is provided to write a file using a format other than the GADF format:

from gammapy.maps import Map, MapAxis

energy_axis = MapAxis.from_bounds(100., 1E5, 12, interp='log')
m = Map.create(binsz=0.1, map_type='wcs', width=10.0,
                  axes=[energy_axis])
# Write a counts cube in a format compatible with the Fermi Science Tools
m.write('ccube.fits', conv='fgst-ccube')

Visualization

All map objects provide a 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.

import matplotlib.pyplot as plt
from gammapy.maps import Map

m = Map.read("$GAMMAPY_DATA/fermi-3fhl-gc/fermi-3fhl-gc-counts.fits.gz")
m.smooth("0.1 deg").plot(cmap="inferno", add_cbar=True, stretch="sqrt")
plt.show()

Examples

Creating counts cubes from event lists

This example shows how to fill a counts cube from an event list:

from gammapy.data import EventList
from gammapy.maps import WcsGeom, WcsNDMap, MapAxis

energy_axis = MapAxis.from_bounds(10., 2E3, 12, interp='log', name='energy', unit='GeV')
m = 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")

m.fill_by_coord({'skycoord': events.radec, 'energy': events.energy})
m.write('ccube.fits', conv='fgst-ccube')

To make a counts map, create an empty map with a geometry of your choice and then fill it using this function

from gammapy.maps import Map
from gammapy.data import EventList
events = EventList.read('$GAMMAPY_DATA/cta-1dc/data/baseline/gps/gps_baseline_110380.fits')
counts = Map.create(frame="galactic", skydir=(0, 0), binsz=0.1, npix=(120, 100))
counts.fill_events(events)
counts.plot()

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), do this

from gammapy.maps import Map
from gammapy.data import EventList
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)
counts.smooth(3).plot()

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.

Generating a Cutout of a Model Cube

This example shows how to extract a cut-out of LAT galactic diffuse model cube using the WcsNDMap.cutout method:

from gammapy.maps import WcsGeom, WcsNDMap
from astropy import units as u
from astropy.coordinates import SkyCoord

m = WcsNDMap.read('$GAMMAPY_DATA/fermi-3fhl-gc/gll_iem_v06_gc.fits.gz')
position = SkyCoord(0, 0, frame="galactic", unit="deg")
m_cutout = m.cutout(position=position, width=(5 * u.deg, 2 * u.deg))
m_cutout.write('cutout.fits', conv='fgst-template')

Using gammapy.maps

Tutorials that show examples using gammapy.maps:

• maps.html

• analysis_3d.html

• simulate_3d.html

• fermi_lat.html

More detailed documentation on the WCS and HPX classes in gammapy.maps can be found in the following sub-pages:

• HEALPix-based Maps
• Image plotting

Reference/API

gammapy.maps Package

Sky maps.

Functions

colormap_hess([transition, width])	Colormap often used in H.E.S.S.
colormap_milagro([transition, width, huestart])	Colormap often used in Milagro collaboration publications.
scale_cube(data, kernels)	Compute scale space cube.

Classes

Geom	Map geometry base class.
HpxGeom(nside[, nest, frame, region, axes, …])	Geometry class for HEALPIX maps.
HpxMap(geom, data[, meta, unit])	Base class for HEALPIX map classes.
HpxNDMap(geom[, data, dtype, meta, unit])	HEALPix map with any number of non-spatial dimensions.
ImageProfile(table)	Image profile class.
ImageProfileEstimator([x_edges, method, …])	Estimate profile from image.
Map(geom, data[, meta, unit])	Abstract map class.
MapAxis(nodes[, interp, name, node_type, unit])	Class representing an axis of a map.
MapCoord(data[, frame, match_by_name])	Represents a sequence of n-dimensional map coordinates.
MapPanelPlotter(figure, xlim, ylim[, npanels])	Map panel plotter class.
WcsGeom(wcs, npix[, cdelt, crpix, axes, …])	Geometry class for WCS maps.
WcsMap(geom, data[, meta, unit])	Base class for WCS map classes.
WcsNDMap(geom[, data, dtype, meta, unit])	HEALPix map with any number of non-spatial dimensions.