How to identify or determine Contiguous Grid Cells of a extreme events using daily NetCDF data

174 Views Asked by At

I've been working with daily precipitation netcdf data with dimensions time:1322, lat:68, lon:136 over my domain. For each time step (date), I want to identify/determine contiguous grid cells, i.e., grid cells in a date in which daily precipitation amount excced my threshod (95th percentile value) that is further connected to another grid cell that is also greater than threshols using Flood Fill Algorithm function

I have written and deploye Flood Fill Algorithm function to visit each cells if the cell value if higher than threshold but I don't know how to determine if the cells is connected to another cells either in up, down, left or right directions that is/are greater than threshold and assign value of 1. Likewise, if it is less than threshold replace the value of the cell with NAN or 0

I need help/suggestion to complete this code to perform the above task and append the result as a new variable to the loaded data for further analysis or another suitable approach to perform this task

here is my code:

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

PS: the approach of hilorywilsmart and solution provided by @Michael Delgado does not fully satify condition neighbourhood in contiguous grid cells analysis in extreme weather events

0

There are 0 best solutions below