WcsGeom#

class gammapy.maps.WcsGeom(wcs, npix, cdelt=None, crpix=None, axes=None)[source]#

Bases: gammapy.maps.geom.Geom

Geometry class for WCS maps.

This class encapsulates both the WCS transformation object and the the image extent (number of pixels in each dimension). Provides methods for accessing the properties of the WCS object and performing transformations between pixel and world coordinates.

Parameters
wcsWCS

WCS projection object

npixtuple

Number of pixels in each spatial dimension

cdelttuple

Pixel size in each image plane. If none then a constant pixel size will be used.

crpixtuple

Reference pixel coordinate in each image plane.

axeslist

Axes for non-spatial dimensions

Attributes Summary

as_energy_true

If the geom contains an energy axis rename it to energy true

axes

List of non-spatial axes.

axes_names

All axes names

center_coord

Map coordinate of the center of the geometry.

center_pix

Pixel coordinate of the center of the geometry.

center_skydir

Sky coordinate of the center of the geometry.

data_shape

Shape of the Numpy data array matching this geometry.

data_shape_axes

Shape of data of the non-spatial axes and unit spatial axes.

footprint

Footprint of the geometry

frame

Coordinate system of the projection.

has_energy_axis

Whether geom has an energy axis

is_allsky

Flag for all-sky maps.

is_flat

Whether the geom non spatial axes have length 1, equivalent to an image.

is_hpx

is_image

Whether the geom is an image without extra dimensions.

is_region

is_regular

Is this geometry is regular in non-spatial dimensions (bool)?

ndim

npix

Tuple with image dimension in pixels in longitude and latitude.

pixel_area

Pixel area in deg^2.

pixel_scales

Pixel scale.

projection

Map projection.

shape_axes

Shape of non-spatial axes.

wcs

WCS projection object.

width

Tuple with image dimension in deg in longitude and latitude.

Methods Summary

bin_volume()

Bin volume (Quantity)

binary_structure(width[, kernel])

Get binary structure

boundary_mask(width)

Create a mask applying binary erosion with a given width from geom edges

contains(coords)

Check if a given map coordinate is contained in the geometry.

contains_pix(pix)

Check if a given pixel coordinate is contained in the geometry.

coord_to_idx(coords[, clip])

Convert map coordinates to pixel indices.

coord_to_pix(coords)

Convert map coordinates to pixel coordinates.

copy(**kwargs)

Copy and overwrite given attributes.

create([npix, binsz, proj, frame, refpix, ...])

Create a WCS geometry object.

crop(crop_width)

Crop the geometry at the edges.

cutout(position, width[, mode, odd_npix])

Create a cutout around a given position.

cutout_slices(geom[, mode])

Compute cutout slices.

data_nbytes([dtype])

Estimate memory usage in megabytes of the Numpy data array matching this geometry depending on the given type.

downsample(factor[, axis_name])

Downsample the spatial dimension of the geometry by a given factor.

drop(axis_name)

Drop an axis from the geom.

energy_mask([energy_min, energy_max, ...])

Create a mask for a given energy range.

from_aligned(geom, skydir, width)

Create an aligned geometry from an existing one

from_hdulist(hdulist[, hdu, hdu_bands])

Load a geometry object from a FITS HDUList.

from_header(header[, hdu_bands, format])

Create a WCS geometry object from a FITS header.

get_coord([idx, mode, frame, sparse, axis_name])

Get map coordinates from the geometry.

get_idx([idx, flat])

Get tuple of pixel indices for this geometry.

get_pix([idx, mode])

Get map pix coordinates from the geometry.

is_aligned(other[, tolerance])

Check if WCS and extra axes are aligned.

is_allclose(other[, rtol_axes, atol_axes, ...])

Compare two data IRFs for equivalency

pad(pad_width, axis_name)

Pad the geometry at the edges.

pix_to_coord(pix)

Convert pixel coordinates to map coordinates.

pix_to_idx(pix[, clip])

Convert pixel coordinates to pixel indices.

region_mask(regions[, inside])

Create a mask from a given list of regions

region_weights(regions[, oversampling_factor])

Compute regions weights

rename_axes(names, new_names)

Rename axes contained in the geometry

replace_axis(axis)

Replace axis with a new one.

resample_axis(axis)

Resample geom to a new axis binning.

separation(center)

Compute sky separation wrt a given center.

slice_by_idx(slices)

Create a new geometry by slicing the non-spatial axes.

solid_angle()

Solid angle array (Quantity in sr).

