使用dask-image圆形足迹执行最大滤波时的内存错误问题

针对30000×50000的7GB GeoTIFF，使用dask_image.ndfilters.maximum_filter()设置半径>100的圆形footprint时调用.compute()触发MemoryError，核心原因是：

• 大半径邻域会让每个计算分块需要加载大量边界扩展数据，内存占用陡增；
• 直接.compute()会将整个结果数组加载到内存，远超单节点内存承载能力。

以下是具体优化方案：

1. 优化输入数组分块策略 #

默认分块可能过大或过小，手动设置匹配邻域半径的分块大小，同时确保分块的边界padding等于邻域半径，避免跨块计算时重复加载过多数据：

• 分块大小建议设为1000×1000~2000×2000（根据内存调整，块太小会增加调度开销）；
• 替代dask_image.imread的默认分块，用rasterio手动读取分块，保证分块逻辑与计算需求匹配。

示例代码：

import rasterio
from rasterio.windows import Window
import dask.array as da

def read_tiff_as_dask(raster_path, block_size=1000):
    with rasterio.open(raster_path) as src:
        height, width = src.height, src.width
        # 生成所有计算分块的窗口
        windows = [
            Window(col_off, row_off, min(block_size, width-col_off), min(block_size, height-row_off))
            for row_off in range(0, height, block_size)
            for col_off in range(0, width, block_size)
        ]
        # 逐个读取分块为延迟加载的dask数组
        chunks = [da.from_delayed(
            lambda w=w: src.read(window=w),
            shape=(src.count, w.height, w.width),
            dtype=src.dtype
        ) for w in windows]
        # 拼接为完整的dask数组
        grid_cols = (width + block_size -1) // block_size
        dask_arr = da.block([[chunks[i*grid_cols +j] for j in range(grid_cols)] 
                             for i in range((height + block_size -1)//block_size)])
        return da.squeeze(dask_arr, axis=0)

# 加载数据，设置分块大小为1000×1000
raster_data = read_tiff_as_dask("data_file.tif", block_size=1000)

2. 避免直接.compute()，分块导出GeoTIFF #

直接导出分块结果到GeoTIFF，无需将整个数组加载到内存。使用rasterio结合dask的map_blocks实现边算边写：

示例代码：

import rasterio

def write_dask_geotiff(dask_arr, output_path, src_path):
    with rasterio.open(src_path) as src:
        profile = src.profile.copy()
        # 同步输出分块与输入分块大小，开启压缩减少文件体积
        profile.update(
            blockxsize=dask_arr.chunks[1][0],
            blockysize=dask_arr.chunks[0][0],
            tiled=True,
            compress='deflate'
        )
        # 分块写入结果
        with rasterio.open(output_path, 'w', **profile) as dst:
            def write_block(block, window):
                dst.write(block, window=window)
                return None
            # 遍历所有分块窗口并写入
            for i, (row_slice, col_slice) in enumerate(zip(*dask_arr.chunks)):
                row_off = row_slice.start
                col_off = col_slice.start
                window = Window(col_off, row_off, col_slice.stop-col_slice.start, row_slice.stop-row_slice.start)
                block = dask_arr[row_slice, col_slice].compute()
                write_block(block, window)

# 计算最大值后直接导出，跳过全量compute()
radius = 300
# 直接用numpy生成圆形footprint，无需转为dask数组
circle_filter = np.zeros((2*radius+1, 2*radius+1), dtype=int)
y, x = np.ogrid[-radius:radius+1, -radius:radius+1]
mask = x**2 + y**2 <= radius**2
circle_filter[mask] = 1

max_value = dask_image.ndfilters.maximum_filter(raster_data, footprint=circle_filter, mode='nearest')
write_dask_geotiff(max_value, "output_max.tif", "data_file.tif")

3. 其他细节优化 #

• footprint无需转为Dask数组：直接使用NumPy数组的圆形掩码即可，dask_image.ndfilters会自动适配，避免额外调度开销；
• 调整集群内存参数：降低memory_target_fraction（比如设为0.8），给系统留足内存余量，同时根据CPU核心数调整worker数量：

  cluster = LocalCluster(
      n_workers=2,
      threads_per_worker=2,
      memory_target_fraction=0.8,
      memory_limit='16GB'
  )
  client = Client(cluster)
  

• 保留原数据类型：如果原GeoTIFF是整数类型，输出时保持一致，避免浮点型带来的内存翻倍。

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