Pyresample

Pyresample is a Python package for resampling (reprojection) of earth observing satellite data. Pyresample handles both resampling of gridded data (e.g. geostationary satellites) and swath data (polar orbiting satellites). Pyresample can use multiple processor cores for resampling. Pyresample supports masked arrays.

Documentation

Installing Pyresample

Pyresample depends on pyproj, numpy(>= 1.3), scipy(>= 0.7), multiprocessing (builtin package for Python > 2.5) and configobj. Optionally pykdtree can be used instead of scipy from v0.8.0.

The correct version of the packages should be installed on your system (refer to numpy and scipy installation instructions) or use easy_install to handle dependencies automatically.

In order to use the pyresample plotting functionality Basemap and matplotlib (>= 0.98) must be installed. These packages are not a prerequisite for using any other pyresample functionality.

Package test

Test the package (requires nose):

$ tar -zxvf pyresample-<version>.tar.gz
$ cd pyresample-<version>
$ nosetests

If all the tests passes the functionality of all pyresample functions on the system has been verified.

Package installation

A sandbox environment can be created for pyresample using Virtualenv

Pyresample is available from pypi.

Install Pyresample using pip:

$ pip install pyresample

Alternatively install from tarball:

$ tar -zxvf pyresample-<version>.tar.gz
$ cd pyresample-<version>
$ python setup.py install

Using pykdtree

As of pyresample v0.8.0 pykdtree can be used as backend instead of scipy. This enables significant speedups for large datasets.

pykdtree is used as a drop-in replacement for scipy. If it’s available it will be used otherwise scipy will be used. To check which backend is active for your pyresample installation do:

>>> import pyresample as pr
>>> pr.kd_tree.which_kdtree()

which returns either ‘pykdtree’ or ‘scipy.spatial’.

Please refere to pykdtree for installation description.

If pykdtree is built with OpenMP support the number of threads is controlled with the standard OpenMP environment variable OMP_NUM_THREADS. The nprocs argument has no effect on pykdtree.

Using numexpr

As of pyresample v1.0.0 numexpr will be used for minor bottleneck optimization if available

Show active plugins

The active drop-in plugins can be show using:

>>> import pyresample as pr
>>> pr.get_capabilities()

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

Geographic filtering

The module pyresample.geo_filter contains classes to filter geo data

GridFilter

Allows for filtering of data based on a geographic mask. The filtering uses a bucket sampling approach.

The following example shows how to select data falling in the upper left and lower right quadrant of a full globe Plate Carrée projection using an 8x8 filter mask

>>> import numpy as np
>>> from pyresample import geometry, geo_filter
>>> lons = np.array([-170, -30, 30, 170])
>>> lats = np.array([20, -40, 50, -80])
>>> swath_def = geometry.SwathDefinition(lons, lats)
>>> data = np.array([1, 2, 3, 4])
>>> filter_area = geometry.AreaDefinition('test', 'test', 'test',
...         {'proj' : 'eqc', 'lon_0' : 0.0, 'lat_0' : 0.0},
...           8, 8,
...          (-20037508.34, -10018754.17, 20037508.34, 10018754.17)
...             )
>>> filter = np.array([[1, 1, 1, 1, 0, 0, 0, 0],
...         [1, 1, 1, 1, 0, 0, 0, 0],
...         [1, 1, 1, 1, 0, 0, 0, 0],
...         [1, 1, 1, 1, 0, 0, 0, 0],
...         [0, 0, 0, 0, 1, 1, 1, 1],
...         [0, 0, 0, 0, 1, 1, 1, 1],
...         [0, 0, 0, 0, 1, 1, 1, 1],
...         [0, 0, 0, 0, 1, 1, 1, 1],
...         ])
>>> grid_filter = geo_filter.GridFilter(filter_area, filter)
>>> swath_def_filtered, data_filtered = grid_filter.filter(swath_def, data)

Input swath_def and data must match as described in Resampling of swath data.

The returned data will always have a 1D geometry_def and if multiple channels are present the filtered data will have the shape (number_of_points, channels).

Resampling of gridded data

Pyresample can be used to resample from an existing grid to another. Nearest neighbour resampling is used.

pyresample.image

A grid can be stored in an object of type ImageContainer along with its area definition. An object of type ImageContainer allows for calculating resampling using preprocessed arrays using the method get_array_from_linesample

Resampling can be done using descendants of ImageContainer and calling their resample method.

An ImageContainerQuick object allows for the grid to be resampled to a new area defintion using an approximate (but fast) nearest neighbour method. Resampling an object of type ImageContainerQuick returns a new object of type ImageContainerQuick.

An ImageContainerNearest object allows for the grid to be resampled to a new area defintion (or swath definition) using an accurate kd-tree method. Resampling an object of type ImageContainerNearest returns a new object of type ImageContainerNearest.

>>> import numpy as np
>>> from pyresample import image, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> msg_area = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
...                                'msg_full',
...                                {'a': '6378169.0', 'b': '6356584.0',
...                                 'h': '35785831.0', 'lon_0': '0',
...                                 'proj': 'geos'},
...                                3712, 3712,
...                                [-5568742.4, -5568742.4,
...                                 5568742.4, 5568742.4])
>>> data = np.ones((3712, 3712))
>>> msg_con_quick = image.ImageContainerQuick(data, msg_area)
>>> area_con_quick = msg_con_quick.resample(area_def)
>>> result_data_quick = area_con_quick.image_data
>>> msg_con_nn = image.ImageContainerNearest(data, msg_area, radius_of_influence=50000)
>>> area_con_nn = msg_con_nn.resample(area_def)
>>> result_data_nn = area_con_nn.image_data

Data is assumed to be a numpy array of shape (rows, cols) or (rows, cols, channels).

Masked arrays can be used as data input. In order to have undefined pixels masked out instead of assigned a fill value set fill_value=None when calling resample_area_*.

Using ImageContainerQuick the risk of image artifacts increases as the distance from source projection center increases.

The constructor argument radius_of_influence to ImageContainerNearest specifices the maximum distance to search for a neighbour for each point in the target grid. The unit is meters.

The constructor arguments of an ImageContainer object can be changed as attributes later

>>> import numpy as np
>>> from pyresample import image, geometry
>>> msg_area = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
...                                'msg_full',
...                                {'a': '6378169.0', 'b': '6356584.0',
...                                 'h': '35785831.0', 'lon_0': '0',
...                                 'proj': 'geos'},
...                                3712, 3712,
...                                [-5568742.4, -5568742.4,
...                                 5568742.4, 5568742.4])
>>> data = np.ones((3712, 3712))
>>> msg_con_nn = image.ImageContainerNearest(data, msg_area, radius_of_influence=50000)
>>> msg_con_nn.radius_of_influence = 45000
>>> msg_con_nn.fill_value = -99
Multi channel images

If the dataset has several channels the last index of the data array specifies the channels

>>> import numpy as np
>>> from pyresample import image, geometry
>>> msg_area = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
...                                'msg_full',
...                                {'a': '6378169.0', 'b': '6356584.0',
...                                 'h': '35785831.0', 'lon_0': '0',
...                                 'proj': 'geos'},
...                                3712, 3712,
...                                [-5568742.4, -5568742.4,
...                                 5568742.4, 5568742.4])
>>> channel1 = np.ones((3712, 3712))
>>> channel2 = np.ones((3712, 3712)) * 2
>>> channel3 = np.ones((3712, 3712)) * 3
>>> data = np.dstack((channel1, channel2, channel3))
>>> msg_con_nn = image.ImageContainerNearest(data, msg_area, radius_of_influence=50000)
Segmented resampling

Pyresample calculates the result in segments in order to reduce memory footprint. This is controlled by the segments contructor keyword argument. If no segments argument is given pyresample will estimate the number of segments to use.

Forcing quick resampling to use 4 resampling segments:

>>> import numpy as np
>>> from pyresample import image, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> msg_area = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
...                                'msg_full',
...                                {'a': '6378169.0', 'b': '6356584.0',
...                                 'h': '35785831.0', 'lon_0': '0',
...                                 'proj': 'geos'},
...                                3712, 3712,
...                                [-5568742.4, -5568742.4,
...                                 5568742.4, 5568742.4])
>>> data = np.ones((3712, 3712))
>>> msg_con_quick = image.ImageContainerQuick(data, msg_area, segments=4)
>>> area_con_quick = msg_con_quick.resample(area_def)
Constructor arguments

The full list of constructor arguments:

ImageContainerQuick:
  • image_data : Dataset. Masked arrays can be used.
  • image_data : Geometry definition.
  • fill_value (optional) : Fill value for undefined pixels. Defaults to 0. If set to None they will be masked out.
  • nprocs (optional) : Number of processor cores to use. Defaults to 1.
  • segments (optional) : Number of segments to split resampling in. Defaults to auto estimation.
