Python 重采样遥感数据 Pyresample (一)

中国地质大学（北京）工学硕士

遥感影像的重采样似乎从来都不是什么大问题，不管是以 ENVI 为代表的商业软件，还是以 GDAL 为代表的开源开发库，甚至我们以为的“业余”的 Photoshop 都提供了完善的、功能强大的影像重采样功能。直到，当我们要去处理诸如 Sentinel-5, Ecostress, Modis等为代表的宽幅扫描遥感卫星的 L1 级甚至更低级别数据产品时，发现先前的工具似乎不太好用了，甚至根本不可用。因为这些数据的像素尺寸是变化的，我们在使用中更喜欢规则化的栅格，这样就需要重采样。虽然 GDAL 也提供了将非规则化的数据点重采样为规则化栅格的工具，但在面对数千x数千的数据点数组尺寸时，对内存、CPU的消耗简直让人沮丧。在经历了一段失败的编程经历之后，发现Pyresample 为我们提供了高效、完备的解决方法。它使用了先进的 kd-tree 算法，及以 Xarray、Dask等为代表最新的并行化计算技术，使得其在应对大尺度的、非规则的遥感影像重采样时显得游刃有余。

Pyresample 简介

Pyresample 是一个用于重采样地理空间影像数据的 python 包。它为 SatPy 库提供了重采样的主要方法，同时它也可以单独使用。重采样或重投影是将地理定位的数据点栅格化为新的目标地理投影和地理区域的过程。

Pyresample 可以处理栅格数据和地理定位的刈幅数据。为了描述这些数据， Pyresample 使用不同的“geometry” 对象，包括 AreaDefinition 和 SwathDefinition 类。

Pyresample 提供多种重采样算法:

Nearest Neighbor
Elliptical Weighted Average (EWA)
Bilinear
Bucket resampling (count hits per bin, averaging, ratios)

对于最近邻法和双线性插值法， pyresample 使用一种 kd-tree 方法，由 pykdtree 库实现的快速 KDTree 算法提供处理引擎。Pyresample 使用 numpy 数组和 numpy 掩膜标注数组；同时它也在另外单独的 Resampler 类中提供XArray 对象的接口 (包括对 dask 数组的支持)，并正处于开发中。Pyresample 提供的工具函数使用 Cartopy 可视化数据非常容易。

Versionchanged: Changed in version 1.15: 放弃对 Python 2 和 Python <3.4 的支持。

安装 Pyresample

Pyresample 依赖 pyproj, numpy(>= 1.10), pyyaml, configobj
和 pykdtree (>= 1.1.1)。

为了使用 pyresample 的绘图功能， Cartopy 和
matplotlib (>= 1.0) 必须安装。这些包并不是使用 pyresample 其他功能的先决条件。

对于需要 dask 和 xarray 支持的功能，这些包必须安装。一些功能，比如将 rasterio 对象转换为 pyresample 对象，要求 rasterio 或其它一些包必须安装。老版本的并行处理接口
(Proj_MP) 使用了 scipy 包的 KDTree 算法实现，当 pyresample 接口中的
nprocs 关键字参数值大于 1 时，会使用该并行处理功能。在可能的情况下，更新的 xarray/dask 接口是更优推荐项。

包测试

测试 pyresample 需要安装所有可选包，包括rasterio, dask, xarray, cartopy, pillow, 和 matplotlib；否则，一些测试可能会失败。

从源包进行测试:

tar -zxf pyresample-<version>.tar.gz
cd pyresample-<version>
pytest pyresample/test/

如果通过了所有测试，意味着 pyresample 所有函数的功能在系统上得到验证。

安装

Pyresample 可通过 pip 从 PyPI 上安装:

pip install pyresample

Pyresample 也可以使用 conda 通过 conda-forge 通道安装:

conda install -c conda-forge pyresample

或者直接从源包安装:

tar -zxvf pyresample-<version>.tar.gz
cd pyresample-<version>
pip install .

为将 pyresample 安装在 “development” 模式下（这样的话，源文件的改变会直接反映在你的 python 环境中），请运行如下命令:

pip install -e .

pykdtree 和 numexpr

Pyresample 使用了 pykdtree 包，它可以编译为支持多线程。如果它被编译为支持多线程，可以通过环境变量 OMP_NUM_THREADS 控制线程数量。
请参考 pykdtree 仓库了解更多信息。

从 v1.0.0 版本起， numexpr 会被用来优化处理瓶颈，如果它在你的系统中可用的话。

