将人口栅格重采样投影匹配降水栅格时结果异常的排查

核心问题根源 #

你遇到的区域统计值异常（总和为0或严重偏离），主要是两个关键错误导致：

• 人口栅格是面权重数据（每个栅格对应区域的总人口），但你使用了针对连续数值（如降水）的双线性插值重采样，破坏了人口总量守恒。
• 部分代码未完全匹配目标降水栅格的空间参数（CRS、分辨率、范围、对齐方式），导致栅格错位。

关键修正要点 #

1. 必须使用求和重采样 #

人口数据重采样时，必须选择**求和（Resampling.sum）**方法：将目标栅格范围内的所有原始1km栅格人口相加，才能保证区域总量准确。双线性、最近邻等方法仅适用于温度、降水这类连续变量。

2. 完整匹配目标栅格空间参数 #

重投影+重采样需同时匹配降水栅格的CRS、分辨率、边界、栅格对齐方式，避免错位导致统计错误。

修正后的代码示例 #

方法一：rioxarray.reproject_match（推荐，自动匹配所有参数） #

import rioxarray
import xarray

def print_raster(raster):
    print(
        f"shape: {raster.rio.shape}\n"
        f"resolution: {raster.rio.resolution()}\n"
        f"bounds: {raster.rio.bounds()}\n"
        f"sum: {raster.sum().item()}\n"
        f"CRS: {raster.rio.crs}\n"
    )

# 加载原始人口栅格与目标降水栅格
xds_pop = rioxarray.open_rasterio('pop_m4_2010.tif')
xds_prcp = rioxarray.open_rasterio('prp_raster.tiff')

# 核心：使用sum重采样保证人口总量守恒
xds_pop_reproj = xds_pop.rio.reproject_match(
    xds_prcp,
    resampling=rioxarray.enums.Resampling.sum
)

# 验证总量一致性
print("原始人口总和：", xds_pop.sum().item())
print("重采样后人口总和：", xds_pop_reproj.sum().item())

# 保存结果
xds_pop_reproj.rio.to_raster("reproj_pop_sum.tif")

方法二：rasterio手动实现 #

修复原代码的参数错误，改用求和重采样：

import rasterio
from rasterio.warp import reproject, Resampling

# 获取目标降水栅格的所有空间参数
with rasterio.open('prp_raster.tiff') as prcp_ds:
    dst_crs = prcp_ds.crs
    dst_transform = prcp_ds.transform
    dst_width = prcp_ds.width
    dst_height = prcp_ds.height

# 处理人口栅格重投影与重采样
with rasterio.open('pop_m4_2010.tif') as src_ds:
    # 更新元数据为目标参数
    kwargs = src_ds.meta.copy()
    kwargs.update({
        'crs': dst_crs,
        'transform': dst_transform,
        'width': dst_width,
        'height': dst_height,
        'dtype': rasterio.float32  # 避免大数值整数溢出
    })

    with rasterio.open('pop_resampled_sum.tif', 'w', **kwargs) as dst_ds:
        for i in range(1, src_ds.count + 1):
            reproject(
                source=rasterio.band(src_ds, i),
                destination=rasterio.band(dst_ds, i),
                src_crs=src_ds.crs,
                dst_crs=dst_crs,
                src_transform=src_ds.transform,
                dst_transform=dst_transform,
                resampling=Resampling.sum
            )

修正NetCDF转TIFF代码（避免多波段问题） #

原代码导出的是全时间序列多波段TIFF，建议选择单月数据导出，便于后续匹配：

import xarray as xr 
import rioxarray as rio

prcp_file = xr.open_dataset('nclimgrid_prcp.nc')
# 选取第一个时间步的降水数据
prp = prcp_file['prcp'].isel(time=0)
prp = prp.rio.set_spatial_dims(x_dim='lon', y_dim='lat')
prp.rio.write_crs("epsg:4326", inplace=True)
prp.rio.to_raster('prp_single_month.tiff')

验证步骤 #

1. 检查重采样前后的人口总和是否一致，不一致则说明重采样方法错误。
2. 在QGIS中叠加县矢量与重采样后的人口栅格，用Zonal Statistics工具计算区域总和，对比已知县人口数据验证准确性。
3. 查看栅格的bounds和transform参数，确认所有栅格空间范围完全对齐。

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