Watershed Delineation in ArcGIS Pro: the 10-step Hydrology workflow

Every watershed study, flood model, RUSLE erosion map and morphometric paper starts the same way: a DEM goes in, a stream network and a set of basin polygons come out. The sequence was formalised by Jenson and Domingue in 1988, and ArcGIS Pro's Spatial Analyst Hydrology toolset still runs it in that order.

The tools are easy. What goes wrong is the settings: a pour point one cell off the channel, a stream threshold copied from a paper written for a different climate, a DEM clipped with no buffer. This guide gives you the full ten-step chain with the parameters I use on Indian basins, the flow-direction codes, threshold guidance by terrain, the morphometric formulas the outputs feed, and an arcpy script that runs steps one to seven unattended.

Before you start: the DEM

Match resolution to the question, not to the best file you can find. A 30 m DEM resolves first-order channels around 60 to 90 m wide; a 1 m LiDAR DEM resolves rills and terraces. For basins over 1000 km², Copernicus GLO-30 or SRTM 30 m are appropriate. For sub-basins under 100 km² and urban drainage, use a 5 to 10 m product.

For Indian river basins I default to Copernicus GLO-30 or CartoDEM, with LiDAR or UAV DEMs where they exist for micro-watershed work.

Prepare it

• Mosaic the tiles (Mosaic To New Raster, 32-bit float, 1 band).
• Project to the local UTM zone with bilinear resampling: UTM 43N (EPSG:32643) for 72°–78° E, 44N (EPSG:32644) for 78°–84° E, 45N (EPSG:32645) for 84°–90° E.
• Clip with Extract by Mask to an extent about 10 percent larger than the basin on every side. Edge cells with no buffer produce truncated watersheds.
• Inspect the histogram for voids and spikes, and write down source, date, resolution, vertical datum (EGM96 or EGM2008 for global products) and EPSG code. That paragraph goes straight into your methods.

Set the environment once

Analysis → Environments: workspace and scratch workspace in a file geodatabase; output coordinate system set to your UTM zone; cell size equal to the DEM; snap raster set to the DEM so every output is co-registered; extent equal to the DEM. Confirm the Spatial Analyst extension is licensed under Project → Licensing, or the toolset will not appear.

The ten steps

1. Fill

A sink is a cell lower than all eight neighbours; D8 cannot assign it a flow direction. Genuine closed depressions are rare in humid Indian terrain, so fill everything: leave Z limit blank. Then check what you filled with the Raster Calculator:

"DEM_fill" - "DEM_raw"

Expect small positive values, usually under 1 to 2 m, at filled sinks and zero elsewhere. Large filled areas mean voids or systematic DEM error. Investigate before going on. Running the Sink tool first, as a diagnostic, is good practice.

2. Flow Direction

Input: the filled DEM. Tick Force all edge cells to flow outward. Flow direction type: D8 for channel and watershed work; MFD for divergent hillslope flow and erosion modelling; D-Infinity for topographic wetness index.

D8 sends all flow from a cell to the steepest of its eight neighbours, with slope computed as elevation difference over distance (1 cell width for cardinal neighbours, √2 for diagonals). The output codes are powers of two:

Validate by symbolising unique values over a hillshade: a south-facing slope should be dominated by 4. Anything else means residual sinks or DEM errors.

3. Flow Accumulation

Input: the flow-direction raster, same flow type as step 2, output as float. Each cell now holds the count of upstream cells draining into it. The histogram is extremely right-skewed, so symbolise with a Standard Deviation stretch (n = 2) or a logarithmic stretch to see the channel network as bright lines. To convert counts to contributing area in square metres, multiply by the cell area:

// 30 m DEM: 30 × 30 = 900 m² per cell
"FlowAcc" * 900

4. Stream definition

Cells at or above a critical support area become streams; everything else becomes NoData.

Con("FlowAcc" >= 1000, 1)

The threshold is the decision reviewers will ask about. For a 30 m DEM:

Run it at three thresholds (for example 500, 1000, 2000), overlay them on Survey of India 1:50,000 sheets, keep the one that matches the blue lines, and report the threshold next to your drainage density. Drainage density without its threshold is meaningless.

5. Stream Link

Inputs: the stream raster and the flow-direction raster. Every segment between junctions gets a unique integer. Symbolise by unique values: segments should end at every confluence and at the outlet.

6. Stream Order

Method STRAHLER for morphometric work (headwaters are order 1; two streams of order n make order n+1; unequal orders keep the higher). SHREVE sums magnitudes at every junction and suits network-topology research. Record the maximum order and the segment count per order; typical Indian sub-basins of 100 to 1000 km² reach Strahler order 4 to 6.

