使用GDAL Warp对齐栅格时输出尺寸偏移问题排查

我尝试在QGIS中使用GDAL Warp工具将一个栅格与参考栅格对齐，待方法验证成熟后，我会通过Python GDAL wrapper将该流程封装到自动化脚本中。我查阅过栅格对齐相关的技术讨论帖，其中提到对齐操作需使用参考影像的范围坐标、像元分辨率，并添加-tap参数强制待对齐影像与参考栅格对齐，我执行了如下命令：

gdalwarp -s_srs EPSG:32630 -t_srs EPSG:32630 -tr 0.25 0.25 -r near -te 511711.9928 6328998.9708 514192.6233 6331359.4866 -of GTiff -tap <path-input> <path-output>

但输出结果仍存在轻微偏移：

参考影像的尺寸为9923×9442，而上述命令输出影像的尺寸为9924×9443。我需要二者尺寸、像元分辨率完全一致，能够完美重叠，以便后续转换为numpy数组开展处理。
我知晓QGIS中自带Align Raster栅格对齐工具，但我需要基于GDAL实现对齐逻辑，以便适配Python脚本开发需求。
以下为两份原始栅格的属性信息，本次对齐的目标是匹配较低分辨率栅格的属性参数：

你的命令存在两个核心配置错误，直接导致输出栅格尺寸不符、出现偏移：

• 手动传入的-te范围值精度不足。你填写的边界坐标仅保留了4位小数，和参考栅格真实仿射变换参数存在微小误差，该误差在-tap对齐规则下，会直接导致输出栅格多生成1行、1列像元。
• 对-tap参数的作用理解有误。-tap仅会强制输出栅格的范围边界对齐到像元分辨率的整数倍网格，不会自动匹配参考栅格的实际像元原点，你手动输入的边界本身和参考栅格原点存在偏差，加-tap只会将偏差放大到整像元级别。
• 冗余参数引入隐式误差：输入、输出栅格坐标系完全一致（均为EPSG:32630），重复指定-s_srs和-t_srs属于无效操作，可能触发不必要的坐标重计算，引入额外偏移。

不要手动复制粘贴参考栅格的范围、分辨率参数，直接通过GDAL接口读取参考栅格的精确仿射变换参数传入，从根源上避免手动取值的精度误差：

命令行版本 #

1. 首先读取参考栅格的完整精确参数：

   gdalinfo <参考栅格文件路径>
   

   从输出结果中提取4类核心值：
   • 左上角原点坐标：Origin = (x_min, y_max)
   • 像元分辨率：Pixel Size = (x_res, y_res)（注意y方向分辨率通常为负值）
   • 栅格尺寸：Size is x_size, y_size
   • 计算精确边界：x_max = x_min + x_size * x_res，y_min = y_max + y_size * y_res
2. 调用gdalwarp时传入精确到双精度浮点数满精度的边界值，保留-tap参数，根据数据类型选择重采样方法（分类数据用near，连续数据用bilinear/cubic）：

   gdalwarp -tr <x_res> <y_res的绝对值> -te <x_min> <y_min> <x_max> <y_max> -tap -r near -of GTiff <待对齐栅格路径> <输出栅格路径>
   

Python脚本版本 #

封装自动化脚本时不需要调用命令行，直接通过GDAL Python API读取参考参数传入Warp接口即可，示例代码如下：

from osgeo import gdal

# 读取参考栅格的精确空间参数
ref_dataset = gdal.Open(ref_raster_path, gdal.GA_ReadOnly)
ref_geotransform = ref_dataset.GetGeoTransform()
x_res = ref_geotransform[1]
y_res = abs(ref_geotransform[5])
x_min = ref_geotransform[0]
y_max = ref_geotransform[3]
x_size = ref_dataset.RasterXSize
y_size = ref_dataset.RasterYSize
x_max = x_min + x_size * x_res
y_min = y_max - y_size * y_res
# 关闭参考数据集释放内存
ref_dataset = None

# 配置Warp参数执行对齐
warp_options = gdal.WarpOptions(
    xRes=x_res,
    yRes=y_res,
    outputBounds=(x_min, y_min, x_max, y_max),
    targetAlignedPixels=True,  # 对应命令行-tap参数
    resampleAlg=gdal.GRA_NearestNeighbour,
    format='GTiff'
)
gdal.Warp(output_raster_path, input_raster_path, options=warp_options)

通过以上方法输出的栅格，行列数、像元分辨率、空间范围与参考栅格完全一致，加载后可完美重叠，直接转换为numpy数组时不会出现索引错位问题。

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