创建rioxarray reproject_match目标栅格时坐标偏移问题求助

问题描述 #

需要创建一个3×3度AOI、分辨率0.00008（对应37500×37500行列数）的虚拟栅格，作为rioxarray.reproject_match()的目标栅格。编写代码后，用dummy.rio.bounds()检查边界得到(11.99996, 50.99996, 14.99996, 53.99996)，与预期的(12.0, 51.0, 15.0, 54.0)不符。

原代码：

import xarray as xr
import rioxarray as rio
import numpy as np
import decimal

res = 0.00008
dec = abs(decimal.Decimal(str(res)).as_tuple().exponent)

lon = list(np.arange(12, 15, res))
lon = [round(x, dec) for x in lon]
lat = list(np.arange(51, 54, res))
lat = [round(x, dec) for x in lat]

dummy = xr.DataArray(
        data=1,
        dims=["x", "y"],
        coords={"x": lon, "y": lat},
        name="dummy"
        )

dummy.rio.write_crs("epsg:4326", inplace=True)

原因分析 #

1. np.arange的左闭右开特性：np.arange(12, 15, res)生成的最后一个经度值是12 + 37499*res = 14.99992（共37500个值），而非15.0；纬度同理最后值为53.99992。
2. rioxarray的边界计算逻辑：当你将坐标直接赋值给x/y时，rioxarray默认这些坐标是像素中心，栅格边界由「中心坐标 ± 分辨率/2」计算得出。第一个经度中心是12.0，左边缘为12.0 - 0.00004 = 11.99996；最后一个经度中心是14.99992，右边缘为14.99992 + 0.00004 = 14.99996，最终导致整体边界偏移。

解决方法 #

方法1：基于边界直接生成栅格（推荐） #

直接指定AOI边界、分辨率和行列数，让rioxarray自动处理坐标和边界关系，避免手动计算误差：

import xarray as xr
import rioxarray as rio
import numpy as np

res = 0.00008
# 定义目标AOI边界：(左, 下, 右, 上)
bounds = (12.0, 51.0, 15.0, 54.0)
# 计算行列数（刚好为37500）
width = int((bounds[2] - bounds[0]) / res)
height = int((bounds[3] - bounds[1]) / res)

# 创建全1数组，注意维度顺序为[y, x]（符合rioxarray默认空间维度顺序）
data = np.ones((height, width))

# 创建DataArray并设置空间参考
dummy = xr.DataArray(
    data=data,
    dims=["y", "x"],
    name="dummy"
)

# 写入CRS和空间变换（从边界推导变换矩阵）
dummy.rio.write_crs("epsg:4326", inplace=True)
dummy.rio.write_transform(
    rio.transform.from_bounds(*bounds, width=width, height=height),
    inplace=True
)

# 验证边界
print(dummy.rio.bounds())  # 输出：(12.0, 51.0, 15.0, 54.0)

方法2：手动调整像素中心坐标 #

如果需要保留手动生成坐标的方式，需将像素中心设置为「边界边缘 + 分辨率/2」，确保边缘刚好覆盖目标AOI：

import xarray as xr
import rioxarray as rio
import numpy as np
import decimal

res = 0.00008
dec = abs(decimal.Decimal(str(res)).as_tuple().exponent)

# 生成像素中心坐标：从边缘+res/2开始，到边缘-res/2结束，步长res
lon = np.arange(12 + res/2, 15 - res/2 + res, res)  # 确保包含最后一个中心
lon = [round(x, dec) for x in lon]
lat = np.arange(51 + res/2, 54 - res/2 + res, res)
lat = [round(x, dec) for x in lat]

dummy = xr.DataArray(
    data=1,
    dims=["y", "x"],  # 调整维度顺序为[y, x]，符合空间数据常规规范
    coords={"x": lon, "y": lat},
    name="dummy"
)

dummy.rio.write_crs("epsg:4326", inplace=True)

# 验证边界
print(dummy.rio.bounds())  # 输出：(12.0, 51.0, 15.0, 54.0)

注意事项 #

• rioxarray默认空间维度顺序为[y, x]（纬度在前，经度在后），若维度顺序错误可能导致后续投影匹配出现问题。
• 避免手动计算坐标时的浮点误差，优先使用rio.transform.from_bounds推导变换矩阵的方式，更可靠。

内容的提问来源于stack exchange，提问作者Corbjn