7. Stream to Feature

Inputs: the stream-order raster and flow direction, with Simplify polylines ticked. Add a Length_m field with Calculate Geometry and summarise it by order. Those totals feed the bifurcation ratio, stream-length ratio and drainage density below.

8. Snap Pour Point

Create a point feature class Outlets in the DEM's projection with a unique integer field PourID, and digitise each outlet: a gauge, a culvert, a reservoir inlet. Snap them to the flow-accumulation raster with a snap distance of 100 m (3 to 4 cells at 30 m). Overlay the result on FlowAcc: every snapped point must sit on a bright channel cell. If not, raise the distance to 250 to 500 m and check it has not jumped to a neighbouring tributary.

9. Watershed

Inputs: flow direction and the snapped pour points (field Value for a raster, PourID for features). Convert to polygons with Raster to Polygon, Simplify polygons off, and add an Area_km2 field. A watershed of fewer than five cells means the pour point never reached the channel: go back to step 8.

10. Basin

No pour points needed: every natural basin that drains to a sink or the DEM edge gets an ID. Useful for reconnaissance of an unfamiliar terrain and for regional drainage characterisation. It also produces many meaningless slivers along the DEM edge; convert to polygons and drop anything under a project threshold such as 1 km².

Automate steps 1 to 7

ArcGIS Pro 3.x with Spatial Analyst; Python 3 with the bundled arcpy.

# Hydrology workflow automation, ArcGIS Pro 3.x
import arcpy
from arcpy.sa import *

arcpy.CheckOutExtension('Spatial')
arcpy.env.workspace = r'C:/Project/Workspace/hydrology.gdb'
arcpy.env.overwriteOutput = True

dem = 'DEM_raw'
threshold = 1000  # cells

dem_fill = Fill(dem);                                    dem_fill.save('DEM_fill')
flowdir  = FlowDirection(dem_fill, 'NORMAL', '', 'D8');  flowdir.save('FlowDir_D8')
flowacc  = FlowAccumulation(flowdir, '', 'FLOAT', 'D8'); flowacc.save('FlowAcc')
streams  = Con(flowacc >= threshold, 1);                 streams.save('Stream_raw')
link     = StreamLink(streams, flowdir);                 link.save('StreamLink')
order    = StreamOrder(streams, flowdir, 'STRAHLER');    order.save('StreamOrder_S')
arcpy.sa.StreamToFeature(order, flowdir, 'StreamNetwork', 'SIMPLIFY')

arcpy.CheckInExtension('Spatial')
print('Hydrology workflow complete.')

From outputs to morphometry

With stream counts and lengths per order, basin area A, perimeter P, basin length L_b and relief H, the standard Horton–Strahler parameters follow directly:

Map it properly

• DEM as a hillshade with a terrain colour ramp at about 50 percent transparency.
• Streams graduated by Strahler order, 0.3 to 1.5 mm, blue.
• Watershed boundary as outline only, 1.0 mm, black or dark red; outlets as red triangles.
• Title with basin, DEM year and metric; north arrow; metric scale bar; legend; inset of the basin within the state and India; graticule on all four sides; DEM source and resolution; EPSG code in the margin; author and date; 300 DPI export.

Troubleshooting

• Watershed is a handful of cells. Pour point not on the channel. Re-snap with 250 to 500 m and verify.
• Parallel streams on flat ground. D8 artefact. Switch to D-Infinity or MFD, or add ±0.01 m of noise to the DEM before Flow Direction.
• Fill changes nothing. Run Sink first; leave Z limit blank.
• Boundary cuts across a known basin. Pour point off the channel, or an edge effect. Re-snap, buffer the DEM extent by 10 percent, re-run.
• Network too dense or too sparse. Wrong threshold for the terrain. Test half and double, compare with SOI sheets.
• Outputs land in the wrong projection. Environment settings were not applied. Set output CRS and snap raster before every run.
• NoData holes in accumulation. Voids in the DEM. Fill them with Focal Statistics (mean) and re-run.
• Stream Order returns blank. Stream cells are not exactly 1. Re-run the Con expression.
• Very slow on a large DEM. Convert to LZW-compressed TIFF and process by sub-basin.
• Jagged watershed edges. Cell-wise boundary. Smooth Polygon, PAEK, tolerance of 1 to 2 cell sizes.

Frequently asked questions

Written by Dr. Aran Castro A J, PhD in Applied Geology, GIS Manager at Geospatial Campus, and author of three books on GIS and coastal science.

← All blog posts