如何利用GDAL重采样实现栅格图层下采样并计算95分位数？

我们有1米分辨率的数字表面模型（DSM），想要生成50米分辨率的版本，要求每个低分辨率栅格单元的值是对应高分辨率区域的95分位数。已知gdalwarp支持中位数（med）、四分位数（q1/q3）重采样，但不支持自定义分位数设置，以下是几种可行的实现方法：

方法一：GDAL Python API 自定义计算 #

通过GDAL的Python接口逐窗口读取高分辨率栅格数据，计算95分位数后写入低分辨率输出栅格：

from osgeo import gdal
import numpy as np

# 配置路径
input_tif = "dsm_1m.tif"
output_tif = "dsm_50m_95p.tif"

# 打开输入栅格
src_ds = gdal.Open(input_tif)
src_band = src_ds.GetRasterBand(1)
src_gt = src_ds.GetGeoTransform()
src_proj = src_ds.GetProjection()
nodata_val = src_band.GetNoDataValue()

# 定义输出栅格参数
output_res = 50  # 50米分辨率
output_cols = int((src_ds.RasterXSize * src_gt[1]) / output_res)
output_rows = int((src_ds.RasterYSize * abs(src_gt[5])) / output_res)

# 创建输出栅格
driver = gdal.GetDriverByName("GTiff")
dst_ds = driver.Create(
    output_tif, output_cols, output_rows, 1, gdal.GDT_Float32,
    options=["COMPRESS=LZW"]
)
dst_ds.SetGeoTransform((src_gt[0], output_res, src_gt[2], src_gt[3], src_gt[4], -output_res))
dst_ds.SetProjection(src_proj)
dst_band = dst_ds.GetRasterBand(1)
dst_band.SetNoDataValue(nodata_val)

# 逐窗口计算95分位数
win_xsize = int(output_res / src_gt[1])
win_ysize = int(output_res / abs(src_gt[5]))

for y in range(output_rows):
    for x in range(output_cols):
        # 计算输入栅格中的窗口位置
        src_x = x * win_xsize
        src_y = y * win_ysize
        # 处理边缘窗口，避免超出栅格范围
        actual_xsize = min(win_xsize, src_ds.RasterXSize - src_x)
        actual_ysize = min(win_ysize, src_ds.RasterYSize - src_y)
        
        # 读取窗口数据并过滤无效值
        data = src_band.ReadAsArray(src_x, src_y, actual_xsize, actual_ysize)
        valid_data = data[data != nodata_val]
        
        if len(valid_data) == 0:
            result = nodata_val
        else:
            result = np.percentile(valid_data, 95)
        
        # 写入当前输出栅格单元
        dst_band.WriteArray(np.array([[result]]), x, y)

# 释放资源
dst_band.FlushCache()
src_ds = None
dst_ds = None

方法二：R语言 raster 包快速实现 #

利用R的raster包的aggregate函数，直接指定分位数计算逻辑：

library(raster)

# 读取高分辨率DSM
dsm_high <- raster("dsm_1m.tif")
# 按50米分辨率聚合，计算95分位数（忽略NA值）
dsm_low_95p <- aggregate(
    dsm_high, 
    fact = 50,  # 聚合因子（50=50米/1米）
    fun = function(x) quantile(x, 0.95, na.rm = TRUE)
)
# 保存结果
writeRaster(dsm_low_95p, "dsm_50m_95p.tif", format = "GTiff", overwrite = TRUE)

方法三：GRASS GIS 命令行工具 #

使用GRASS的r.resamp.stats工具直接指定分位数参数：

1. 导入高分辨率栅格：

r.in.gdal input=dsm_1m.tif output=dsm_1m

2. 设置工作区域为50米分辨率：

g.region raster=dsm_1m res=50

3. 计算95分位数下采样：

r.resamp.stats input=dsm_1m output=dsm_50m_95p method=quantile quantile=0.95

4. 导出结果为TIFF：

r.out.gdal input=dsm_50m_95p output=dsm_50m_95p.tif format=GTiff

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