如何使用Rasterio对多幅栅格影像求和并解决水陆边界异常块问题

• 现有100余幅空间范围一致的二值栅格影像，像元值仅为0或1，需对所有影像对应位置的像元值求和，输出单幅结果栅格。
• 采用逐次累加的代码实现后，水陆交界位置出现异常：预期水体区域累加值更高，但紧邻水体的陆地区域出现大块异常值区块，与预期效果偏差较大。
• 异常效果示意图：
  
• 原有实现代码：

def Summing(images_to_sum):
    """
    Creates one image out of all images in given list.
    """
    path_out = r"C:\Documents\UFV.jp2"
    images_used = 2

    def Summ(path_1, path_2):
        """
        Sums up the two images from given paths.
        """
        with rio.open(path_1) as src_1:
            array_1 = src_1.read()
            profile = src_1.profile

            with rio.open(path_2) as src_2:
                array_2 = src_2.read()

                # Sum the two
                result = array_1 + array_2

        with rio.open(path_out, 'w', **profile) as dst:
            dst.write(result)

    # Creating first image, so others can be added
    Summ(images_to_sum[0], images_to_sum[1])
    images_to_sum.pop(0)
    images_to_sum.pop(0)

    # Adding all the other images for this granule
    for image in images_to_sum:
        Summ(path_out, image)
        images_used += 1

核心问题出在逐次写入中间结果的逻辑，两个关键bug直接导致异常：

1. 同路径读写冲突
   循环累加阶段，每次传入Summ函数的第一个参数就是最终输出路径path_out。此时读取path_out的上下文管理器（with rio.open(path_1) as src_1）还未退出，文件句柄没有释放，就直接以写入模式打开同一路径的文件覆盖写入，在Windows系统下会触发文件锁冲突，导致读取到的数组内容损坏，多次循环后损坏值累积，在像素值跳变剧烈的水陆交界区域会形成明显异常块。
2. JP2有损压缩误差累积
   JP2格式默认采用有损压缩算法，每次写入中间结果都会引入像素值偏差。代码中每完成一次两图相加就写入一次JP2，经过近百次循环，压缩误差反复叠加，在高对比度的边缘位置会产生肉眼可见的异常区块。
3. 潜在风险：代码直接复用原始栅格的配置文件profile，没有修改数据类型。如果原始二值栅格是1位深度/布尔类型，直接做加法会出现值溢出翻转，不过这类问题通常会在全图范围出现，不是边缘异常的核心诱因。

完全摒弃逐次写入中间文件的逻辑，采用一次读取、内存累加、单次写入的流程，从根源上避免文件读写冲突和压缩误差问题。如果单幅栅格文件过大无法全部载入内存，可采用分块读取累加的方式，中间结果始终保存在内存中，不要落地为压缩格式文件。
优化后代码如下：

import rasterio as rio
import numpy as np

def Summing(images_to_sum, path_out=r"C:\Documents\UFV.tif"):
    # 先读取第一幅影像，获取基准元信息
    with rio.open(images_to_sum[0]) as src_base:
        profile = src_base.profile.copy()
        # 校验所有影像的行列数、仿射参数、投影是否一致
        for img_path in images_to_sum[1:]:
            with rio.open(img_path) as src_check:
                if (src_check.width != src_base.width or
                    src_check.height != src_base.height or
                    src_check.transform != src_base.transform or
                    src_check.crs != src_base.crs):
                    raise ValueError(f"影像{img_path}与基准影像空间参考/范围/分辨率不匹配")
        # 修改输出数据类型为uint16，足够存储100+景二值栅格的求和结果（最大值100+）
        profile.update(dtype='uint16', nodata=None)
        # 初始化累加数组
        sum_array = np.zeros((src_base.count, src_base.height, src_base.width), dtype=np.uint16)
        # 遍历所有影像累加
        for img_path in images_to_sum:
            with rio.open(img_path) as src:
                sum_array += src.read().astype(np.uint16)
    # 所有影像读取累加完成后，一次性写入结果，推荐用无损的GTiff格式存储避免压缩损失
    with rio.open(path_out, 'w', **profile) as dst:
        dst.write(sum_array)

如果必须输出JP2格式，写入时需要在profile中添加无损压缩参数：

profile.update(
    driver='JP2OpenJPEG',
    dtype='uint16',
    nodata=None,
    QUALITY=100,
    REVERSIBLE='YES' # 开启可逆无损压缩
)

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