checking the connectivity of cells in a raster

12-14-2017 01:04 PM
New Contributor II

Hello, can anybody help with a python script that checks if cell A is connected to its neighboring cells and if it flooded then the cell connecting cells should be flooded too. I tried the code below but i cant seem to figure out how to in-cooperate my data into the script. My data include a DEM from which i will check inundated zone:


import heapq
import numpy as np
from scipy import ndimage

def flood_fill(test_array, four_way=False):
input_array = np.copy(test_array)
input_rows, input_cols = input_array.shape
# Set h_max to a value larger than the array maximum to ensure
#   that the while loop will terminate
h_max = np.max(input_array*2.0)
# Build mask of cells with data not on the edge of the image
# Use 3x3 square structuring element
inside_mask = ndimage.binary_erosion(
 structure=ndimage.generate_binary_structure(2, 2).astype(
# Initialize output array as max value test_array except edges
output_array = np.copy(input_array)
output_array[inside_mask] = h_max
# Build priority queue and place edge pixels into priority queue
# Last value is flag to indicate if cell is an edge cell
put = heapq.heappush
get = heapq.heappop
fill_heap = [
    (output_array[t_row, t_col], int(t_row), int(t_col), True)
    for t_row, t_col in np.transpose(np.where(edge_mask))]

def neighbors(row, col, four_way=False):
    """Return indices of adjacent cells"""
    if four_way:
        return [
            (row - 1, col), (row, col + 1), 
            (row + 1, col), (row, col - 1)]
        return [
            (row - 1, col), (row - 1, col + 1), 
            (row, col + 1), (row + 1, col + 1), 
            (row + 1, col), (row + 1, col - 1), 
            (row, col - 1), (row - 1, col - 1)]

# Iterate until priority queue is empty
while True:
    h_crt, t_row, t_col, edge_flag = get(fill_heap)
except IndexError:
for n_row, n_col in neighbors(t_row, t_col, four_way):
    # Skip cell if outside array edges
    if edge_flag:
            if not inside_mask[n_row, n_col]:
except IndexError:

Tags (2)
0 Kudos
6 Replies
MVP Esteemed Contributor

what constitutes flooded?

0 Kudos
Esri Esteemed Contributor

Maybe the code in this old post helps: 

MVP Esteemed Contributor

Grief! forgot about Bill's avenue incarnation and its translation!

New Contributor II

But i dont have the hard coded values every data is from the user

0 Kudos
Esri Esteemed Contributor

The value for nFloodMax in the script could be made variable as a parameter for the script or the tool that you create. See Understanding script tool parameters—Help | ArcGIS Desktop  and Accessing parameters in a script tool—Help | ArcGIS Desktop 

0 Kudos
New Contributor II

I have sea level rise scenarios, so if a cell is below the sea level it is  flooded so i want to check if the other neighboring are connected to the  flooded cell.

0 Kudos