squash(axis_name)

Squash geom axis.

to_bands_hdu([hdu_bands, format])

to_binsz(binsz)

Change pixel size of the geometry.

to_cube(axes)

Append non-spatial axes to create a higher-dimensional geometry.

to_even_npix()

Create a new geom object with an even number of pixel and a maximum size.

to_header()

to_image()

Create 2D image geometry (drop non-spatial dimensions).

to_odd_npix([max_radius])

Create a new geom object with an odd number of pixel and a maximum size.

upsample(factor[, axis_name])

Upsample the spatial dimension of the geometry by a given factor.

Attributes Documentation

as_energy_true#

If the geom contains an energy axis rename it to energy true

axes#

List of non-spatial axes.

axes_names#

All axes names

center_coord#

Map coordinate of the center of the geometry.

Returns
coordtuple
center_pix#

Pixel coordinate of the center of the geometry.

Returns
pixtuple
center_skydir#

Sky coordinate of the center of the geometry.

Returns
pixSkyCoord
data_shape#

Shape of the Numpy data array matching this geometry.

data_shape_axes#

Shape of data of the non-spatial axes and unit spatial axes.

footprint#

Footprint of the geometry

frame#

Coordinate system of the projection.

Galactic (“galactic”) or Equatorial (“icrs”).

has_energy_axis#

Whether geom has an energy axis

is_allsky#

Flag for all-sky maps.

is_flat#

Whether the geom non spatial axes have length 1, equivalent to an image.

is_hpx = False#
is_image#

Whether the geom is an image without extra dimensions.

is_region = False#
is_regular#

Is this geometry is regular in non-spatial dimensions (bool)?

  • False for multi-resolution or irregular geometries.

  • True if all image planes have the same pixel geometry.

ndim#
npix#

Tuple with image dimension in pixels in longitude and latitude.

pixel_area#

Pixel area in deg^2.

pixel_scales#

Pixel scale.

Returns angles along each axis of the image at the CRPIX location once it is projected onto the plane of intermediate world coordinates.

Returns
angle: Angle
projection#

Map projection.

shape_axes#

Shape of non-spatial axes.

wcs#

WCS projection object.

width#

Tuple with image dimension in deg in longitude and latitude.

Methods Documentation

bin_volume()[source]#

Bin volume (Quantity)

binary_structure(width, kernel='disk')[source]#

Get binary structure

Parameters
widthQuantity, str or float

If a float is given it interpreted as width in pixels. If an (angular) quantity is given it converted to pixels using geom.wcs.wcs.cdelt. The width corresponds to radius in case of a disk kernel, and the side length in case of a box kernel.

kernel{‘disk’, ‘box’}

Kernel shape

Returns
structurendarray

Binary structure

boundary_mask(width)[source]#

Create a mask applying binary erosion with a given width from geom edges

Parameters
widthtuple of Quantity

Angular sizes of the margin in (lon, lat) in that specific order. If only one value is passed, the same margin is applied in (lon, lat).

Returns
mask_mapWcsNDMap of boolean type

Boundary mask

contains(coords)[source]#

Check if a given map coordinate is contained in the geometry.

Parameters
coordstuple or MapCoord

Tuple of map coordinates.

Returns
containmentndarray

Bool array.

contains_pix(pix)#

Check if a given pixel coordinate is contained in the geometry.

Parameters
pixtuple

Tuple of pixel coordinates.

Returns
containmentndarray

Bool array.

coord_to_idx(coords, clip=False)#

Convert map coordinates to pixel indices.

Parameters
coordstuple or MapCoord

Coordinate values in each dimension of the map. This can either be a tuple of numpy arrays or a MapCoord object. If passed as a tuple then the ordering should be (longitude, latitude, c_0, …, c_N) where c_i is the coordinate vector for axis i.

clipbool

Choose whether to clip indices to the valid range of the geometry. If false then indices for coordinates outside the geometry range will be set -1.

Returns
pixtuple

Tuple of pixel indices in image and band dimensions. Elements set to -1 correspond to coordinates outside the map.

coord_to_pix(coords)[source]#

Convert map coordinates to pixel coordinates.

Parameters
coordstuple

Coordinate values in each dimension of the map. This can either be a tuple of numpy arrays or a MapCoord object. If passed as a tuple then the ordering should be (longitude, latitude, c_0, …, c_N) where c_i is the coordinate vector for axis i.

Returns
pixtuple

Tuple of pixel coordinates in image and band dimensions.

copy(**kwargs)#

