使用 Python 将栅格从 CRS 重新投影到另一个的正确方法是什么？答案

【问题标题】：What is the correct way to reproject a raster from a CRS to another using Python?使用 Python 将栅格从 CRS 重新投影到另一个的正确方法是什么？
【发布时间】：2021-09-08 06:28:22
【问题描述】：

我有一个使用“epsg:32632”作为 CRS 的土地覆盖数据栅格（特别是 https://finder.creodias.eu 中的 /eodata/auxdata/S2GLC/2017/S2GLC_T32TMS_2017）。我想在“epsg:21781”上重新投影这个光栅。这是我用 xarray 打开光栅时的样子。

fn = 'data/S2GLC_T32TMS_2017/S2GLC_T32TMS_2017.tif'
da = xr.open_rasterio(fn).sel(band=1, drop=True)
da
<xarray.DataArray (y: 10980, x: 10980)>
[120560400 values with dtype=uint8]
Coordinates:
  * y        (y) float64 5.2e+06 5.2e+06 5.2e+06 ... 5.09e+06 5.09e+06 5.09e+06
  * x        (x) float64 4e+05 4e+05 4e+05 ... 5.097e+05 5.097e+05 5.098e+05
Attributes:
    transform:      (10.0, 0.0, 399960.0, 0.0, -10.0, 5200020.0)
    crs:            +init=epsg:32632
    res:            (10.0, 10.0)
    is_tiled:       0
    nodatavals:     (nan,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area
    INTERLEAVE:     BAND

我通常的工作流程是转换所有点坐标，创建目标网格并使用最近邻进行插值。看起来像这样的东西：

import numpy as np
import xarray as xr
import pyproj
from scipy.interpolate import griddata 

y = da.y.values
x = da.x.values
xx, yy = np.meshgrid(x,y)
# (n,2) point coordinates in the original CRS
src_coords = np.column_stack([xx.flatten(), yy.flatten()]) 

transformer = pyproj.transformer.Transformer.from_crs('epsg:32632', 'epsg:21781')
xx, yy = transformer.transform(src_coords[:,0], src_coords[:,1])
# (n,2) point coordinates in the destination CRS, which are not on a regular grid
dst_coords = np.column_stack([xx.flatten(), yy.flatten()]) 


# I define my destination **regular** grid coordinates
x = np.linspace(620005,719995,10)
y = np.linspace(199995,100005,10)
xx, yy = np.meshgrid(x,y)
dst_grid = np.column_stack([xx.flatten(), yy.flatten()])

# I interpolate onto the grid
reprojected_array = griddata(
    src_coords, da.values.flatten(), dst_coords, method='nearest'
).reshape(dst_shape)

虽然这种方法相当透明并且（显然）没有错误，但在处理数十亿个点时可能需要很长时间。最近，我发现了rasterio 的reproject 函数，我被它的速度所震撼。这就是我实现它的方式：

source = da.values
destination = np.zeros(dst_shape, np.int16)

res, aff = reproject(
    source,
    destination,
    src_transform=src_transform, # affine transformation from original data
    src_crs=src_crs,
    dst_transform=dst_transform, # affine transformation that corresponds to the grid defined in the other approach
    dst_crs=dst_crs,
    resampling=Resampling.nearest) # using nearest neighbors just like with scope's griddata

当然，我想比较结果，期望它们相同，但实际上并非如此，如图所示。

分辨率为 10 米，因此差异不大，但在与“epsg:21781”坐标中的精确卫星数据仔细比较后，看起来旧方法产生了更好的结果。

所以我的问题是：

为什么这些结果不同？
一种方法比另一种更好吗？是否存在应该选择其中一种的特定条件？

【问题讨论】：

标签： python geospatial gdal coordinate-transformation rasterio

【解决方案1】：

Griddata 在 Euclidean 距离内找到最近的点，在你给它的任何地图投影上。因此，来自管道的最近邻居像
4326 个数据点 --> 重新投影 --> 最接近欧几里得网格数据
查询点

依赖于“重新项目”。你能试试+proj=sinu +lon_0=middle lon 对于数据和查询？

人们真正想要的是一个具有大圆距离的最近邻引擎，不是欧几里得距离。对于小网格或赤道附近的网格，差异可能微不足道，但在芬兰则更少——cos 61° / cos 60° 约为 97 %。

TL;DR
是pyproj.transformer.Transformer.from_crs('epsg:32632', 'epsg:21781') “正确的” ？不知道。我没有看到测试套件，还有几个问题：
warp.reproject() generates the wrong result
roundtrip test \

“最近邻”在数据点之间定义不明确/敏感，例如沿线 x 或 y = int + 0.5 在 int 网格上。这很容易用 KDTree 进行测试。

xarray 使常规（笛卡尔）网格变得容易，但 afaik 不这样做 curvilinear (2d) grids.

【讨论】：