如何在Rasterio中实现raster1与raster2的差值计算（匹配raster1范围）

问题场景 #

需要计算raster1.tif与raster2.tif的差值，其中raster1.tif范围更大且完全覆盖raster2.tif，要求输出栅格的范围与raster1.tif一致。

原代码如下：

import rasterio 

with rasterio.open('raster1.tif') as src1, rasterio.open('raster2.tif') as src2:
    data1 = src1.read(1)
    data2 = src2.read(1)
    data = data1 - data2

    with rasterio.open("output.tif", 'w', **src1.meta) as dst:
        dst.write(data,1)

运行后报错：

ValueError: operands could not be broadcast together with shapes (790,1554) (57,67) 

错误原因是两个栅格的行列数（形状）不匹配，无法直接进行数值运算。

补充说明：已在ArcGIS Pro中通过栅格计算器实现需求，逻辑为：

Con(IsNull('raster2.tif'), 'raster1.tif', 'raster1.tif' - 'raster2.tif')

解决方案 #

核心思路是将raster2.tif的数据对齐到raster1.tif的栅格网格（包括坐标系、分辨率、尺寸），再实现与ArcGIS Con工具等效的逻辑：当raster2无数据时保留raster1的值，否则计算两者差值。

实现代码如下：

import rasterio
from rasterio.enums import Resampling
import numpy as np

with rasterio.open('raster1.tif') as src1, rasterio.open('raster2.tif') as src2:
    # 获取raster1的元数据、数据和无数据值
    meta1 = src1.meta
    data1 = src1.read(1)
    nodata1 = meta1['nodata']
    nodata2 = src2.nodata

    # 计算raster2边界在raster1栅格中的窗口范围
    window = src1.window(*src2.bounds)
    # 读取raster2对应窗口的数据，并重采样到raster1的分辨率
    data2 = src2.read(
        1,
        window=window,
        out_shape=(src1.window_height(window), src1.window_width(window)),
        resampling=Resampling.bilinear  # 可根据需求替换为Resampling.nearest等
    )

    # 创建与raster1同形状的数组，初始填充raster1的无数据值
    data2_full = np.full_like(data1, nodata1)
    # 获取窗口对应的行列切片范围
    row_slice = slice(window.row_off, window.row_off + window.height)
    col_slice = slice(window.col_off, window.col_off + window.width)
    # 将raster2的数据填充到对应位置
    data2_full[row_slice, col_slice] = data2

    # 实现ArcGIS Con逻辑：raster2无数据时取raster1，否则计算差值
    mask = (data2_full == nodata1) | (data2_full == nodata2) | np.isnan(data2_full)
    result = np.where(mask, data1, data1 - data2_full)

    # 写入输出栅格
    with rasterio.open("output.tif", 'w', **meta1) as dst:
        dst.write(result, 1)

代码关键点说明 #

• 窗口计算：通过src1.window(*src2.bounds)定位raster2在raster1栅格中的精确位置，确保数据对齐。
• 重采样：使用out_shape和resampling参数将raster2的分辨率匹配到raster1，避免尺寸不兼容问题。
• 无数据处理：通过np.where实现条件判断，覆盖raster2未覆盖区域及自身无数据的情况，与ArcGIS Con逻辑完全等效。

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