ImageContainerNearest:
  • image_data : Dataset. Masked arrays can be used.
  • image_data : Geometry definition.
  • radius_of_influence : Cut off radius in meters when considering neighbour pixels.
  • epsilon (optional) : The distance to a found value is guaranteed to be no further than (1 + eps) times the distance to the correct neighbour.
  • fill_value (optional) : Fill value for undefined pixels. Defaults to 0. If set to None they will be masked out.
  • reduce_data (optional) : Apply geographic reduction of dataset before resampling. Defaults to True
  • nprocs (optional) : Number of processor cores to use. Defaults to 1.
  • segments (optional) : Number of segments to split resampling in. Defaults to auto estimation.
Preprocessing of grid resampling

For preprocessing of grid resampling see Preprocessing of grids

Resampling of swath data

Pyresample can be used to resample a swath dataset to a grid, a grid to a swath or a swath to another swath. Resampling can be done using nearest neighbour method, Guassian weighting, weighting with an arbitrary radial function.

pyresample.image

The ImageContainerNearest class can be used for nearest neighbour resampling of swaths as well as grids.

>>> import numpy as np
>>> from pyresample import image, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> swath_con = image.ImageContainerNearest(data, swath_def, radius_of_influence=5000)
>>> area_con = swath_con.resample(area_def)
>>> result = area_con.image_data

For other resampling types or splitting the process in two steps use the functions in pyresample.kd_tree described below.

pyresample.kd_tree

This module contains several functions for resampling swath data.

Note distance calculation is approximated with cartesian distance.

Masked arrays can be used as data input. In order to have undefined pixels masked out instead of assigned a fill value set fill_value=None when calling the resample_* function.

resample_nearest

Function for resampling using nearest neighbour method.

Example showing how to resample a generated swath dataset to a grid using nearest neighbour method:

>>> import numpy as np
>>> from pyresample import kd_tree, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> result = kd_tree.resample_nearest(swath_def, data,
... area_def, radius_of_influence=50000, epsilon=0.5)

If the arguments swath_def and area_def where switched (and data matched the dimensions of area_def) the grid of area_def would be resampled to the swath defined by swath_def.

Note the keyword arguments:

  • radius_of_influence: The radius around each grid pixel in meters to search for neighbours in the swath.
  • epsilon: The distance to a found value is guaranteed to be no further than (1 + eps) times the distance to the correct neighbour. Allowing for uncertanty decreases execution time.

If data is a masked array the mask will follow the neighbour pixel assignment.

If there are multiple channels in the dataset the data argument should be of the shape of the lons and lat arrays with the channels along the last axis e.g. (rows, cols, channels). Note: the convention of pyresample < 0.7.4 is to pass data in the form of (number_of_data_points, channels) is still accepted.

>>> import numpy as np
>>> from pyresample import kd_tree, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> channel1 = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> channel2 = np.fromfunction(lambda y, x: y*x, (50, 10)) * 2
>>> channel3 = np.fromfunction(lambda y, x: y*x, (50, 10)) * 3
>>> data = np.dstack((channel1, channel2, channel3))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> result = kd_tree.resample_nearest(swath_def, data,
... area_def, radius_of_influence=50000)

For nearest neighbour resampling the class image.ImageContainerNearest can be used as well as kd_tree.resample_nearest

resample_gauss

Function for resampling using nearest Gussian weighting. The Gauss weigh function is defined as exp(-dist^2/sigma^2). Note the pyresample sigma is not the standard deviation of the gaussian. Example showing how to resample a generated swath dataset to a grid using Gaussian weighting:

>>> import numpy as np
>>> from pyresample import kd_tree, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> result = kd_tree.resample_gauss(swath_def, data,
... area_def, radius_of_influence=50000, sigmas=25000)

If more channels are present in data the keyword argument sigmas must be a list containing a sigma for each channel.

If data is a masked array any pixel in the result data that has been “contaminated” by weighting of a masked pixel is masked.

Using the function utils.fwhm2sigma the sigma argument to the gauss resampling can be calculated from 3 dB FOV levels.

resample_custom

Function for resampling using arbitrary radial weight functions.

Example showing how to resample a generated swath dataset to a grid using an arbitrary radial weight function:

>>> import numpy as np
>>> from pyresample import kd_tree, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> wf = lambda r: 1 - r/100000.0
>>> result  = kd_tree.resample_custom(swath_def, data,
...  area_def, radius_of_influence=50000, weight_funcs=wf)

If more channels are present in data the keyword argument weight_funcs must be a list containing a radial function for each channel.

If data is a masked array any pixel in the result data that has been “contaminated” by weighting of a masked pixel is masked.

Uncertainty estimates

Uncertainty estimates in the form of weighted standard deviation can be obtained from the resample_custom and resample_gauss functions. By default the functions return the result of the resampling as a single numpy array. If the functions are given the keyword argument with_uncert=True then the following list of numpy arrays will be returned instead: (result, stddev, count). result is the usual result. stddev is the weighted standard deviation for each element in the result. count is the number of data values used in the weighting for each element in the result.

The principle is to view the calculated value for each element in the result as a weighted average of values sampled from a statistical variable. An estimate of the standard deviation of the distribution is calculated using the unbiased weighted estimator given as stddev = sqrt((V1 / (V1 ** 2 + V2)) * sum(wi * (xi - result) ** 2)) where result is the result of the resampling. xi is the value of a contributing neighbour and wi is the corresponding weight. The coefficients are given as V1 = sum(wi) and V2 = sum(wi ** 2). The standard deviation is only calculated for elements in the result where more than one neighbour has contributed to the weighting. The count numpy array can be used for filtering at a higher number of contributing neigbours.

Usage only differs in the number of return values from resample_gauss and resample_custom. E.g.:

>>> result, stddev, count = pr.kd_tree.resample_gauss(swath_def, ice_conc, area_def,
...                                                   radius_of_influence=20000,
...                                                   sigmas=pr.utils.fwhm2sigma(35000),
...                                                   fill_value=None, with_uncert=True)
Below is shown a plot of the result of the resampling using a real data set:
_images/uncert_conc_nh.png
The corresponding standard deviations:
_images/uncert_stddev_nh.png
And the number of contributing neighbours for each element:
_images/uncert_count_nh.png

Notice the standard deviation is only calculated where there are more than one contributing neighbour.

Resampling from neighbour info

The resampling can be split in two steps:

First get arrays containing information about the nearest neighbours to each grid point. Then use these arrays to retrive the resampling result.

This approch can be useful if several datasets based on the same swath are to be resampled. The computational heavy task of calculating the neighbour information can be done once and the result can be used to retrieve the resampled data from each of the datasets fast.

>>> import numpy as np
>>> from pyresample import kd_tree, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> valid_input_index, valid_output_index, index_array, distance_array = \
...                        kd_tree.get_neighbour_info(swath_def,
...                                                       area_def, 50000,
...                                                   neighbours=1)
>>> res = kd_tree.get_sample_from_neighbour_info('nn', area_def.shape, data,
...                                              valid_input_index, valid_output_index,
...                                              index_array)

Note the keyword argument neighbours=1. This specifies only to consider one neighbour for each grid point (the nearest neighbour). Also note distance_array is not a required argument for get_sample_from_neighbour_info when using nearest neighbour resampling

Segmented resampling

Whenever a resampling function takes the keyword argument segments the number of segments to split the resampling process in can be specified. This affects the memory footprint of pyresample. If the value of segments is left to default pyresample will estimate the number of segments to use.

Speedup using pykdtree

pykdtree can be used instead of scipy to gain significant speedup for large datasets. See Using multiple processor cores.

pyresample.bilinear

Compared to nearest neighbour resampling, bilinear interpolation produces smoother results near swath edges of polar satellite data and edges of geostationary satellites.

The algorithm is implemented from http://www.ahinson.com/algorithms_general/Sections/InterpolationRegression/InterpolationIrregularBilinear.pdf

Below is shown a comparison between image generated with nearest neighbour resampling (top) and with bilinear interpolation (bottom):

_images/nearest_overview.png _images/bilinear_overview.png

Click images to see the full resolution versions.

The perceived sharpness of the bottom image is lower, but there is more detail present.

resample_bilinear

Function for resampling using bilinear interpolation for irregular source grids.