Geometry 定义

pyresample.geometry 模块包含使用网格点或像素描述不同地理区域的类。一些类将地理区域表达为等间距（尺寸）的像素，其它的类处理由非一致像素描述的地理区域。描述一个地理区域的最好对象取决于用例和与之有关的信息。pyresample 中可用的类描述如下。

注意，传递给 pyresample.geometry 的类的所有经度和纬度必须以度为单位。另外，经度必须在[-180;+180] 的有效范围内。

Versionchanged: Changed in version 1.8.0: 为了提高性能， Geometry 对象不再检查经度和纬度坐标的有效性。经度被期望介于 -180 和 180 度，纬度介于 -90 到 90 度。这对于所有的初始化阶段需要经纬度数组的地理几何定义同样适用。

一个 AreaDefinition , 或 area , 是 pyresample 描述一个均匀间隔的像素组成的地理区域的主要方式，它也是唯一的使用地理坐标系的地理区域对象。 Areas 使用 PROJ.4 方法描述其坐标系 (CRS)。如果一个区的投影没有用经度/纬度坐标的形式描述，那么它会采用以米为单位的 X/Y 坐标形式。参见 PROJ.4 文档了解有关投影和坐标系的更多信息。

下列参数用来初始化一个 AreaDefinition 对象:

area_id : ID
description : 描述
proj_id : 投影 ID (已废弃)
projection : Proj4 参数（字典或字符串形式）
width : 栅格列数
height : 栅格行数
area_extent : (lower_left_x, lower_left_y, upper_right_x, upper_right_y)

其中

lower_left_x : 左下角像素的左下角 x 坐标
lower_left_y : 左下角像素的左下角 y 坐标
upper_right_x : 右上角像素的右上角 x 坐标
upper_right_y : 右上角像素的右上角 y 坐标

示例:

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

你也可以使用 PROJ.4 字符串指定投影坐标系。

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

或者一个 EPSG code :

>>> projection = '+init=EPSG:3409'  # 对 2.0+ 版本pyproj，使用 'EPSG:3409'
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)

NOTE : 对于 2.0+ 版本的 pyproj 请使用 'EPSG:XXXX' 语法， '+init=EPSG:XXXX' 不再支持。

创建一个 AreaDefinition 对象可能会很复杂，如果你并不知道要描述的区域的所有信息。
Pyresample 提供了用于创建、保存区域的多个工具。参见
Geometry Utilities 文档了解更多信息。

如果每个像素的经度和纬度值已知,那么即可跳过创建 AreaDefinition 的复杂性，直接创建
GridDefinition 对象即可。
注意，尽管 grid 定义更简单，但它几乎为所有的操作带来更多的内存和计算成本。传递给 GridDefinition 的经度和纬度数组应是均匀间隔的；否则，你应该使用 SwathDefinition (见下)。

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

一个刈幅通过每个像素的经度和纬度坐标来定义。坐标代表像素的中心位置。刈幅不假定像素尺寸和间隔的均一性。这意味着操作它将花费更长的时间，当然也很精确。

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

两个刈幅可以拼接在一起，如果它们的列维数一致。

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

所有的几何定义对象都提供了访问像素经度和纬度坐标值的接口。
get_lonlats() 方法用来得到像素的经纬度坐标数据，它会自动执行可能需要的投影变换。

AreaDefinition 用成员变量 projection_x_coords 和 projection_y_coords 来提供整套坐标值的访问接口。注意，对于地理坐标系 (+proj=latlong)，坐标值以度为单位，即 projection_x_coords 会是一个经度向量， projection_y_coords 是一个纬度向量。

Versionchanged: Changed in version 1.5.1: 将 proj_x_coords 更名为 projection_x_coords,将 proj_y_coords
更名为 projection_y_coords。

得到经度和纬度数组:

>>> area_id = 'ease_sh'
>>> description = 'Antarctic EASE grid'
>>> proj_id = 'ease_sh'
>>> projection = '+proj=laea +lat_0=-90 +lon_0=0 +a=6371228.0 +units=m'
>>> width = 425
>>> height = 425
>>> area_extent = (-5326849.0625,-5326849.0625,5326849.0625,5326849.0625)
>>> area_def = AreaDefinition(area_id, description, proj_id, projection,
...                           width, height, area_extent)
>>> lons, lats = area_def.get_lonlats()

得到地心 X, Y, Z 坐标:

>>> area_def = AreaDefinition(area_id, description, proj_id, projection,