Geometry definitions

The pyresample.geometry module contains classes for describing different geographic areas using a mesh of points or pixels. Some classes represent geographic areas made of of evenly spaced/sized pixels, others handle the cases where the region is described by non-uniform pixels. The best object for describing a region depends on the use case and the information known about it. The different classes available in pyresample are described below.

Note that all longitudes and latitudes provided to pyresample.geometry classes must be in degrees. Additionally, longitudes must be in the [-180;+180[ validity range.

Changed in version 1.8.0: Geometry objects no longer check the validity of the provided longitude and latitude coordinates to improve performance. Longitude arrays are expected to be between -180 and 180 degrees, latitude -90 to 90 degrees. This also applies to all geometry definitions that are provided longitude and latitude arrays on initialization. Use check_and_wrap() to preprocess your arrays.

AreaDefinition

An AreaDefinition, or area, is the primary way of specifying a uniformly spaced geographic region in pyresample. It is also one of the only geometry objects that understands geographic projections. Areas use the PROJ.4 method for describing projected coordinate reference systems (CRS). If the projection for an area is not described by longitude/latitude coordinates then it is typically described in X/Y coordinates in meters. See the PROJ.4 documentation for more information on projections and coordinate reference systems.

The following arguments are needed to initialize an area:

  • area_id: ID of area
  • description: Description
  • proj_id: ID of projection (being deprecated)
  • projection: Proj4 parameters as a dict or string
  • width: Number of grid columns
  • height: Number of grid rows
  • area_extent: (lower_left_x, lower_left_y, upper_right_x, upper_right_y)

where

  • lower_left_x: projection x coordinate of lower left corner of lower left pixel
  • lower_left_y: projection y coordinate of lower left corner of lower left pixel
  • upper_right_x: projection x coordinate of upper right corner of upper right pixel
  • upper_right_y: projection y coordinate of upper right corner of upper right pixel

Example:

>>> from pyresample.geometry import AreaDefinition
>>> 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'}
>>> width = 425
>>> height = 425
>>> area_extent = (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)
>>> area_def
Area ID: ease_sh
Description: Antarctic EASE grid
Projection ID: ease_sh
Projection: {'R': '6371228', 'lat_0': '-90', 'lon_0': '0', 'no_defs': 'None', 'proj': 'laea', 'type': 'crs', 'units': 'm', 'x_0': '0', 'y_0': '0'}
Number of columns: 425
Number of rows: 425
Area extent: (-5326849.0625, -5326849.0625, 5326849.0625, 5326849.0625)

You can also specify the projection using a PROJ.4 string

>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)

or an EPSG code:

>>> projection = '+init=EPSG:3409'  # Use 'EPSG:3409' with pyproj 2.0+
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)

Note

With pyproj 2.0+ please use the new 'EPSG:XXXX' syntax as the old '+init=EPSG:XXXX' is no longer supported.

Creating an AreaDefinition can be complex if you don’t know everything about the region being described. Pyresample provides multiple utilities for creating areas as well as storing them on disk for repeated use. See the Geometry Utilities documentation for more information.

GridDefinition

If the longitude and latitude values for an area are known, the complexity of an AreaDefinition can be skipped by using a GridDefinition object instead. Note that although grid definitions are simpler to define they come at the cost of much higher memory and CPU usage for almost all operations. The longitude and latitude arrays passed to GridDefinition are expected to be evenly spaced. If they are not then a SwathDefinition should be used (see below).

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

SwathDefinition

A swath is defined by the longitude and latitude coordinates for the pixels it represents. The coordinates represent the center point of each pixel. Swaths make no assumptions about the uniformity of pixel size and spacing. This means that operations using then may take longer, but are also accurately represented.

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

Two swaths can be concatenated if their column count matches

>>> lons1 = np.ones((500, 20))
>>> lats1 = np.ones((500, 20))
>>> swath_def1 = SwathDefinition(lons=lons1, lats=lats1)
>>> lons2 = np.ones((300, 20))
>>> lats2 = np.ones((300, 20))
>>> swath_def2 = SwathDefinition(lons=lons2, lats=lats2)
>>> swath_def3 = swath_def1.concatenate(swath_def2)

Geographic coordinates and boundaries

All geometry definition objects provide access to longitude and latitude coordinates. The get_lonlats() method can be used to get this data and will perform any additional calculations needed to get the coordinates.

AreaDefinition exposes the full set of projection coordinates as projection_x_coords and projection_y_coords properties. Note that for lon/lat projections (+proj=latlong) these coordinates will be in longitude/latitude degrees, where projection_x_coords will be longitude and projection_y_coords will be latitude.

Changed in version 1.5.1: Renamed proj_x_coords to projection_x_coords and proj_y_coords to projection_y_coords.

Get longitude and latitude arrays:

>>> 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'
>>> width = 425
>>> height = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)
>>> lons, lats = area_def.get_lonlats()

Get geocentric X, Y, Z coordinates:

>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, 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 extracted using the projection_x_coord or projection_y_coords property of a geographic coordinate

>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)
>>> proj_x_range = area_def.projection_x_coords

Spherical geometry operations

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

Geometries can be checked for overlap:

>>> import numpy as np
>>> 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'
>>> width = 425
>>> height = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, 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 = SwathDefinition(lons, lats)
>>> print(swath_def.overlaps(area_def))
True

The fraction of overlap can be calculated

>>> overlap_fraction = swath_def.overlap_rate(area_def)
>>> overlap_fraction = round(overlap_fraction, 10)
>>> print(overlap_fraction)
0.0584395313

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

>>> overlap_polygon = swath_def.intersection(area_def)
>>> print(overlap_polygon)
[(-40.0, -70.1), (-11.1, -58.3), (72.3, -50.0), (90.3, -59.5)]

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

>>> print((0, -90) in area_def)
True