Automatically determine downstream grid cell ID

419
1
07-22-2023 12:29 PM
TaylorJordan
New Contributor

Hello,

I'll start by explaining a little about the project I'm working on.  I'm trying to analyze a watershed (~60 sq. km) using a lumped rainfall-runoff model.  To improve the resolution of the watershed, I've created a grid of polygons with a cell size of 100mX100m (each with an elevation assigned to it).  The grid cells act as subcatchments in the model and each subcatchment needs to have an outlet assigned to it.  With the size of my watershed and number of grid cells, I need the outlet IDs to be automatically determined as opposed to doing it all myself. 

I'm sure there's a way for ArcGIS Pro to automatically tell me the name/ID of the neighboring cell that receives runoff.  The Flow Direction tool doesn't quite tell me what I need, since it just determines which direction each cell drains.  Does anyone know of a way to do this?  Ideally, the final result would be a column in the Attribute Table for my polygon grid called "Outlet", which would have the name of the neighboring cell that that cell drains to.

Hopefully that all makes sense.  I can't seem to find a solution to this in the rainfall-runoff model I'm using, so I'm hoping to process it in ArcGIS first and then import into my model.  That way, the outlet will already be defined for each cell to use in the model.  For reference, I'm using ArcGIS Pro 3.0.2.

I appreciate any help.

Thanks!

Taylor

1 Reply
JohannesLindner
MVP Frequent Contributor

Hmmm, this works, but I really hope I missed an easier solution.

 

# extract your watershed from the larger digital elevation model
ws_dem = arcpy.sa.ExtractByMask("DEM", "Cutter")

# resample to 100x100
ws_dem_resampled = arcpy.management.Resample(ws_dem, "ws_dem_resampled", "100 100")

# fill and flow dir
ws_dem_filled = arcpy.sa.Fill(ws_dem_resampled)
flow_dir = arcpy.sa.FlowDirection(ws_dem_filled)

# create an id raster (numbered sequentially)
cell_points = arcpy.conversion.RasterToPoint(ws_dem_filled, "cell_points", "Value")
id_raster = arcpy.Raster(arcpy.conversion.PointToRaster(cell_points, "pointid", "id_raster", cellsize=flow_dir))

# calculate downstream id
directions = {
    1: (1, 0),
    2: (1, 1),
    4: (0, 1),
    8: (-1, 1),
    16: (-1, 0),
    32: (-1, -1),
    64: (0, -1),
    128: (1, -1)
}
ds_id_raster = arcpy.Raster(id_raster.getRasterInfo())
for x, y in ds_id_raster:
    d = flow_dir[x, y]
    i = id_raster[x,y ]
    try:
        dy, dx = directions[d]
        i_ = id_raster[x+dx, y+dy]
    except:
        continue
    ds_id_raster[x, y] = i_
ds_id_raster.save("ds_id_raster")

# convert raster to polygons
cells = arcpy.conversion.RasterToPolygon(id_raster, "Cells", "NO_SIMPLIFY")
# extract values into the points and join
arcpy.sa.ExtractMultiValuesToPoints(cell_points, [(id_raster, "Cell"), (ds_id_raster, "CellDS"), (flow_dir, "FlowDir"), (ws_dem_filled, "Elevation")])
arcpy.management.JoinField(cells, "gridcode", cell_points, "pointid", ["Cell", "CellDS", "FlowDir", "Elevation"])

 

This script will:

  • resample your dem

JohannesLindner_0-1690219954987.pngJohannesLindner_1-1690219966161.png

 

  • calculate flow direction

JohannesLindner_2-1690220026698.png

 

  •  Give a unique id to each raster cell

JohannesLindner_3-1690220137889.png

  • Use the flow direction to get the downstream cell id

JohannesLindner_4-1690220232474.png

 

  • Extract all this information into a polygon grid

JohannesLindner_5-1690220337986.png

 

 


Have a great day!
Johannes