Copy and overwrite given attributes.

Parameters
**kwargsdict

Keyword arguments to overwrite in the map geometry constructor.

Returns
copyGeom

Copied map geometry.

classmethod create(npix=None, binsz=0.5, proj='CAR', frame='icrs', refpix=None, axes=None, skydir=None, width=None)[source]#

Create a WCS geometry object.

Pixelization of the map is set with binsz and one of either npix or width arguments. For maps with non-spatial dimensions a different pixelization can be used for each image plane by passing a list or array argument for any of the pixelization parameters. If both npix and width are None then an all-sky geometry will be created.

Parameters
npixint or tuple or list

Width of the map in pixels. A tuple will be interpreted as parameters for longitude and latitude axes. For maps with non-spatial dimensions, list input can be used to define a different map width in each image plane. This option supersedes width.

widthfloat or tuple or list or string

Width of the map in degrees. A tuple will be interpreted as parameters for longitude and latitude axes. For maps with non-spatial dimensions, list input can be used to define a different map width in each image plane.

binszfloat or tuple or list

Map pixel size in degrees. A tuple will be interpreted as parameters for longitude and latitude axes. For maps with non-spatial dimensions, list input can be used to define a different bin size in each image plane.

skydirtuple or SkyCoord

Sky position of map center. Can be either a SkyCoord object or a tuple of longitude and latitude in deg in the coordinate system of the map.

frame{“icrs”, “galactic”}, optional

Coordinate system, either Galactic (“galactic”) or Equatorial (“icrs”).

axeslist

List of non-spatial axes.

projstring, optional

Any valid WCS projection type. Default is ‘CAR’ (Plate-Carrée projection). See WCS supported projections # noqa: E501

refpixtuple

Reference pixel of the projection. If None this will be set to the center of the map.

Returns
geomWcsGeom

A WCS geometry object.

Examples

>>> from gammapy.maps import WcsGeom
>>> from gammapy.maps import MapAxis
>>> axis = MapAxis.from_bounds(0,1,2)
>>> geom = WcsGeom.create(npix=(100,100), binsz=0.1)
>>> geom = WcsGeom.create(npix=(100,100), binsz="0.1deg")
>>> geom = WcsGeom.create(npix=[100,200], binsz=[0.1,0.05], axes=[axis])
>>> geom = WcsGeom.create(npix=[100,200], binsz=["0.1deg","0.05deg"], axes=[axis])
>>> geom = WcsGeom.create(width=[5.0,8.0], binsz=[0.1,0.05], axes=[axis])
>>> geom = WcsGeom.create(npix=([100,200],[100,200]), binsz=0.1, axes=[axis])
crop(crop_width)[source]#

Crop the geometry at the edges.

Parameters
crop_width{sequence, array_like, int}

Number of values cropped from the edges of each axis.

Returns
geomGeom

Cropped geometry.

cutout(position, width, mode='trim', odd_npix=False)[source]#

Create a cutout around a given position.

Parameters
positionSkyCoord

Center position of the cutout region.

widthtuple of Angle

Angular sizes of the region in (lon, lat) in that specific order. If only one value is passed, a square region is extracted.

mode{‘trim’, ‘partial’, ‘strict’}

Mode option for Cutout2D, for details see Cutout2D.

odd_npixbool

Force width to odd number of pixels.

Returns
cutoutWcsNDMap

Cutout map

cutout_slices(geom, mode='partial')[source]#

Compute cutout slices.

Parameters
geomWcsGeom

Parent geometry

mode{“trim”, “partial”, “strict”}

Cutout slices mode.

Returns
slicesdict

Dictionary containing “parent-slices” and “cutout-slices”.

data_nbytes(dtype='float32')#

Estimate memory usage in megabytes of the Numpy data array matching this geometry depending on the given type.

Parameters
dtypedata-type

The desired data-type for the array. Default is “float32”

Returns
memoryQuantity

Estimated memory usage in megabytes (MB)

downsample(factor, axis_name=None)[source]#

Downsample the spatial dimension of the geometry by a given factor.

Parameters
factorint

Downsampling factor.

axis_namestr

Axis to downsample.

Returns
geomGeom

Downsampled geometry.

drop(axis_name)#

Drop an axis from the geom.

Parameters
axis_namestr

Name of the axis to remove.

Returns
geomGeom

New geom with the axis removed.

energy_mask(energy_min=None, energy_max=None, round_to_edge=False)#

Create a mask for a given energy range.

The energy bin must be fully contained to be included in the mask.