>>> import numpy as np
>>> from pyresample import bilinear, geometry
>>> target_def = geometry.AreaDefinition('areaD',
...                                      'Europe (3km, HRV, VTC)',
...                                      'areaD',
...                                      {'a': '6378144.0', 'b': '6356759.0',
...                                       'lat_0': '50.00', 'lat_ts': '50.00',
...                                       'lon_0': '8.00', 'proj': 'stere'},
...                                      800, 800,
...                                      [-1370912.72, -909968.64,
...                                       1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> source_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> result = bilinear.resample_bilinear(data, source_def, target_def,
...                                     radius=50e3, neighbours=32,
...                                     nprocs=1, fill_value=0,
...                                     reduce_data=True, segments=None,
...                                     epsilon=0)

The target_area needs to be an area definition with proj4_string attribute.

Keyword arguments which are passed to kd_tree:

  • radius: radius around each target pixel in meters to search for neighbours in the source data
  • neighbours: number of closest locations to consider when selecting the four data points around the target pixel
  • nprocs: number of processors to use for finding the closest pixels
  • fill_value: fill invalid pixel with this value. If fill_value=None is used, masked arrays will be returned
  • reduce_data: do/don’t do preliminary data reduction before calculating the neigbour info
  • segments: number of segments to use in neighbour search
  • epsilon: maximum uncertainty allowed in neighbour search

The example above shows the default value for each keyword argument.

Resampling from bilinear coefficients

As for nearest neighbour resampling, also bilinear interpolation can be split in two steps.

  • Calculate interpolation coefficients, input data reduction matrix and mapping matrix
  • Use this information to resample several datasets between these two areas/swaths

Only the first step is computationally expensive operation, so by re-using this information the overall processing time is reduced significantly. This is also done internally by the resample_bilinear function, but separating these steps makes it possible to cache the coefficients if the same transformation is done over and over again. This is very typical in operational geostationary satellite image processing.

>>> import numpy as np
>>> from pyresample import bilinear, geometry
>>> target_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)',
...                                      'areaD',
...                                      {'a': '6378144.0', 'b': '6356759.0',
...                                       'lat_0': '50.00', 'lat_ts': '50.00',
...                                       'lon_0': '8.00', 'proj': 'stere'},
...                                      800, 800,
...                                      [-1370912.72, -909968.64,
...                                       1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> source_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> t_params, s_params, input_idxs, idx_ref = \
...     bilinear.get_bil_info(source_def, target_def, radius=50e3, nprocs=1)
>>> res = bilinear.get_sample_from_bil_info(data.ravel(), t_params, s_params,
...                                         input_idxs, idx_ref)

pyresample.ewa

Pyresample makes it possible to resample swath data to a uniform grid using an Elliptical Weighted Averaging algorithm or EWA for short. This algorithm behaves differently than the KDTree based resampling algorithms that pyresample provides. The KDTree-based algorithms process each output grid pixel by searching for all “nearby” input pixels and applying a certain interpolation (nearest neighbor, gaussian, etc). The EWA algorithm processes each input pixel mapping it to one or more output pixels. Once each input pixel has been analyzed the intermediate results are averaged to produce the final gridded result.

The EWA algorithm also has limitations on how the input data is structured compared to the generic KDTree algorithms. EWA assumes that data in the array is organized geographically; adjacent data in the array is adjacent data geographically. The algorithm uses this to configure parameters based on the size and location of the swath pixels.

The EWA algorithm consists of two steps: ll2cr and fornav. The algorithm was originally part of the MODIS Swath to Grid Toolbox (ms2gt) created by the NASA National Snow & Ice Data Center (NSIDC). Its default parameters work best with MODIS L1B data, but it has been proven to produce high quality images from VIIRS and AVHRR data with the right parameters.

Note

This code was originally part of the Polar2Grid project. This documentation and the API documentation for this algorithm may still use references or concepts from Polar2Grid until everything can be updated.

Gridding

The first step is called ‘ll2cr’ which stands for “longitude/latitude to column/row”. This step maps the pixel location (lon/lat space) into area (grid) space. Areas in pyresample are defined by a PROJ.4 projection specification. An area is defined by the following parameters:

  • Grid Name
  • PROJ.4 String (either lat/lon or metered projection space)
  • Width (number of pixels in the X direction)
  • Height (number of pixels in the Y direction)
  • Cell Width (pixel size in the X direction in grid units)
  • Cell Height (pixel size in the Y direction in grid units)
  • X Origin (upper-left X coordinate in grid units)
  • Y Origin (upper-left Y coordinate in grid units)
Resampling

The second step of EWA remapping is called “fornav”, short for “forward navigation”. This EWA algorithm processes one input scan line at a time. The algorithm weights the effect of an input pixel on an output pixel based on its location in the scan line and other calculated coefficients. It can also handle swaths that are not scan based by specifying rows_per_scan as the number of rows in the entire swath. How the algorithm treats the data can be configured with various keyword arguments, see the API documentation for more information. Both steps provide additional information to inform the user how much data was used in the result. The first returned value of ll2cr tells you how many of the input swath pixels overlap the grid. The first returned value of fornav tells you how many grid points have valid data values in them.

Example

Note

EWA resampling in pyresample is still in an alpha stage. As development continues, EWA resampling may be called differently.

>>> import numpy as np
>>> from pyresample.ewa import ll2cr, fornav
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> # ll2cr converts swath longitudes and latitudes to grid columns and rows
>>> swath_points_in_grid, cols, rows = ll2cr(swath_def, area_def)
>>> # if the data is scan based, specify how many data rows make up one scan
>>> rows_per_scan = 5
>>> # fornav resamples the swath data to the gridded area
>>> num_valid_points, gridded_data = fornav(cols, rows, area_def, data, rows_per_scan=rows_per_scan)

Using multiple processor cores

Multi core processing

Bottlenecks of pyresample can be executed in parallel. Parallel computing can be executed if the pyresample function has the nprocs keyword argument. nprocs specifies the number of processes to be used for calculation. If a class takes the constructor argument nprocs this sets nprocs for all methods of this class

Example of resampling in parallel using 4 processes:

>>> import numpy
>>> from pyresample import kd_tree, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = numpy.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = numpy.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = numpy.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> swath_def = geometry.SwathDefinition(lons=lons, lats=lats)
>>> result = kd_tree.resample_nearest(swath_def, data.ravel(),
... area_def, radius_of_influence=50000, nprocs=4)

Note: Do not use more processes than available processor cores. As there is a process creation overhead there might be neglible performance improvement using say 8 compared to 4 processor cores. Test on the actual system to determine the most sensible number of processes to use.

Here is an example of the performance for a varying number of processors on a 64-bit ubuntu 14.04, 32 GB RAM, 2 x Intel Xeon E5-2630 with 6 physical cores each:
_images/time_vs_nproc_1-12.png

Preprocessing of grids

When resampling is performed repeatedly to the same grid significant execution time can be save by preprocessing grid information.

Preprocessing for grid resampling

Using the function generate_quick_linesample_arrays or generate_nearest_neighbour_linesample_arrays from pyresample.utils arrays containing the rows and cols indices used to calculate the result in image.resample_area_quick or resample_area_nearest_neighbour can be obtained. These can be fed to the method get_array_from_linesample of an ImageContainer object to obtain the resample result.

>>> import numpy
>>> from pyresample import utils, image, geometry
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> msg_area = geometry.AreaDefinition('msg_full', 'Full globe MSG image 0 degrees',
...                                'msg_full',
...                                {'a': '6378169.0', 'b': '6356584.0',
...                                 'h': '35785831.0', 'lon_0': '0',
...                                 'proj': 'geos'},
...                                3712, 3712,
...                                [-5568742.4, -5568742.4,
...                                 5568742.4, 5568742.4])
>>> data = numpy.ones((3712, 3712))
>>> msg_con = image.ImageContainer(data, msg_area)
>>> row_indices, col_indices = \
...            utils.generate_nearest_neighbour_linesample_arrays(msg_area, area_def, 50000)
>>> result = msg_con.get_array_from_linesample(row_indices, col_indices)

The numpy arrays returned by generate_*_linesample_arrays can be and used with the ImageContainer.get_array_from_linesample method when the same resampling is to be performed again thus eliminating the need for calculating the reprojection.

Numpy arrays can be saved and loaded using numpy.save and numpy.load.

Plotting with pyresample and Basemap

Pyresample supports basic integration with Basemap (http://matplotlib.sourceforge.net/basemap).

Displaying data quickly

Pyresample has some convenience functions for displaying data from a single channel. The function plot.show_quicklook shows a Basemap image of a dataset for a specified AreaDefinition. The function plot.save_quicklook saves the Basemap image directly to file.

Example usage:

>>> import numpy as np
>>> import pyresample as pr
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = pr.utils.load_area('areas.cfg', 'ease_sh')
>>> swath_def = pr.geometry.SwathDefinition(lons, lats)
>>> result = pr.kd_tree.resample_nearest(swath_def, tb37v, area_def,
...                                      radius_of_influence=20000, fill_value=None)
>>> pr.plot.save_quicklook('tb37v_quick.png', area_def, result, label='Tb 37v (K)')
Assuming lons, lats and tb37v are initialized with real data the result might look something like this:
_images/tb37v_quick.png

The data passed to the functions is a 2D array matching the AreaDefinition.

The Plate Carree projection

The Plate Carree projection (regular lon/lat grid) is named eqc in Proj.4 and cyl in Basemap. pyresample uses the Proj.4 name. Assuming the file areas.cfg has the following area definition:

REGION: pc_world {
   NAME:    Plate Carree world map
   PCS_ID:  pc_world
   PCS_DEF: proj=eqc
   XSIZE: 640
   YSIZE: 480
   AREA_EXTENT:  (-20037508.34, -10018754.17, 20037508.34, 10018754.17)
};

Example usage:

>>> import numpy as np
>>> import pyresample as pr
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = pr.utils.load_area('areas.cfg', 'pc_world')
>>> swath_def = pr.geometry.SwathDefinition(lons, lats)
>>> result = pr.kd_tree.resample_nearest(swath_def, tb37v, area_def, radius_of_influence=20000, fill_value=None)
>>> pr.plot.save_quicklook('tb37v_pc.png', area_def, result, num_meridians=0, num_parallels=0, label='Tb 37v (K)')
Assuming lons, lats and tb37v are initialized with real data the result might look something like this:
_images/tb37v_pc.png
The Globe projections

From v0.7.12 pyresample can use the geos, ortho and nsper projections with Basemap. Assuming the file areas.cfg has the following area definition for an ortho projection area:

REGION: ortho {
  NAME:    Ortho globe
  PCS_ID:  ortho_globe
  PCS_DEF: proj=ortho, a=6370997.0, lon_0=40, lat_0=-40
  XSIZE: 640
  YSIZE: 480
  AREA_EXTENT:  (-10000000, -10000000, 10000000, 10000000)
};

Example usage:

>>> import numpy as np
>>> import pyresample as pr
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = pr.utils.load_area('areas.cfg', 'ortho')
>>> swath_def = pr.geometry.SwathDefinition(lons, lats)
>>> result = pr.kd_tree.resample_nearest(swath_def, tb37v, area_def, radius_of_influence=20000, fill_value=None)
>>> pr.plot.save_quicklook('tb37v_ortho.png', area_def, result, num_meridians=0, num_parallels=0, label='Tb 37v (K)')
Assuming lons, lats and tb37v are initialized with real data the result might look something like this:
_images/tb37v_ortho.png

Getting a Basemap object

In order to make more advanced plots than the preconfigured quicklooks a Basemap object can be generated from an AreaDefintion using the plot.area_def2basemap(area_def, **kwargs) function.

Example usage:

>>> import numpy as np
>>> import matplotlib.pyplot as plt
>>> import pyresample as pr
>>> lons = np.zeros(1000)
>>> lats = np.arange(-80, -90, -0.01)
>>> tb37v = np.arange(1000)
>>> area_def = pr.utils.load_area('areas.cfg', 'ease_sh')
>>> swath_def = pr.geometry.SwathDefinition(lons, lats)
>>> result = pr.kd_tree.resample_nearest(swath_def, tb37v, area_def,
...                                      radius_of_influence=20000, fill_value=None)
>>> bmap = pr.plot.area_def2basemap(area_def)
>>> bmng = bmap.bluemarble()
>>> col = bmap.imshow(result, origin='upper')
>>> plt.savefig('tb37v_bmng.png', bbox_inches='tight')
Assuming lons, lats and tb37v are initialized with real data the result might look something like this:
_images/tb37v_bmng.png

Any keyword arguments (not concerning the projection) passed to plot.area_def2basemap will be passed directly to the Basemap initialization.

For more information on how to plot with Basemap please refer to the Basemap and matplotlib documentation.

Limitations

The pyresample use of Basemap is basically a conversion from a pyresample AreaDefintion to a Basemap object which allows for correct plotting of a resampled dataset using the basemap.imshow function.

Currently only the following set of Proj.4 arguments can be interpreted in the conversion: {‘proj’, ‘a’, ‘b’, ‘ellps’, ‘lon_0’, ‘lat_0’, ‘lon_1’, ‘lat_1’, ‘lon_2’, ‘lat_2’, ‘lat_ts’}

Any other Proj.4 parameters will be ignored. If the ellipsoid is not defined in terms of ‘ellps’, ‘a’ or (‘a’, ‘b’) it will default to WGS84.

The xsize and ysize in an AreaDefinition will only be used during resampling when the image data for use in basemap.imshow is created. The actual size and shape of the final plot is handled by matplotlib.

Reduction of swath data

Given a swath and a cartesian grid or grid lons and lats pyresample can reduce the swath data to only the relevant part covering the grid area. The reduction is coarse in order not to risk removing relevant data.

From data_reduce the function swath_from_lonlat_grid can be used to reduce the swath data set to the area covering the lon lat grid

>>> import numpy as np
>>> from pyresample import geometry, data_reduce
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> grid_lons, grid_lats = area_def.get_lonlats()
>>> reduced_lons, reduced_lats, reduced_data = \
...                            data_reduce.swath_from_lonlat_grid(grid_lons, grid_lats,
...                            lons, lats, data,
...                            radius_of_influence=3000)

radius_of_influence is used to calculate a buffer zone around the grid where swath data points are not reduced.

The function get_valid_index_from_lonlat_grid returns a boolean array of same size as the swath indicating the relevant swath data points compared to the grid

>>> import numpy as np
>>> from pyresample import geometry, data_reduce
>>> area_def = geometry.AreaDefinition('areaD', 'Europe (3km, HRV, VTC)', 'areaD',
...                                {'a': '6378144.0', 'b': '6356759.0',
...                                 'lat_0': '50.00', 'lat_ts': '50.00',
...                                 'lon_0': '8.00', 'proj': 'stere'},
...                                800, 800,
...                                [-1370912.72, -909968.64,
...                                 1029087.28, 1490031.36])
>>> data = np.fromfunction(lambda y, x: y*x, (50, 10))
>>> lons = np.fromfunction(lambda y, x: 3 + x, (50, 10))
>>> lats = np.fromfunction(lambda y, x: 75 - y, (50, 10))
>>> grid_lons, grid_lats = area_def.get_lonlats()
>>> valid_index = data_reduce.get_valid_index_from_lonlat_grid(grid_lons, grid_lats,
...                                            lons, lats,
...                                            radius_of_influence=3000)

pyresample API

pyresample.geometry

Classes for geometry operations

class geometry.AreaDefinition(area_id, name, proj_id, proj_dict, x_size, y_size, area_extent, nprocs=1, lons=None, lats=None, dtype=<Mock object>)

Holds definition of an area.

Parameters:
  • area_id (str) – ID of area
  • name (str) – Name of area
  • proj_id (str) – ID of projection
  • proj_dict (dict) – Dictionary with Proj.4 parameters
  • x_size (int) – x dimension in number of pixels
  • y_size (int) – y dimension in number of pixels
  • area_extent (list) – Area extent as a list (LL_x, LL_y, UR_x, UR_y)
  • nprocs (int, optional) – Number of processor cores to be used
  • lons (numpy array, optional) – Grid lons
  • lats (numpy array, optional) – Grid lats
area_id

str – ID of area

name

str – Name of area

proj_id

str – ID of projection

proj_dict

dict – Dictionary with Proj.4 parameters

x_size

int – x dimension in number of pixels

y_size

int – y dimension in number of pixels

shape

tuple – Corresponding array shape as (rows, cols)

size

int – Number of points in grid

area_extent

tuple – Area extent as a tuple (LL_x, LL_y, UR_x, UR_y)

area_extent_ll

tuple – Area extent in lons lats as a tuple (LL_lon, LL_lat, UR_lon, UR_lat)

pixel_size_x

float – Pixel width in projection units

pixel_size_y

float – Pixel height in projection units

pixel_upper_left

list – Coordinates (x, y) of center of upper left pixel in projection units

pixel_offset_x

float – x offset between projection center and upper left corner of upper left pixel in units of pixels.

pixel_offset_y

float – y offset between projection center and upper left corner of upper left pixel in units of pixels..

proj4_string

str – Projection defined as Proj.4 string

lons

object – Grid lons

lats

object – Grid lats

cartesian_coords

object – Grid cartesian coordinates

projection_x_coords

object – Grid projection x coordinate

projection_y_coords

object – Grid projection y coordinate

colrow2lonlat(cols, rows)

Return longitudes and latitudes for the given image columns and rows. Both scalars and arrays are supported. To be used with scarse data points instead of slices (see get_lonlats).

get_lonlat(row, col)

Retrieves lon and lat values of single point in area grid

Parameters:
  • row (int) –
  • col (int) –
Returns:

(lon, lat)

Return type:

tuple of floats

get_lonlats(nprocs=None, data_slice=None, cache=False, dtype=None)

Returns lon and lat arrays of area.

Parameters:
  • nprocs (int, optional) – Number of processor cores to be used. Defaults to the nprocs set when instantiating object
  • data_slice (slice object, optional) – Calculate only coordinates for specified slice
  • cache (bool, optional) – Store result the result. Requires data_slice to be None
Returns:

(lons, lats) – Grids of area lons and and lats

Return type:

tuple of numpy arrays

get_proj_coords(data_slice=None, cache=False, dtype=None)

Get projection coordinates of grid

Parameters:
  • data_slice (slice object, optional) – Calculate only coordinates for specified slice
  • cache (bool, optional) – Store result the result. Requires data_slice to be None
Returns:

(target_x, target_y) – Grids of area x- and y-coordinates in projection units

Return type:

tuple of numpy arrays

get_xy_from_lonlat(lon, lat)

Retrieve closest x and y coordinates (column, row indices) for the specified geolocation (lon,lat) if inside area. If lon,lat is a point a ValueError is raised if the return point is outside the area domain. If lon,lat is a tuple of sequences of longitudes and latitudes, a tuple of masked arrays are returned.

Input:

lon : point or sequence (list or array) of longitudes lat : point or sequence (list or array) of latitudes

Returns:

(x, y) : tuple of integer points/arrays

get_xy_from_proj_coords(xm_, ym_)

Retrieve closest x and y coordinates (column, row indices) for a location specified with projection coordinates (xm_,ym_) in meters. A ValueError is raised, if the return point is outside the area domain. If xm_,ym_ is a tuple of sequences of projection coordinates, a tuple of masked arrays are returned.

Input:

xm_ : point or sequence (list or array) of x-coordinates in m (map projection) ym_ : point or sequence (list or array) of y-coordinates in m (map projection)

Returns:

(x, y) : tuple of integer points/arrays

lonlat2colrow(lons, lats)

Return image columns and rows for the given longitudes and latitudes. Both scalars and arrays are supported. Same as get_xy_from_lonlat, renamed for convenience.

outer_boundary_corners

Returns the lon,lat of the outer edges of the corner points

proj4_string

Returns projection definition as Proj.4 string

class geometry.BaseDefinition(lons=None, lats=None, nprocs=1)

Base class for geometry definitions

corners

Returns the corners of the current area.

get_area()

Get the area of the convex area defined by the corners of the current area.

get_area_extent_for_subset(row_LR, col_LR, row_UL, col_UL)

Retrieves area_extent for a subdomain rows are counted from upper left to lower left columns are counted from upper left to upper right

Parameters:
row_LR : int
row of the lower right pixel
col_LR : int
col of the lower right pixel
row_UL : int
row of the upper left pixel
col_UL : int
col of the upper left pixel
Returns:
area_extent : list
Area extent as a list (LL_x, LL_y, UR_x, UR_y) of the subset
Author:

Ulrich Hamann

get_boundary_lonlats()

Returns Boundary objects

get_cartesian_coords(nprocs=None, data_slice=None, cache=False)

Retrieve cartesian coordinates of geometry definition

Parameters:
  • nprocs (int, optional) – Number of processor cores to be used. Defaults to the nprocs set when instantiating object
  • data_slice (slice object, optional) – Calculate only cartesian coordnates for the defined slice
  • cache (bool, optional) – Store result the result. Requires data_slice to be None
Returns:

cartesian_coords

Return type:

numpy array

get_lonlat(row, col)

Retrieve lon and lat of single pixel

Parameters:
  • row (int) –
  • col (int) –
Returns:

(lon, lat)

Return type:

tuple of floats

get_lonlats(data_slice=None, **kwargs)

Base method for lon lat retrieval with slicing

intersection(other)

Returns the corners of the intersection polygon of the current area with other.

Parameters:other (object) – Instance of subclass of BaseDefinition
Returns:(corner1, corner2, corner3, corner4)
Return type:tuple of points
overlap_rate(other)

Get how much the current area overlaps an other area.

Parameters:other (object) – Instance of subclass of BaseDefinition
Returns:overlap_rate
Return type:float
overlaps(other)

Tests if the current area overlaps the other area. This is based solely on the corners of areas, assuming the boundaries to be great circles.

Parameters:other (object) – Instance of subclass of BaseDefinition
Returns:overlaps
Return type:bool
class geometry.Boundary(side1, side2, side3, side4)

Container for geometry boundary. Labelling starts in upper left corner and proceeds clockwise

class geometry.CoordinateDefinition(lons, lats, nprocs=1)

Base class for geometry definitions defined by lons and lats only

class geometry.GridDefinition(lons, lats, nprocs=1)

Grid defined by lons and lats

Parameters:
  • lons (numpy array) –
  • lats (numpy array) –
  • nprocs (int, optional) – Number of processor cores to be used for calculations.
shape

tuple – Grid shape as (rows, cols)

size

int – Number of elements in grid

lons

object – Grid lons

lats

object – Grid lats

cartesian_coords

object – Grid cartesian coordinates

exception geometry.IncompatibleAreas

Error when the areas to combine are not compatible.

class geometry.StackedAreaDefinition(*definitions, **kwargs)

Definition based on muliple vertically stacked AreaDefinitions.

append(definition)

Append another definition to the area.

get_lonlats(nprocs=None, data_slice=None, cache=False, dtype=None)

Return lon and lat arrays of the area.

proj4_string

Returns projection definition as Proj.4 string

squeeze()

Generate a single AreaDefinition if possible.

class geometry.SwathDefinition(lons, lats, nprocs=1)

Swath defined by lons and lats

Parameters:
  • lons (numpy array) –
  • lats (numpy array) –
  • nprocs (int, optional) – Number of processor cores to be used for calculations.
shape

tuple – Swath shape

size

int – Number of elements in swath

ndims

int – Swath dimensions

lons

object – Swath lons

lats

object – Swath lats

cartesian_coords

object – Swath cartesian coordinates

geometry.combine_area_extents_vertical(area1, area2)

Combine the area extents of areas 1 and 2.

geometry.concatenate_area_defs(area1, area2, axis=0)

Append area2 to area1 and return the results

pyresample.image

Handles resampling of images with assigned geometry definitions

class image.ImageContainer(image_data, geo_def, fill_value=0, nprocs=1)

Holds image with geometry definition. Allows indexing with linesample arrays.

Parameters:
  • image_data (numpy array) – Image data
  • geo_def (object) – Geometry definition
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • nprocs (int, optional) – Number of processor cores to be used
image_data

numpy array – Image data

geo_def

object – Geometry definition

fill_value

int or None – Resample result fill value

nprocs

int – Number of processor cores to be used for geometry operations

get_array_from_linesample(row_indices, col_indices)

Samples from image based on index arrays.

Parameters:
  • row_indices (numpy array) – Row indices. Dimensions must match col_indices
  • col_indices (numpy array) – Col indices. Dimensions must match row_indices
Returns:

image_data – Resampled image data

Return type:

numpy_array

get_array_from_neighbour_info(*args, **kwargs)

Base method for resampling from preprocessed data.

resample(target_geo_def)

Base method for resampling

class image.ImageContainerNearest(image_data, geo_def, radius_of_influence, epsilon=0, fill_value=0, reduce_data=True, nprocs=1, segments=None)

Holds image with geometry definition. Allows nearest neighbour resampling to new geometry definition.

Parameters:
  • image_data (numpy array) – Image data
  • geo_def (object) – Geometry definition
  • radius_of_influence (float) – Cut off distance in meters
  • epsilon (float, optional) – Allowed uncertainty in meters. Increasing uncertainty reduces execution time
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • reduce_data (bool, optional) – Perform coarse data reduction before resampling in order to reduce execution time
  • nprocs (int, optional) – Number of processor cores to be used for geometry operations
  • segments (int or None) – Number of segments to use when resampling. If set to None an estimate will be calculated
image_data

numpy array – Image data

geo_def

object – Geometry definition

radius_of_influence

float – Cut off distance in meters

epsilon

float – Allowed uncertainty in meters

fill_value

int or None – Resample result fill value

reduce_data

bool – Perform coarse data reduction before resampling

nprocs

int – Number of processor cores to be used

segments

int or None – Number of segments to use when resampling

resample(target_geo_def)

Resamples image to area definition using nearest neighbour approach

Parameters:target_geo_def (object) – Target geometry definition
Returns:image_container – ImageContainerNearest object of resampled geometry
Return type:object
class image.ImageContainerQuick(image_data, geo_def, fill_value=0, nprocs=1, segments=None)

Holds image with area definition. ‘ Allows quick resampling within area.

Parameters:
  • image_data (numpy array) – Image data
  • geo_def (object) – Area definition as AreaDefinition object
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • nprocs (int, optional) – Number of processor cores to be used for geometry operations
  • segments (int or None) – Number of segments to use when resampling. If set to None an estimate will be calculated
image_data

numpy array – Image data

geo_def

object – Area definition as AreaDefinition object

fill_value

int or None – Resample result fill value If fill_value is None a masked array is returned with undetermined pixels masked

nprocs

int – Number of processor cores to be used

segments

int or None – Number of segments to use when resampling

resample(target_area_def)

Resamples image to area definition using nearest neighbour approach in projection coordinates.

Parameters:target_area_def (object) – Target area definition as AreaDefinition object
Returns:image_container – ImageContainerQuick object of resampled area
Return type:object

pyresample.grid

Resample image from one projection to another using nearest neighbour method in cartesian projection coordinate systems

grid.get_image_from_linesample(row_indices, col_indices, source_image, fill_value=0)

Samples from image based on index arrays.

Parameters:
  • row_indices (numpy array) – Row indices. Dimensions must match col_indices
  • col_indices (numpy array) – Col indices. Dimensions must match row_indices
  • source_image (numpy array) – Source image
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
Returns:

image_data – Resampled image

Return type:

numpy array

grid.get_image_from_lonlats(lons, lats, source_area_def, source_image_data, fill_value=0, nprocs=1)

Samples from image based on lon lat arrays using nearest neighbour method in cartesian projection coordinate systems.

Parameters:
  • lons (numpy array) – Lons. Dimensions must match lats
  • lats (numpy array) – Lats. Dimensions must match lons
  • source_area_def (object) – Source definition as AreaDefinition object
  • source_image_data (numpy array) – Source image data
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • nprocs (int, optional) – Number of processor cores to be used
Returns:

image_data – Resampled image data

Return type:

numpy array

grid.get_linesample(lons, lats, source_area_def, nprocs=1)

Returns index row and col arrays for resampling

Parameters:
  • lons (numpy array) – Lons. Dimensions must match lats
  • lats (numpy array) – Lats. Dimensions must match lons
  • source_area_def (object) – Source definition as AreaDefinition object
  • nprocs (int, optional) – Number of processor cores to be used
Returns:

(row_indices, col_indices) – Arrays for resampling area by array indexing

Return type:

tuple of numpy arrays

grid.get_resampled_image(target_area_def, source_area_def, source_image_data, fill_value=0, nprocs=1, segments=None)

Resamples image using nearest neighbour method in cartesian projection coordinate systems.

Parameters:
  • target_area_def (object) – Target definition as AreaDefinition object
  • source_area_def (object) – Source definition as AreaDefinition object
  • source_image_data (numpy array) – Source image data
  • fill_value ({int, None} optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • nprocs (int, optional) – Number of processor cores to be used
  • segments ({int, None} optional) – Number of segments to use when resampling. If set to None an estimate will be calculated.
Returns:

image_data – Resampled image data

Return type:

numpy array

pyresample.kd_tree

Handles reprojection of geolocated data. Several types of resampling are supported

kd_tree.get_neighbour_info(source_geo_def, target_geo_def, radius_of_influence, neighbours=8, epsilon=0, reduce_data=True, nprocs=1, segments=None)

Returns neighbour info

Parameters:
  • source_geo_def (object) – Geometry definition of source
  • target_geo_def (object) – Geometry definition of target
  • radius_of_influence (float) – Cut off distance in meters
  • neighbours (int, optional) – The number of neigbours to consider for each grid point
  • epsilon (float, optional) – Allowed uncertainty in meters. Increasing uncertainty reduces execution time
  • reduce_data (bool, optional) – Perform initial coarse reduction of source dataset in order to reduce execution time
  • nprocs (int, optional) – Number of processor cores to be used
  • segments (int or None) – Number of segments to use when resampling. If set to None an estimate will be calculated
Returns:

  • (valid_input_index, valid_output_index,
  • index_array, distance_array) (tuple of numpy arrays) – Neighbour resampling info

kd_tree.get_sample_from_neighbour_info(resample_type, output_shape, data, valid_input_index, valid_output_index, index_array, distance_array=None, weight_funcs=None, fill_value=0, with_uncert=False)

Resamples swath based on neighbour info

Parameters:
  • resample_type ({'nn', 'custom'}) – ‘nn’: Use nearest neighbour resampling ‘custom’: Resample based on weight_funcs
  • output_shape ((int, int)) – Shape of output as (rows, cols)
  • data (numpy array) – Source data
  • valid_input_index (numpy array) – valid_input_index from get_neighbour_info
  • valid_output_index (numpy array) – valid_output_index from get_neighbour_info
  • index_array (numpy array) – index_array from get_neighbour_info
  • distance_array (numpy array, optional) – distance_array from get_neighbour_info Not needed for ‘nn’ resample type
  • weight_funcs (list of function objects or function object, optional) – List of weight functions f(dist) to use for the weighting of each channel 1 to k. If only one channel is resampled weight_funcs is a single function object. Must be supplied when using ‘custom’ resample type
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
Returns:

result – Source data resampled to target geometry

Return type:

numpy array

kd_tree.resample_custom(source_geo_def, data, target_geo_def, radius_of_influence, weight_funcs, neighbours=8, epsilon=0, fill_value=0, reduce_data=True, nprocs=1, segments=None, with_uncert=False)

Resamples data using kd-tree custom radial weighting neighbour approach

Parameters:
  • source_geo_def (object) – Geometry definition of source
  • data (numpy array) – Array of single channel data points or (source_geo_def.shape, k) array of k channels of datapoints
  • target_geo_def (object) – Geometry definition of target
  • radius_of_influence (float) – Cut off distance in meters
  • weight_funcs (list of function objects or function object) – List of weight functions f(dist) to use for the weighting of each channel 1 to k. If only one channel is resampled weight_funcs is a single function object.
  • neighbours (int, optional) – The number of neigbours to consider for each grid point
  • epsilon (float, optional) – Allowed uncertainty in meters. Increasing uncertainty reduces execution time
  • fill_value ({int, None}, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • reduce_data (bool, optional) – Perform initial coarse reduction of source dataset in order to reduce execution time
  • nprocs (int, optional) – Number of processor cores to be used
  • segments ({int, None}) – Number of segments to use when resampling. If set to None an estimate will be calculated
Returns:

  • data (numpy array (default)) – Source data resampled to target geometry
  • data, stddev, counts (numpy array, numpy array, numpy array (if with_uncert == True)) – Source data resampled to target geometry. Weighted standard devaition for all pixels having more than one source value Counts of number of source values used in weighting per pixel

kd_tree.resample_gauss(source_geo_def, data, target_geo_def, radius_of_influence, sigmas, neighbours=8, epsilon=0, fill_value=0, reduce_data=True, nprocs=1, segments=None, with_uncert=False)

Resamples data using kd-tree gaussian weighting neighbour approach

Parameters:
  • source_geo_def (object) – Geometry definition of source
  • data (numpy array) – Array of single channel data points or (source_geo_def.shape, k) array of k channels of datapoints
  • target_geo_def (object) – Geometry definition of target
  • radius_of_influence (float) – Cut off distance in meters
  • sigmas (list of floats or float) – List of sigmas to use for the gauss weighting of each channel 1 to k, w_k = exp(-dist^2/sigma_k^2). If only one channel is resampled sigmas is a single float value.
  • neighbours (int, optional) – The number of neigbours to consider for each grid point
  • epsilon (float, optional) – Allowed uncertainty in meters. Increasing uncertainty reduces execution time
  • fill_value ({int, None}, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • reduce_data (bool, optional) – Perform initial coarse reduction of source dataset in order to reduce execution time
  • nprocs (int, optional) – Number of processor cores to be used
  • segments (int or None) – Number of segments to use when resampling. If set to None an estimate will be calculated
  • with_uncert (bool, optional) – Calculate uncertainty estimates
Returns:

  • data (numpy array (default)) – Source data resampled to target geometry
  • data, stddev, counts (numpy array, numpy array, numpy array (if with_uncert == True)) – Source data resampled to target geometry. Weighted standard devaition for all pixels having more than one source value Counts of number of source values used in weighting per pixel

kd_tree.resample_nearest(source_geo_def, data, target_geo_def, radius_of_influence, epsilon=0, fill_value=0, reduce_data=True, nprocs=1, segments=None)

Resamples data using kd-tree nearest neighbour approach

Parameters:
  • source_geo_def (object) – Geometry definition of source
  • data (numpy array) – 1d array of single channel data points or (source_size, k) array of k channels of datapoints
  • target_geo_def (object) – Geometry definition of target
  • radius_of_influence (float) – Cut off distance in meters
  • epsilon (float, optional) – Allowed uncertainty in meters. Increasing uncertainty reduces execution time
  • fill_value (int or None, optional) – Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
  • reduce_data (bool, optional) – Perform initial coarse reduction of source dataset in order to reduce execution time
  • nprocs (int, optional) – Number of processor cores to be used
  • segments (int or None) – Number of segments to use when resampling. If set to None an estimate will be calculated
Returns:

data – Source data resampled to target geometry

Return type:

numpy array

kd_tree.which_kdtree()

Returns the name of the kdtree used for resampling

pyresample.bilinear

Code for resampling using bilinear algorithm for irregular grids.

The algorithm is taken from

http://www.ahinson.com/algorithms_general/Sections/InterpolationRegression/InterpolationIrregularBilinear.pdf

bilinear.get_bil_info(source_geo_def, target_area_def, radius=50000.0, neighbours=32, nprocs=1, masked=False, reduce_data=True, segments=None, epsilon=0)

Calculate information needed for bilinear resampling.

source_geo_def : object
Geometry definition of source data
target_area_def : object
Geometry definition of target area
radius : float, optional
Cut-off distance in meters
neighbours : int, optional
Number of neighbours to consider for each grid point when searching the closest corner points
nprocs : int, optional
Number of processor cores to be used for getting neighbour info
masked : bool, optional
If true, return masked arrays, else return np.nan values for invalid points (default)
reduce_data : bool, optional
Perform initial coarse reduction of source dataset in order to reduce execution time
segments : int or None
Number of segments to use when resampling. If set to None an estimate will be calculated
epsilon : float, optional
Allowed uncertainty in meters. Increasing uncertainty reduces execution time
Returns:
  • t__ (numpy array) – Vertical fractional distances from corner to the new points
  • s__ (numpy array) – Horizontal fractional distances from corner to the new points
  • input_idxs (numpy array) – Valid indices in the input data
  • idx_arr (numpy array) – Mapping array from valid source points to target points
bilinear.get_sample_from_bil_info(data, t__, s__, input_idxs, idx_arr, output_shape=None)

Resample data using bilinear interpolation.

Parameters:
  • data (numpy array) – 1d array to be resampled
  • t (numpy array) – Vertical fractional distances from corner to the new points
  • s (numpy array) – Horizontal fractional distances from corner to the new points
  • input_idxs (numpy array) – Valid indices in the input data
  • idx_arr (numpy array) – Mapping array from valid source points to target points
  • output_shape (tuple, optional) – Tuple of (y, x) dimension for the target projection. If None (default), do not reshape data.
Returns:

result – Source data resampled to target geometry

Return type:

numpy array

bilinear.resample_bilinear(data, source_geo_def, target_area_def, radius=50000.0, neighbours=32, nprocs=1, fill_value=0, reduce_data=True, segments=None, epsilon=0)

Resample using bilinear interpolation.

data : numpy array
Array of single channel data points or (source_geo_def.shape, k) array of k channels of datapoints
source_geo_def : object
Geometry definition of source data
target_area_def : object
Geometry definition of target area
radius : float, optional
Cut-off distance in meters
neighbours : int, optional
Number of neighbours to consider for each grid point when searching the closest corner points
nprocs : int, optional
Number of processor cores to be used for getting neighbour info
fill_value : {int, None}, optional
Set undetermined pixels to this value. If fill_value is None a masked array is returned with undetermined pixels masked
reduce_data : bool, optional
Perform initial coarse reduction of source dataset in order to reduce execution time
segments : int or None
Number of segments to use when resampling. If set to None an estimate will be calculated
epsilon : float, optional
Allowed uncertainty in meters. Increasing uncertainty reduces execution time
Returns:data – Source data resampled to target geometry
Return type:numpy array

pyresample.utils

Utility functions for pyresample

exception utils.AreaNotFound

Exception raised when specified are is no found in file

utils.fwhm2sigma(fwhm)

Calculate sigma for gauss function from FWHM (3 dB level)

Parameters:fwhm (float) – FWHM of gauss function (3 dB level of beam footprint)
Returns:sigma – sigma for use in resampling gauss function
Return type:float
utils.generate_nearest_neighbour_linesample_arrays(source_area_def, target_area_def, radius_of_influence, nprocs=1)

Generate linesample arrays for nearest neighbour grid resampling

Parameters:
  • source_area_def (object) – Source area definition as AreaDefinition object
  • target_area_def (object) – Target area definition as AreaDefinition object
  • radius_of_influence (float) – Cut off distance in meters
  • nprocs (int, optional) – Number of processor cores to be used
Returns:

(row_indices, col_indices)

Return type:

tuple of numpy arrays

utils.generate_quick_linesample_arrays(source_area_def, target_area_def, nprocs=1)

Generate linesample arrays for quick grid resampling

Parameters:
  • source_area_def (object) – Source area definition as AreaDefinition object
  • target_area_def (object) – Target area definition as AreaDefinition object
  • nprocs (int, optional) – Number of processor cores to be used
Returns:

(row_indices, col_indices)

Return type:

tuple of numpy arrays

utils.get_area_def(area_id, area_name, proj_id, proj4_args, x_size, y_size, area_extent)

Construct AreaDefinition object from arguments

Parameters:
  • area_id (str) – ID of area
  • proj_id (str) – ID of projection
  • area_name (str) – Description of area
  • proj4_args (list or str) – Proj4 arguments as list of arguments or string
  • x_size (int) – Number of pixel in x dimension
  • y_size (int) – Number of pixel in y dimension
  • area_extent (list) – Area extent as a list of ints (LL_x, LL_y, UR_x, UR_y)
Returns:

area_def – AreaDefinition object

Return type:

object

utils.load_area(area_file_name, *regions)

Load area(s) from area file

Parameters:
  • area_file_name (str) – Path to area definition file
  • regions (str argument list) – Regions to parse. If no regions are specified all regions in the file are returned
Returns:

area_defs – If one area name is specified a single AreaDefinition object is returned If several area names are specified a list of AreaDefinition objects is returned

Return type:

object or list

Raises:

AreaNotFound: – If a specified area name is not found

utils.parse_area_file(area_file_name, *regions)

Parse area information from area file

Parameters:
  • area_file_name (str) – Path to area definition file
  • regions (str argument list) – Regions to parse. If no regions are specified all regions in the file are returned
Returns:

area_defs – List of AreaDefinition objects

Return type:

list

Raises:

AreaNotFound: – If a specified area is not found

utils.proj4_str_to_dict(proj4_str)

Convert PROJ.4 compatible string definition to dict

Note: Key only parameters will be assigned a value of True.

utils.wrap_longitudes(lons)

Wrap longitudes to the [-180:+180[ validity range (preserves dtype)

Parameters:lons (numpy array) – Longitudes in degrees
Returns:lons – Longitudes wrapped into [-180:+180[ validity range
Return type:numpy array

pyresample.data_reduce

Reduce data sets based on geographical information

data_reduce.get_valid_index_from_cartesian_grid(cart_grid, lons, lats, radius_of_influence)

Calculates relevant data indices using coarse data reduction of swath data by comparison with cartesian grid

Parameters:
  • chart_grid (numpy array) – Grid of area cartesian coordinates
  • lons (numpy array) – Swath lons
  • lats (numpy array) – Swath lats
  • data (numpy array) – Swath data
  • radius_of_influence (float) – Cut off distance in meters
Returns:

valid_index – Boolean array of same size as lons and lats indicating relevant indices

Return type:

numpy array

data_reduce.get_valid_index_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, radius_of_influence)

Find relevant indices from grid boundaries using the winding number theorem

data_reduce.get_valid_index_from_lonlat_grid(grid_lons, grid_lats, lons, lats, radius_of_influence)

Calculates relevant data indices using coarse data reduction of swath data by comparison with lon lat grid

Parameters:
  • chart_grid (numpy array) – Grid of area cartesian coordinates
  • lons (numpy array) – Swath lons
  • lats (numpy array) – Swath lats
  • data (numpy array) – Swath data
  • radius_of_influence (float) – Cut off distance in meters
Returns:

valid_index – Boolean array of same size as lon and lat indicating relevant indices

Return type:

numpy array

data_reduce.swath_from_cartesian_grid(cart_grid, lons, lats, data, radius_of_influence)

Makes coarse data reduction of swath data by comparison with cartesian grid

Parameters:
  • chart_grid (numpy array) – Grid of area cartesian coordinates
  • lons (numpy array) – Swath lons
  • lats (numpy array) – Swath lats
  • data (numpy array) – Swath data
  • radius_of_influence (float) – Cut off distance in meters
Returns:

(lons, lats, data) – Reduced swath data and coordinate set

Return type:

list of numpy arrays

data_reduce.swath_from_lonlat_boundaries(boundary_lons, boundary_lats, lons, lats, data, radius_of_influence)

Makes coarse data reduction of swath data by comparison with lon lat boundary

Parameters:
  • boundary_lons (numpy array) – Grid of area lons
  • boundary_lats (numpy array) – Grid of area lats
  • lons (numpy array) – Swath lons
  • lats (numpy array) – Swath lats
  • data (numpy array) – Swath data
  • radius_of_influence (float) – Cut off distance in meters
Returns:

(lons, lats, data) – Reduced swath data and coordinate set

Return type:

list of numpy arrays

data_reduce.swath_from_lonlat_grid(grid_lons, grid_lats, lons, lats, data, radius_of_influence)

Makes coarse data reduction of swath data by comparison with lon lat grid

Parameters:
  • grid_lons (numpy array) – Grid of area lons
  • grid_lats (numpy array) – Grid of area lats
  • lons (numpy array) – Swath lons
  • lats (numpy array) – Swath lats
  • data (numpy array) – Swath data
  • radius_of_influence (float) – Cut off distance in meters
Returns:

(lons, lats, data) – Reduced swath data and coordinate set

Return type:

list of numpy arrays

pyresample.plot

plot.area_def2basemap(area_def, **kwargs)

Get Basemap object from AreaDefinition

Parameters:
  • area_def (object) – geometry.AreaDefinition object
  • **kwargs (Keyword arguments) – Additional initialization arguments for Basemap
Returns:

bmap

Return type:

Basemap object

plot.ellps2axis(ellps_name)

Get semi-major and semi-minor axis from ellipsis definition

Parameters:ellps_name (str) – Standard name of ellipsis
Returns:(a, b)
Return type:semi-major and semi-minor axis
plot.save_quicklook(filename, area_def, data, vmin=None, vmax=None, label='Variable (units)', num_meridians=45, num_parallels=10, coast_res='c', backend='AGG', cmap='jet')

Display default quicklook plot

Parameters:
  • filename (str) – path to output file
  • area_def (object) – geometry.AreaDefinition object
  • data (numpy array | numpy masked array) – 2D array matching area_def. Use masked array for transparent values
  • vmin (float, optional) – Min value for luminescence scaling
  • vmax (float, optional) – Max value for luminescence scaling
  • label (str, optional) – Label for data
  • num_meridians (int, optional) – Number of meridians to plot on the globe
  • num_parallels (int, optional) – Number of parallels to plot on the globe
  • coast_res ({'c', 'l', 'i', 'h', 'f'}, optional) – Resolution of coastlines
  • backend (str, optional) – matplotlib backend to use’
plot.show_quicklook(area_def, data, vmin=None, vmax=None, label='Variable (units)', num_meridians=45, num_parallels=10, coast_res='c', cmap='jet')

Display default quicklook plot

Parameters:
  • area_def (object) – geometry.AreaDefinition object
  • data (numpy array | numpy masked array) – 2D array matching area_def. Use masked array for transparent values
  • vmin (float, optional) – Min value for luminescence scaling
  • vmax (float, optional) – Max value for luminescence scaling
  • label (str, optional) – Label for data
  • num_meridians (int, optional) – Number of meridians to plot on the globe
  • num_parallels (int, optional) – Number of parallels to plot on the globe
  • coast_res ({'c', 'l', 'i', 'h', 'f'}, optional) – Resolution of coastlines
Returns:

bmap

Return type:

Basemap object

pyresample.ewa

Code for resampling using the Elliptical Weighted Averaging (EWA) algorithm.

The logic and original code for this algorithm were translated from the software package “MODIS Swath 2 Grid Toolbox” or “ms2gt” created by the NASA National Snow & Ice Data Center (NSIDC):

Since the project has slowed down, Terry Haran has maintained the package and made updates available:

The ms2gt C executables “ll2cr” and “fornav” were rewritten for the Polar2Grid software package created by the Space Science Engineering Center (SSEC)/Cooperative Institute for Meteorological Satellite Studies. They were rewritten as a combination of C++ and Cython to make them more python friendly by David Hoese and were then copied and modified here in pyresample. The rewrite of “ll2cr” also included an important switch from using the “mapx” library to using the more popular and capable pyproj (PROJ.4) library.

The EWA algorithm consists of two parts “ll2cr” and “fornav” and are described below.

ll2cr

The “ll2cr” process is the first step in the EWA algorithm. It stands for “latitude/longitude to column/row”. Its main purpose is to convert input longitude and latitude coordinates to column and row coordinates of the destination grid. These coordinates are then used in the next step “fornav”.

fornav

The “fornav” or “Forward Navigation” step of the EWA algorithm is where the actual Elliptical Weighted Averaging algorithm is run. The algorithm maps input swath pixels to output grid pixels by averaging multiple input pixels based on an elliptical region and other coefficients, some of which are determined at run time.

For more information on these steps see the documentation for the corresponding modules.

ewa.fornav(cols, rows, area_def, data_in, rows_per_scan=None, fill=None, out=None, weight_count=10000, weight_min=0.01, weight_distance_max=1.0, weight_delta_max=10.0, weight_sum_min=-1.0, maximum_weight_mode=False)

Remap data in to output grid using elliptical weighted averaging.

This algorithm works under the assumption that the data is observed one scan line at a time. However, good results can still be achieved for non-scan based data is provided if rows_per_scan is set to the number of rows in the entire swath or by setting it to None.

Parameters:
  • cols (numpy array) – Column location for each input swath pixel (from ll2cr)
  • rows (numpy array) – Row location for each input swath pixel (from ll2cr)
  • area_def (AreaDefinition) – Grid definition to be mapped to
  • data_in (numpy array or tuple of numpy arrays) – Swath data to be remapped to output grid
  • rows_per_scan (int or None, optional) – Number of data rows for every observed scanline. If None then the entire swath is treated as one large scanline.
  • fill (float/int or None, optional) – If data_in is made of numpy arrays then this represents the fill value used to mark invalid data pixels. This value will also be used in the output array(s). If None, then np.nan will be used for float arrays and -999 will be used for integer arrays.
  • out (numpy array or tuple of numpy arrays, optional) – Specify a numpy array to be written to for each input array. This can be used as an optimization by providing np.memmap arrays or other array-like objects.
  • weight_count (int, optional) – number of elements to create in the gaussian weight table. Default is 10000. Must be at least 2
  • weight_min (float, optional) – the minimum value to store in the last position of the weight table. Default is 0.01, which, with a weight_distance_max of 1.0 produces a weight of 0.01 at a grid cell distance of 1.0. Must be greater than 0.
  • weight_distance_max (float, optional) – distance in grid cell units at which to apply a weight of weight_min. Default is 1.0. Must be greater than 0.
  • weight_delta_max (float, optional) – maximum distance in grid cells in each grid dimension over which to distribute a single swath cell. Default is 10.0.
  • weight_sum_min (float, optional) – minimum weight sum value. Cells whose weight sums are less than weight_sum_min are set to the grid fill value. Default is EPSILON.
  • maximum_weight_mode (bool, optional) – If False (default), a weighted average of all swath cells that map to a particular grid cell is used. If True, the swath cell having the maximum weight of all swath cells that map to a particular grid cell is used. This option should be used for coded/category data, i.e. snow cover.
Returns:

(valid grid points, output arrays) – The valid_grid_points tuple holds the number of output grid pixels that were written with valid data. The second element in the tuple is a tuple of output grid numpy arrays for each input array. If there was only one input array provided then the returned tuple is simply the singe points integer and single output grid array.

Return type:

tuple of integer tuples and numpy array tuples

ewa.ll2cr(swath_def, area_def, fill=<Mock object>, copy=True)

Map input swath pixels to output grid column and rows.

Parameters:
  • swath_def (SwathDefinition) – Navigation definition for swath data to remap
  • area_def (AreaDefinition) – Grid definition to be mapped to
  • fill (float, optional) – Fill value used in longitude and latitude arrays
  • copy (bool, optional) – Create a copy of the longitude and latitude arrays (default: True)
Returns:

(swath_points_in_grid, cols, rows) – Number of points from the input swath overlapping the destination area and the column and row arrays to pass to fornav.

Return type:

tuple of integer, numpy array, numpy array

Note

ll2cr uses the pyproj library which is limited to 64-bit float navigation arrays in order to not do additional copying or casting of data types.