如何利用逐日NetCDF数据识别极端事件的连续网格单元

我正在处理一套逐日降水NetCDF数据，维度为time:1322、lat:68、lon:136。需求是针对每个时间步（单日），通过洪水填充算法识别符合以下条件的连续网格单元：

• 当日该网格降水量超过对应网格的95百分位数阈值
• 与同满足阈值条件的四邻域（上下左右）网格相连

需要将这类连续单元赋值为1，低于阈值的单元赋值为0，最终将结果作为新变量加入原数据用于后续分析。

现有代码已实现洪水填充的递归框架，但无法正确判断单元与邻域的连通性，且之前hilorywilsmart和@Michael Delgado的方案无法满足极端天气事件分析的邻域条件。

现有代码 #

import numpy as np
import xarray as xr

# Load the daily rainfall data in netcdf format
path_gpcc = '/path_to_my_file'
ds = xr.open_dataset(path_gpcc)

# Extract the rainfall variable
rainfall = ds.precip

# Define the Flood Fill Algorithm function
def flood_fill(t, x, y):
    """
    Recursive function to perform Flood Fill Algorithm
    """
    # Set current cell to visited
    visited[t, x, y] = True
    
    # Check if current cell exceeds the threshold and has not been visited
    if rainfall[t, x, y] >= threshold[x, y] and not visited[t, x, y]:
        # Check neighboring cells
        if x > 0:
            flood_fill(t, x-1, y)
        if x < rainfall.shape[1]-1:
            flood_fill(t, x+1, y)
        if y > 0:
            flood_fill(t, x, y-1)
        if y < rainfall.shape[2]-1:
            flood_fill(t, x, y+1)
            
# Calculate the 90th percentile value of the rainfall for each grid cell at all time steps
threshold = np.percentile(rainfall, 90, axis=0)

# Create an array to keep track of visited cells for all time steps
visited = np.zeros_like(rainfall, dtype=bool)

# List to store the contiguous grid cells
contiguous_cells = []

# Loop through each time step in the data period
for t in range(rainfall.shape[0]):
    # Loop through each cell in the grid for the current time step
    for x in range(rainfall.shape[1]):
        for y in range(rainfall.shape[2]):
            # Check if cell has been visited
            if not visited[t, x, y]:
                # Start Flood Fill Algorithm
                flood_fill(t, x, y)
                # Add contiguous grid cells to list
                contiguous_cells.append(np.where(visited))
                # Reset visited array
                visited = np.zeros_like(rainfall, dtype=bool)

# Create a new variable with the contiguous grid cells for all time steps
contiguous_cells_var = xr.Variable(dims=('time', 'lat', 'lon'), data=np.zeros_like(rainfall, dtype=int))
for cells in contiguous_cells:
    contiguous_cells_var[cells] = 1

修正后的代码及说明 #

关键修改点 #

1. 调整洪水填充函数逻辑顺序：先判断单元是否符合条件（未访问+超过阈值），再标记已访问，避免提前标记导致逻辑错误
2. 按时间步独立处理：每个时间步单独初始化访问标记和结果数组，避免跨时间步污染
3. 修正阈值计算为需求的95百分位数
4. 直接生成结果数组，简化存储逻辑
5. 将结果转为匹配原数据维度的xarray变量，加入原数据集

import numpy as np
import xarray as xr

# 加载数据
path_gpcc = '/path_to_my_file'
ds = xr.open_dataset(path_gpcc)
rainfall = ds.precip

# 计算每个网格的95百分位数阈值（按时间轴统计）
threshold = np.percentile(rainfall, 95, axis=0)

# 初始化结果数组，默认值为0，后续将连通单元设为1
result = np.zeros_like(rainfall, dtype=int)

def flood_fill(t, x, y, visited_t, rainfall_t, threshold_t, result_t):
    """
    针对单个时间步的洪水填充函数
    参数:
        t: 当前时间步索引
        x: 纬度索引
        y: 经度索引
        visited_t: 当前时间步的访问标记数组（2D）
        rainfall_t: 当前时间步的降水数据（2D）
        threshold_t: 网格阈值数组（2D）
        result_t: 当前时间步的结果数组（2D）
    """
    # 边界判断：超出网格范围则返回
    if x < 0 or x >= rainfall_t.shape[0] or y < 0 or y >= rainfall_t.shape[1]:
        return
    # 如果已访问或未超过阈值，返回
    if visited_t[x, y] or rainfall_t[x, y] < threshold_t[x, y]:
        return
    
    # 标记为已访问，并将结果设为1
    visited_t[x, y] = True
    result_t[x, y] = 1
    
    # 递归遍历四邻域
    flood_fill(t, x-1, y, visited_t, rainfall_t, threshold_t, result_t)  # 上
    flood_fill(t, x+1, y, visited_t, rainfall_t, threshold_t, result_t)  # 下
    flood_fill(t, x, y-1, visited_t, rainfall_t, threshold_t, result_t)  # 左
    flood_fill(t, x, y+1, visited_t, rainfall_t, threshold_t, result_t)  # 右

# 遍历每个时间步
for t in range(rainfall.shape[0]):
    # 获取当前时间步的降水数据和阈值
    rainfall_t = rainfall[t].values
    threshold_t = threshold.values
    # 初始化当前时间步的访问标记和结果数组
    visited_t = np.zeros_like(rainfall_t, dtype=bool)
    result_t = np.zeros_like(rainfall_t, dtype=int)
    
    # 遍历当前时间步的所有网格
    for x in range(rainfall_t.shape[0]):
        for y in range(rainfall_t.shape[1]):
            # 如果未访问且超过阈值，触发洪水填充
            if not visited_t[x, y] and rainfall_t[x, y] >= threshold_t[x, y]:
                flood_fill(t, x, y, visited_t, rainfall_t, threshold_t, result_t)
    
    # 将当前时间步的结果存入总结果数组
    result[t] = result_t

# 将结果转为xarray变量，匹配原数据的维度和坐标
contiguous_cells = xr.DataArray(
    result,
    dims=('time', 'lat', 'lon'),
    coords={'time': ds.time, 'lat': ds.lat, 'lon': ds.lon},
    name='contiguous_extreme_rain'
)

# 将新变量加入原数据集
ds['contiguous_extreme_rain'] = contiguous_cells

# 可选：保存结果到新的NetCDF文件
# ds.to_netcdf('/path_to_save_result.nc')

递归深度问题解决方案 #

如果连通区域过大，递归可能触发RecursionError，可改用迭代版洪水填充（栈实现）：

def flood_fill_iterative(t, x, y, visited_t, rainfall_t, threshold_t, result_t):
    stack = [(x, y)]
    while stack:
        cx, cy = stack.pop()
        if cx < 0 or cx >= rainfall_t.shape[0] or cy < 0 or cy >= rainfall_t.shape[1]:
            continue
        if visited_t[cx, cy] or rainfall_t[cx, cy] < threshold_t[cx, cy]:
            continue
        visited_t[cx, cy] = True
        result_t[cx, cy] = 1
        # 将四邻域加入栈
        stack.append((cx-1, cy))
        stack.append((cx+1, cy))
        stack.append((cx, cy-1))
        stack.append((cx, cy+1))

替换原递归函数即可，适合处理大范围连通区域。

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