Parameters
energy_min, energy_maxQuantity

Energy range

Returns
maskndarray

Energy mask

classmethod from_aligned(geom, skydir, width)[source]#

Create an aligned geometry from an existing one

Parameters
geomWcsGeom

A reference WCS geometry object.

skydirtuple or SkyCoord

Sky position of map center. Can be either a SkyCoord object or a tuple of longitude and latitude in deg in the coordinate system of the map.

widthfloat or tuple or list or string

Width of the map in degrees. A tuple will be interpreted as parameters for longitude and latitude axes. For maps with non-spatial dimensions, list input can be used to define a different map width in each image plane.

Returns
geomWcsGeom

An aligned WCS geometry object with specified size and center.

classmethod from_hdulist(hdulist, hdu=None, hdu_bands=None)#

Load a geometry object from a FITS HDUList.

Parameters
hdulistHDUList

HDU list containing HDUs for map data and bands.

hdustr

Name or index of the HDU with the map data.

hdu_bandsstr

Name or index of the HDU with the BANDS table. If not defined this will be inferred from the FITS header of the map HDU.

Returns
geomGeom

Geometry object.

classmethod from_header(header, hdu_bands=None, format='gadf')[source]#

Create a WCS geometry object from a FITS header.

Parameters
headerHeader

The FITS header

hdu_bandsBinTableHDU

The BANDS table HDU.

format{‘gadf’, ‘fgst-ccube’,’fgst-template’}

FITS format convention.

Returns
wcsWcsGeom

WCS geometry object.

get_coord(idx=None, mode='center', frame=None, sparse=False, axis_name=None)[source]#

Get map coordinates from the geometry.

Parameters
mode{‘center’, ‘edges’}

Get center or edge coordinates for the spatial axes.

framestr or Frame

Coordinate frame

sparsebool

Compute sparse coordinates

axis_namestr

If mode = “edges”, the edges will be returned for this axis.

Returns
coordMapCoord

Map coordinate object.

get_idx(idx=None, flat=False)[source]#

Get tuple of pixel indices for this geometry.

Returns all pixels in the geometry by default. Pixel indices for a single image plane can be accessed by setting idx to the index tuple of a plane.

Parameters
idxtuple, optional

A tuple of indices with one index for each non-spatial dimension. If defined only pixels for the image plane with this index will be returned. If none then all pixels will be returned.

localbool

Flag to return local or global pixel indices. Local indices run from 0 to the number of pixels in a given image plane.

flatbool, optional

Return a flattened array containing only indices for pixels contained in the geometry.

Returns
idxtuple

Tuple of pixel index vectors with one vector for each dimension.

get_pix(idx=None, mode='center')[source]#

Get map pix coordinates from the geometry.

Parameters
mode{‘center’, ‘edges’}

Get center or edge pix coordinates for the spatial axes.

Returns
coordtuple

Map pix coordinate tuple.

is_aligned(other, tolerance=1e-06)[source]#

Check if WCS and extra axes are aligned.

Parameters
otherWcsGeom

Other geom.

tolerancefloat

Tolerance for the comparison.

Returns
alignedbool

Whether geometries are aligned

is_allclose(other, rtol_axes=1e-06, atol_axes=1e-06, rtol_wcs=1e-06)[source]#

Compare two data IRFs for equivalency

Parameters
otherWcsGeom

Geom to compare against

rtol_axesfloat

Relative tolerance for the axes comparison.

atol_axesfloat

Relative tolerance for the axes comparison.

rtol_wcsfloat

Relative tolerance for the wcs comparison.

Returns
is_allclosebool

Whether the geometry is all close.

pad(pad_width, axis_name)#

Pad the geometry at the edges.

Parameters
pad_width{sequence, array_like, int}

Number of values padded to the edges of each axis.

axis_namestr

Name of the axis to pad.

Returns
geomGeom

Padded geometry.

pix_to_coord(pix)[source]#

Convert pixel coordinates to map coordinates.

Parameters
pixtuple

Tuple of pixel coordinates.

Returns
coordstuple

Tuple of map coordinates.

pix_to_idx(pix, clip=False)[source]#

Convert pixel coordinates to pixel indices.

Returns -1 for pixel coordinates that lie outside of the map.

Parameters
pixtuple

Tuple of pixel coordinates.

clipbool

Choose whether to clip indices to the valid range of the geometry. If false then indices for coordinates outside the geometry range will be set -1.

Returns
idxtuple

Tuple of pixel indices.

region_mask(regions, inside=True)[source]#

