If I have the polygon defining boundaries of a river network, is there any way to directly segment the tributaries out using geospatial data analysis?

To elaborate, I am looking for a workflow more than a specific software, but a code in python or R would be preferred. Shapely has potential but I really don't want a commercial algorithm/tool i.e. Esri products if that makes sense. QGIS is fine but it is a high dependency for what I have in mind.

enter image description here

  • 1
    Add any gis package - just one - you intend to use. No software specification is a reason of interesting questions being closed. +1
    – FelixIP
    Apr 12 at 0:24
  • BTW do you have stream network itself? Does it run more or less in the middle of flood plain? If yes, solution is simple
    – FelixIP
    Apr 12 at 0:39
  • No, stream networks are DEM driven (additional data requirement) and very inaccurate (centerline or even thalweg doesn't mean much in general case). But it's possible to do Medial Axis Skeletonization or Straight Skeletonization to derive pseudo stream networks. How would that help? I understand we can trace the graph topology and create junctions etc, but how can that be used to recover the main polygon segments?
    – MathX
    Apr 12 at 0:58

3 Answers 3


Good news you are confident with skeleton creation, good luck. Assuming that you have one:

  • Calculate distance of your nodes to outlet
  • Flip segments if Travel[from]<Travel[to]
  • Path from remotest source node is your 'main' stream

enter image description here

Triangulate vertices of flood plain, extract edges that intersect your streams. Delete ones that intersect main stream:

enter image description here

Use remaining sections and flood plain polygon to built multiple polygons. The ones that don't intersect 'main' stream are tributaries flooding areas.

enter image description here

For flood plain with islands, you need all stream branches downstream from remotest node. In terms of graphs we are talking all simple paths to outlet, and possibly OMG with parallel edges.

  • Nice, thanks! I was thinking about triangulation but didn't know how to use them. I guess I can find the triangles that intersect any stream I want and just merge them? The edge selection seems a little more complicated.
    – MathX
    Apr 12 at 2:25
  • Selection of triangles only is not going to work, there are potentially hundreds of tiny ones along flood plain banks, that are not going to be selected. BTW note that sections shown (and hidden) are perfect candidates for cross-sections, with some weeding of course. The most challenging thing in skeleton for me was finding 'source' points, no ideal algorithm so far:(
    – FelixIP
    Apr 12 at 2:29

The following algorithm is just based on your test example you asked input for in the comment... so it will need iterations to improve for more/more complex situations, but it shows the main idea I would follow to split the main river from its tributaries.

"POLYGON((0 0,0 8,-4 8,-4 10,0 10,0 14,2 14,2 9,4 9,4 8,2 8,2 2,10 2,10 0,2 0,2 -5,0 -5,0 0))"


  • determine the centerline of the river polygon
  • determine the main centerline of the river: keep the main path + the longest start and end branch.
  • triangulate the river polygon
  • determine the tributary triangles: all triangles that don't intersect with the main centerline
  • determine the tributary polygons: dissolve/union the tributary triangles
  • cleanup the "tributary polygons": with more complex examples, I think some (typically smaller) pieces will actually belong to the main river. Some heuristics will be needed to find (most of) these cases.
  • Apply difference of the "tributary polygons" on the main river polygon, to only retain the main river.

The sample script doesn't implement everything in a generic way. Mainly for determining the main river from the centerline, I just use a variant in the pygeoops library that still tries to retain all branches rather than only the start and end branch, but this would just be a minor change in the current code.

Sample script:

from matplotlib import pyplot as plt
import pygeoops
from pygeoops._centerline import _remove_short_branches
import shapely
import shapely.plotting as plotter

# Sample data.
wkt = "POLYGON((0 0,0 8,-4 8,-4 10,0 10,0 14,2 14,2 9,4 9,4 8,2 8,2 2,10 2,10 0,2 0,2 -5,0 -5,0 0))"
river = shapely.from_wkt(wkt)

# Calculate centerline.
river_center = pygeoops.centerline(river)
# Determine main river. Not the ideal algorith, but it's only a POC, and it shows the idea.
main_river_center = _remove_short_branches(
    river_center, min_branch_length=5, remove_one_by_one=True

# Split river to delaunay triangles + clip to the river polygon.
river_triangles = shapely.delaunay_triangles(river)
river_triangles_clip = pygeoops.collection_extract(
    shapely.intersection(shapely.get_parts(river_triangles), river),
river_triangles_clip = [poly for poly in river_triangles_clip if poly is not None]

# Form the tributaries from the triangles not intersecting the main river centerline.
river_triangles_clip = [
    for poly in river_triangles_clip
    if not shapely.intersects(poly, main_river_center)
tributaries = shapely.get_parts(shapely.union_all(river_triangles_clip))

# Remove the tributaries from the river.
main_river = pygeoops.difference_all(river, tributaries, keep_geom_type=True)

# Plot
plotter.plot_polygon(main_river, color="green")



enter image description here

  • 1
    It's never fully automatic. Length limit has a tendency to remove last edge before outlet if there are tributaries nearby
    – FelixIP
    Apr 12 at 6:31
  • 1
    Thanks for the suggestion but creating centerline is not an issue. I was not looking for ways to get a centerline, I already have my own code for that. The questions was partitioning a polygon to main stem and branches.
    – MathX
    Apr 12 at 15:15
  • Lets say you use this method to get a centerline, how would you partition (0 0,0 8,-4 8,-4 10,0 10,0 14,2 14,2 9,4 9,4 8,2 8,2 2,10 2,10 0,2 0,2 -5,0 -5,0 0) to the main stem and side branches automatically?
    – MathX
    Apr 12 at 16:21
  • 1
    This is very nice, thnks! I will be experimenting with it definitely.
    – MathX
    Apr 15 at 19:58

Here is another possible approach you can explore. You say you have the centreline or can at least build it. The segmentation you seem to want appears to be what I would call Hack Order. RivEX can assign Hack order to a network and you can take advantage of this.

1] Here is a network colour coded by Hack order.

Hack order

2] Network along with polygon bank edge data displaying.


3] Run the Euclidean Allocation tool to build a raster where the pixel value is the nearest Hack order. Set cell size to be quite small (I had used 5m) and ensure the river bank dataset is a mask, the result is:


4] Convert raster back into polygons.

5] Depending upon your usage you might accept output as-is but if that is not acceptable there is some minor clean up editing required at some junctions. One useful "side affect" of this approach is that if your centreline falls outside the bank polygon data the mask constrains the output and forces the results to be within the bank data.

  • Thanks for the suggestion. Stream ordering regardless of Horton, Hack or others have some merits but the biggest issue is having a raster output. Back and forth rasterizing and polygonization is usually not efficient or very desirable. It will depend on the resolution of the raster and pixelated output that may not be exactly the same coordinates as the original polygon, unless a very high resolution raster cell size is employed. I am hoping more vectorized approaches are out there.
    – MathX
    Apr 15 at 20:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.