优化栅格DEM洪水填充算法：快速识别8连通淹没区域

我拥有包含高程信息的栅格数字高程模型（DEM），以及若干带指定高程的点（格式为Point(x, y, h)，其中h为高程值）。需求为：从给定点位置出发，找出所有8连通且高程低于h的单元格。

示例 #

给定矩阵M和Point(2,2,2.1)，期望得到结果矩阵MS：

import numpy as np

M = np.array([
    [2, 6, 6, 6, 6],
    [6, 2, 6, 6, 6],
    [6, 6, 2, 6, 6],
    [6, 6, 6, 6, 6],
    [6, 6, 6, 6, 6]
])

MS = np.array([
    [True, False, False, False, False],
    [False, True, False, False, False],
    [False, False, True, False, False],
    [False, False, False, False, False],
    [False, False, False, False, False]
], dtype=bool)

已实现的迭代式洪水填充算法 #

import numpy as np

def flood_fill(image, x, y, threshold):
    """
    Calculates inundated cells. Cells get inundated if the threshold value is higher
    than the starting point elevation.

    Parameters
    ----------
    image: numpy.ndarray
        array of elevations
    x : int
        starting coordinate x in pixel coordinates
    y : int
        starting coordinate y in pixel coordinates
    threshold : float
        threshold elevation.

    Returns
    -------
    filled : numpy.ndarray
        array with boolean values where True means the cells is flooded
    """
    
    filled = np.zeros_like(image, dtype=bool)
    toFill = []
    toFill.append((x,y))
    
    while len(toFill) > 0:
        (x,y) = toFill.pop()
        
        if x < 0 or y < 0 or x >= image.shape[0] or y >= image.shape[1]:
            continue
        
        if filled[x, y] or image[x, y] > threshold:
            continue
        
        filled[x, y] = True
        toFill.append((x - 1, y))
        toFill.append((x + 1, y))
        toFill.append((x, y - 1))
        toFill.append((x, y + 1))
        toFill.append((x-1, y - 1))
        toFill.append((x-1, y + 1))
        toFill.append((x+1, y - 1))
        toFill.append((x+1, y + 1))
        
    return filled

性能问题与需求 #

该算法可正常运行，但处理效率较低：处理5000×3000的DEM单次耗时15秒，而我需处理多个点；递归版本则因栈溢出无法处理大尺寸图像。现寻求算法加速方案，或具备相同功能的Python库。

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