Geometry definitions

The module pyresample.geometry contains classes for describing different kinds of types of remote sensing data geometries. The use of the different classes is described below.

Remarks

All longitudes and latitudes provided to pyresample.geometry must be in degrees. Longitudes must additionally be in the [-180;+180[ validity range.

As of version 1.1.1, the pyresample.geometry contructors will check the range of longitude values, send a warning if some of them fall outside validity range, and automatically correct the invalid values into [-180;+180[.

Use function utils.wrap_longitudes for wrapping longitudes yourself.

AreaDefinition

The cartographic definition of grid areas used by Pyresample is contained in an object of type AreaDefintion. The following arguments are needed to initialize an area:

  • area_id ID of area
  • name: Description
  • proj_id: ID of projection (being deprecated)
  • proj_dict: Proj4 parameters as dict
  • x_size: Number of grid columns
  • y_size: Number of grid rows
  • area_extent: (x_ll, y_ll, x_ur, y_ur)

where

  • x_ll: projection x coordinate of lower left corner of lower left pixel
  • y_ll: projection y coordinate of lower left corner of lower left pixel
  • x_ur: projection x coordinate of upper right corner of upper right pixel
  • y_ur: projection y coordinate of upper right corner of upper right pixel

Creating an area definition:

>>> from pyresample import geometry
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
 >>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> proj_dict = {'a': '6371228.0', 'units': 'm', 'lon_0': '0',
...              'proj': 'laea', 'lat_0': '-90'}
>>> area_def = geometry.AreaDefinition(area_id, description, proj_id,
...                                                                    proj_dict, x_size, y_size, area_extent)
>>> print(area_def)
Area ID: ease_sh
Description: Antarctic EASE grid
Projection ID: ease_sh
Projection: {'a': '6371228.0', 'lat_0': '-90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

pyresample.utils

The utils module of pyresample has convenience functions for constructing area defintions. The function get_area_def can construct an area definition based on area extent and a proj4-string or a list of proj4 arguments.

>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                                       x_size, y_size, area_extent)
>>> print(area_def)
Area ID: ease_sh
Description: Antarctic EASE grid
Projection ID: ease_sh
Projection: {'a': '6371228.0', 'lat_0': '-90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

The load_area function can be used to parse area definitions from a configuration file. Assuming the file areas.yaml exists with the following content

ease_sh:
  description: Antarctic EASE grid
  projection:
    a: 6371228.0
    units: m
    lon_0: 0
    proj: laea
    lat_0: -90
  shape:
    height: 425
    width: 425
  area_extent:
    lower_left_xy: [-5326849.0625, -5326849.0625]
    upper_right_xy: [5326849.0625, 5326849.0625]
    units: m
ease_nh:
  description: Arctic EASE grid
  projection:
    a: 6371228.0
    units: m
    lon_0: 0
    proj: laea
    lat_0: 90
  shape:
    height: 425
    width: 425
  area_extent:
    lower_left_xy: [-5326849.0625, -5326849.0625]
    upper_right_xy: [5326849.0625, 5326849.0625]
    units: m

An area definition dict can be read using

>>> from pyresample import utils
>>> area = utils.load_area('areas.yaml', 'ease_nh')
>>> print(area)
Area ID: ease_nh
Description: Arctic EASE grid
Projection: {'a': '6371228.0', 'lat_0': '90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

Note: In the configuration file, the section name maps to area_id.

Note

The lower_left_xy and upper_right_xy items give the coordinates of the outer edges of the corner pixels on the x and y axis respectively. When the projection coordinates are longitudes and latitudes, it is expected to provide the extent in longitude, latitude order.

Several area definitions can be read at once using the region names in an argument list

>>> from pyresample import utils
>>> nh_def, sh_def = utils.load_area('areas.yaml', 'ease_nh', 'ease_sh')
>>> print(sh_def)
Area ID: ease_sh
Description: Antarctic EASE grid
Projection: {'a': '6371228.0', 'lat_0': '-90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

Note

For backwards compatibility, we still support the legacy area file format:

Assuming the file areas.cfg exists with the following content

REGION: ease_sh {
       NAME:           Antarctic EASE grid
       PCS_ID:         ease_sh
       PCS_DEF:        proj=laea, lat_0=-90, lon_0=0, a=6371228.0, units=m
       XSIZE:          425
       YSIZE:          425
       AREA_EXTENT:    (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
};

REGION: ease_nh {
       NAME:           Arctic EASE grid
       PCS_ID:         ease_nh
       PCS_DEF:        proj=laea, lat_0=90, lon_0=0, a=6371228.0, units=m
       XSIZE:          425
       YSIZE:          425
       AREA_EXTENT:    (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
};

An area definition dict can be read using

>>> from pyresample import utils
>>> area = utils.load_area('areas.cfg', 'ease_nh')
>>> print(area)
Area ID: ease_nh
Description: Arctic EASE grid
Projection ID: ease_nh
Projection: {'a': '6371228.0', 'lat_0': '90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

Note: In the configuration file REGION maps to area_id and PCS_ID maps to proj_id.

Several area definitions can be read at once using the region names in an argument list

>>> from pyresample import utils
>>> nh_def, sh_def = utils.load_area('areas.cfg', 'ease_nh', 'ease_sh')
>>> print(sh_def)
Area ID: ease_sh
Description: Antarctic EASE grid
Projection ID: ease_sh
Projection: {'a': '6371228.0', 'lat_0': '-90', 'lon_0': '0', 'proj': 'laea', 'units': 'm'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

GridDefinition

If the lons and lats grid values are known the area definition information can be skipped for some types of resampling by using a GridDefinition object instead an AreaDefinition object.

>>> import numpy as np
>>> from pyresample import geometry
>>> lons = np.ones((100, 100))
>>> lats = np.ones((100, 100))
>>> grid_def = geometry.GridDefinition(lons=lons, lats=lats)

SwathDefinition

A swath is defined by the lon and lat values of the data points

>>> import numpy as np
>>> from pyresample import geometry
>>> lons = np.ones((500, 20))
>>> lats = np.ones((500, 20))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)

Two swaths can be concatenated if their coloumn count matches

>>> import numpy as np
>>> from pyresample import geometry
>>> lons1 = np.ones((500, 20))
>>> lats1 = np.ones((500, 20))
>>> swath_def1 = geometry.SwathDefinition(lons=lons1, lats=lats1)
>>> lons2 = np.ones((300, 20))
>>> lats2 = np.ones((300, 20))
>>> swath_def2 = geometry.SwathDefinition(lons=lons2, lats=lats2)
>>> swath_def3 = swath_def1.concatenate(swath_def2)

Geographic coordinates and boundaries

A *definition object allows for retrieval of geographic coordinates using array slicing (slice stepping is currently not supported).

All *definition objects exposes the coordinates lons, lats and cartesian_coords. AreaDefinition exposes the full set of projection coordinates as projection_x_coords and projection_y_coords. Note that in the case of projection coordinates expressed in longitude and latitude, projection_x_coords will be longitude and projection_y_coords will be latitude.

Get full coordinate set:

>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                               x_size, y_size, area_extent)
>>> lons, lats = area_def.get_lonlats()

Get slice of coordinate set:

>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                               x_size, y_size, area_extent)
>>> cart_subset = area_def.get_cartesian_coords()[100:200, 350:]

If only the 1D range of a projection coordinate is required it can be extraxted using the proj_x_coord or proj_y_coords property of a geographic coordinate

>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                                       x_size, y_size, area_extent)
>>> proj_x_range = area_def.proj_x_coords

Spherical geometry operations

Some basic spherical operations are available for *definition objects. The spherical geometry operations are calculated based on the corners of a GeometryDefinition (2D SwathDefinition or Grid/AreaDefinition) and assuming the edges are great circle arcs.

It can be tested if geometries overlaps

>>> import numpy as np
>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                                       x_size, y_size, area_extent)
>>> lons = np.array([[-40, -11.1], [9.5, 19.4], [65.5, 47.5], [90.3, 72.3]])
>>> lats = np.array([[-70.1, -58.3], [-78.8, -63.4], [-73, -57.6], [-59.5, -50]])
>>> swath_def = geometry.SwathDefinition(lons, lats)
>>> print(swath_def.overlaps(area_def))
True

The fraction of overlap can be calculated

>>> import numpy as np
>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                                       x_size, y_size, area_extent)
>>> lons = np.array([[-40, -11.1], [9.5, 19.4], [65.5, 47.5], [90.3, 72.3]])
>>> lats = np.array([[-70.1, -58.3], [-78.8, -63.4], [-73, -57.6], [-59.5, -50]])
>>> swath_def = geometry.SwathDefinition(lons, lats)
>>> overlap_fraction = swath_def.overlap_rate(area_def)

And the polygon defining the (great circle) boundaries over the overlapping area can be calculated

>>> import numpy as np
>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                                       x_size, y_size, area_extent)
>>> lons = np.array([[-40, -11.1], [9.5, 19.4], [65.5, 47.5], [90.3, 72.3]])
>>> lats = np.array([[-70.1, -58.3], [-78.8, -63.4], [-73, -57.6], [-59.5, -50]])
>>> swath_def = geometry.SwathDefinition(lons, lats)
>>> overlap_polygon = swath_def.intersection(area_def)

It can be tested if a (lon, lat) point is inside a GeometryDefinition

>>> import numpy as np
>>> from pyresample import utils
>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> x_size = 425
>>> y_size = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = utils.get_area_def(area_id, description, proj_id, projection,
...                                       x_size, y_size, area_extent)
>>> print((0, -90) in area_def)
True