Repository for Ideas & Research
Open Source GIS, Hydrologic Modeling, Optimization

How to calculate the longest flow path in GRASS GIS can calculate the flow length downstream (FLDS) and flow length upstream (FLUS) raster maps. In the FLDS map, each cell has the flow length from that cell to the downstream-most outlet cell flowing outside the watershed. Outlet cells along the watershed boundary get assigned 0. Similarly, in the FLUS map, each cell has the flow length from that cell to the upstream-most cell where the watershed divide starts. Those cells along the watershed divide get assigned 0. These two raster maps have the same maximum value at the outlet cell in FLUS and the upstream-most cell on the longest flow path (LFP) in FLDS. If we add both maps, cells on the output map will have the same maximum value greater than non-LFP cells and we use that fact to calculate LFP.

This tutorial uses the North Carolina sample data set.

1   Calculate the longest flow path raster map

First, create the drainage direction map (drain_directions) from the elevation map (elevation).

g.region rast=elevation
r.watershed elevation=elevation drainage=drain_directions

Create the basin map (basin) at x=642455 and y=222614 and apply it as mask.

r.water.outlet input=drain_directions output=basin coordinates=642455,222614
g.region rast=basin
r.mask raster=basin

Create a raster map (outlet) containing the outlet point.

echo 642455,222614 | input=- output=outlet separator=, input=outlet output=outlet use=cat

Calculate FLDS (flds) and FLUS (flus). -o stream_rast=outlet direction=drain_directions method=downstream distance=flds -o stream_rast=outlet direction=drain_directions method=upstream distance=flus

Combine FLDS and FLUS to create a new raster map (fldsus).

r.mapcalc expression="fldsus=flds+flus"

Find the maximum flow length that each cell on FLDSUS got assigned. -r map=fldsus | sed '/^max=/!d; s/^max=//'

Create the LFP raster map (lfp) allowing small floating-point errors in FLDSUS. The maximum cell value obtained from the last step (or LFP distance) was 21452.825. Subtract a small number to avoid a potential floating-point error.

r.mapcalc expression="lfp=if(fldsus>=21452.825-0.0005,1,null())"

Now, lfp is the longest flow path raster map.

2   Calculate the longest flow path vector map

Simply converting the lfp raster map to vector won’t produce the correct vector map because of this issue. To properly trace the longest flow path in vector format, we need to find the headwater raster cells and start tracing flow directions from there. Using the maximum flow length found above (21452.825), find the headwater cells from FLDS and create a new vector map with the headwater points.

r.mapcalc expression="heads=if(flds>=21452.825-0.0005,1,null())" input=heads output=heads type=point

Trace flow directions from the headwater points to calculate the longest flow path in vector format.

r.path input=drain_directions vector_path=path start_points=heads

The output vector map (path) includes the longest flow path, but it passes through the outlet and flows beyond the watershed. We need to split the path at the outlet first and select only the upstream segment within the watershed. To split the path at the outlet, create a new vector map with the snapped outlet point. input=outlet output=outlet type=point

Find the coordinates of the snapped outlet. Read snapped x=642455 and y=222615 and split the path at this point. -p map=outlet option=coor
v.edit map=path tool=break coords=642455,222615

Select the upstream segment only that touches the headwater points. ainput=path binput=heads output=lfp

Finally, we got the longest flow path in lfp.


This example produces two lines with the same flow length.


3   r.lfp addon

I have implemented the procedure above as a GRASS addon module r.lfp.