Create a mask from a given list of regions

The mask is filled such that a pixel inside the region is filled with “True”. To invert the mask, e.g. to create a mask with exclusion regions the tilde (~) operator can be used (see example below).

Parameters
regionsstr, Region or list of Region

Region or list of regions (pixel or sky regions accepted). A region can be defined as a string ind DS9 format as well. See http://ds9.si.edu/doc/ref/region.html for details.

insidebool

For inside=True, pixels in the region to True (the default). For inside=False, pixels in the region are False.

Returns
mask_mapWcsNDMap of boolean type

Boolean region mask

Examples

Make an exclusion mask for a circular region:

from regions import CircleSkyRegion
from astropy.coordinates import SkyCoord, Angle
from gammapy.maps import WcsNDMap, WcsGeom

pos = SkyCoord(0, 0, unit='deg')
geom = WcsGeom.create(skydir=pos, npix=100, binsz=0.1)

region = CircleSkyRegion(
    SkyCoord(3, 2, unit='deg'),
    Angle(1, 'deg'),
)

# the Gammapy convention for exclusion regions is to take the inverse
mask = ~geom.region_mask([region])

Note how we made a list with a single region, since this method expects a list of regions.

region_weights(regions, oversampling_factor=10)[source]#

Compute regions weights

Parameters
regionsstr, Region or list of Region

Region or list of regions (pixel or sky regions accepted). A region can be defined as a string ind DS9 format as well. See http://ds9.si.edu/doc/ref/region.html for details.

oversampling_factorint

Over-sampling factor to compute the region weights

Returns
mapWcsNDMap of boolean type

Weights region mask

rename_axes(names, new_names)#

Rename axes contained in the geometry

Parameters
nameslist or str

Names of the axes.

new_nameslist or str

New names of the axes (list must be of same length than names).

Returns
geomGeom

Renamed geometry.

replace_axis(axis)#

Replace axis with a new one.

Parameters
axisMapAxis

New map axis.

Returns
mapGeom

Geom with replaced axis.

resample_axis(axis)#

Resample geom to a new axis binning.

This method groups the existing bins into a new binning.

Parameters
axisMapAxis

New map axis.

Returns
mapGeom

Geom with resampled axis.

separation(center)[source]#

Compute sky separation wrt a given center.

Parameters
centerSkyCoord

Center position

Returns
separationAngle

Separation angle array (2D)

slice_by_idx(slices)#

Create a new geometry by slicing the non-spatial axes.

Parameters
slicesdict

Dict of axes names and integers or slice object pairs. Contains one element for each non-spatial dimension. For integer indexing the corresponding axes is dropped from the map. Axes not specified in the dict are kept unchanged.

Returns
geomGeom

Sliced geometry.

solid_angle()[source]#

Solid angle array (Quantity in sr).

The array has the same dimension as the WcsGeom object.

To return solid angles for the spatial dimensions only use:

WcsGeom.to_image().solid_angle()
squash(axis_name)#

Squash geom axis.

Parameters
axis_namestr

Axis to squash.

Returns
geomGeom

Geom with squashed axis.

to_bands_hdu(hdu_bands=None, format='gadf')#
to_binsz(binsz)[source]#

Change pixel size of the geometry.

Parameters
binszfloat or tuple or list

New pixel size in degree.

Returns
geomWcsGeom

Geometry with new pixel size.

to_cube(axes)[source]#

Append non-spatial axes to create a higher-dimensional geometry.

This will result in a new geometry with N+M dimensions where N is the number of current dimensions and M is the number of axes in the list.

Parameters
axeslist

Axes that will be appended to this geometry.

Returns
geomGeom

Map geometry.

to_even_npix()[source]#

Create a new geom object with an even number of pixel and a maximum size.

Returns
geomWcsGeom

Geom with odd number of pixels

to_header()[source]#
to_image()[source]#

Create 2D image geometry (drop non-spatial dimensions).

Returns
geomGeom

Image geometry.

to_odd_npix(max_radius=None)[source]#

Create a new geom object with an odd number of pixel and a maximum size.

This is useful for PSF kernel creation.

Parameters
max_radiusQuantity

Max. radius of the geometry (half the width)

Returns
geomWcsGeom

Geom with odd number of pixels

upsample(factor, axis_name=None)[source]#

Upsample the spatial dimension of the geometry by a given factor.

Parameters
factorint

Upsampling factor.

axis_namestr

Axis to upsample.

Returns
geomGeom

Upsampled geometry.