How to Get Rivers and Water Bodies From OpenStreetMap in Python

On this page

Problem statement

OpenStreetMap maps water two ways at once. Rivers and streams are lines, waterway=river or stream, drawn along the centre and pointing downstream. The water surface itself is a polygon, natural=water, often a large multipolygon relation for a wide river or a lake. A map that uses only the lines shows the Potomac as a thin thread. A map that uses only the polygons loses every small stream too narrow to have a surface mapped.

This page extracts both from the OpenStreetMap extract of Washington, D.C. dated 1 January 2026. It measures waterway length by type and by named river, assembles the water polygons, and compares OSM's water area inside the District with the US Census Bureau's own figure.

Quick answer

Read the waterway lines and the water polygons separately:

fp = (osmium.FileProcessor(str(PBF)).with_locations().with_filter(osmium.filter.EntityFilter(osmium.osm.WAY))
      .with_filter(osmium.filter.KeyFilter("waterway"))
      .with_filter(osmium.filter.GeoInterfaceFilter(tags=["waterway", "name", "tunnel"])))
W = gpd.GeoDataFrame.from_features(fp, crs=4326)
W = gpd.clip(W[W.geom_type == "LineString"], DC).to_crs(32618)
W["km"] = W.length / 1000
fp = (osmium.FileProcessor(str(PBF)).with_areas().with_filter(osmium.filter.TagFilter(("natural", "water")))
      .with_filter(osmium.filter.GeoInterfaceFilter(tags=["water", "name"])))
A = gpd.GeoDataFrame.from_features(fp, crs=4326)
A = gpd.clip(A[A.geom_type.isin(["Polygon", "MultiPolygon"])], DC).to_crs(32618)
print(f"waterway lines: {len(W):,}, {W.km.sum():.0f} km; water polygons: {len(A):,}, {A.area.sum() / 1e6:.2f} km²")

waterway lines: 312, 136 km; water polygons: 319, 21.07 km²

Inside the District there are 312 waterway lines, 136 km in all, and 319 water polygons whose areas sum to 21.07 km² — more than the water they cover, as the Census comparison below shows.

Map of water in D.C.: pale blue natural=water polygons of the Potomac and Anacostia, dark river lines running through them, blue streams such as Rock Creek's tributaries, and short red culvert sections. — Lines inside polygons.

Waterways by type, buried streams, water area against the Census

1. Waterways by type and by river

print(W.groupby("waterway").km.sum().sort_values(ascending=False).round(1).to_dict())
rivers = W[W.name.notna()].groupby("name").km.sum().sort_values(ascending=False)
print(rivers.head(6).round(1).to_string())

{'stream': 60.4, 'river': 56.2, 'pressurised': 9.5, 'canal': 7.3, 'ditch': 0.4, 'link': 0.4, 'dock': 0.4, 'drain': 0.4, 'dam': 0.3, 'floating_barrier': 0.1, 'derelict_canal': 0.1, 'fish_pass': 0.1, 'weir': 0.0, 'lock_gate': 0.0, 'sluice_gate': 0.0}
name
Rock Creek                   15.6
Potomac River                13.6
Anacostia River              11.5
Chesapeake and Ohio Canal     7.1
Oxon Run                      7.0
Washington City Tunnel        6.4

Streams (60.4 km) and rivers (56.2 km) make up most of the length. The pressurised class, 9.5 km, is the water-supply tunnels, among them the Washington City Tunnel. By name, Rock Creek is the longest waterway inside the District, at 15.6 km, ahead of the stretches of the Potomac (13.6 km) and the Anacostia (11.5 km) within it.

Waterway length by type.

2. Streams that run underground

Many urban streams were buried in pipes. OSM marks those sections with tunnel=culvert or tunnel=yes:

piped = W.tunnel.notna() & (W.tunnel != "no")
print(f"waterway length in tunnels or culverts: {W.loc[piped, 'km'].sum():.1f} km of {W.km.sum():.0f} km ({W.loc[piped, 'km'].sum() / W.km.sum():.0%})")
print(W.loc[piped, "tunnel"].value_counts().to_dict())

waterway length in tunnels or culverts: 2.3 km of 136 km (2%)
{'culvert': 59, 'yes': 16, 'flooded': 1}

Only 2.3 km of the 136 km, 2%, is tagged as running underground, mostly through culverts. Streams buried long ago are often not mapped at all.

3. Water area against the Census

The Census cartographic boundaries carry each county's water area, AWATER. Compare it with the area of OSM's water polygons inside the boundary:

row = C[C.GEOID == "11001"].iloc[0]
osm_km2 = A.union_all().area / 1e6
print(f"Census AWATER: {row.AWATER / 1e6:.2f} km²; OSM natural=water inside D.C.: {osm_km2:.2f} km²")
print("largest OSM water polygons:", A.assign(km2=A.area / 1e6).nlargest(4, "km2")[["name", "water", "km2"]].round(2).to_dict("records"))

Census AWATER: 18.71 km²; OSM natural=water inside D.C.: 18.61 km²
largest OSM water polygons: [{'name': nan, 'water': 'river', 'km2': 9.77}, {'name': nan, 'water': 'river', 'km2': 3.25}, {'name': 'Anacostia River', 'water': 'river', 'km2': 2.45}, {'name': nan, 'water': 'river', 'km2': 2.45}]

Merged into one surface, OSM's water polygons cover 18.61 km² of the District, within 0.1 km² of the Census Bureau's 18.71 km². Summed one by one they come to 21.07 km², because two polygons of 2.45 km² cover the same water: an unnamed part of the river multipolygon and a separate polygon named Anacostia River. Area from water polygons has to be taken from their union, not their sum.

Which way the water flows

The data used on this page

The Geofabrik OpenStreetMap extract for the District of Columbia dated 1 January 2026 (© OpenStreetMap contributors, ODbL), and the US Census Bureau's 2023 cartographic boundary file for counties at 1:500,000, whose District outline includes its share of the Potomac and Anacostia.

from pathlib import Path

import geopandas as gpd
import osmium
import requests

Path("data").mkdir(exist_ok=True)
def cached(path, url):
    p = Path(path)
    if not p.exists():
        r = requests.get(url, timeout=600)
        r.raise_for_status()
        p.write_bytes(r.content)
    return p
PBF = cached("data/district-of-columbia-260101.osm.pbf", "https://download.geofabrik.de/north-america/us/district-of-columbia-260101.osm.pbf")
C = gpd.read_file(cached("data/cb_2023_us_county_500k.zip", "https://www2.census.gov/geo/tiger/GENZ2023/shp/cb_2023_us_county_500k.zip"))
DC = C[C.GEOID == "11001"].to_crs(4326)
print(f"D.C. boundary (1:500k): {DC.to_crs(32618).area.iloc[0] / 1e6:.1f} km²")

D.C. boundary (1:500k): 176.9 km²

At 1:500,000, the District covers 176.9 km², water included.

Example 1: which way does the water flow?

A waterway's direction is the direction of flow. Check that the Anacostia's ways run downstream, toward its mouth at the Potomac in the south-west:

ana = W[W.name == "Anacostia River"]
starts = ana.geometry.apply(lambda g: g.coords[0]); ends = ana.geometry.apply(lambda g: g.coords[-1])
south_west = sum((e[0] - s[0] < 0) or (e[1] - s[1] < 0) for s, e in zip(starts, ends))
print(f"Anacostia River ways: {len(ana)}; running toward the south or west: {south_west}")

Anacostia River ways: 8; running toward the south or west: 8

All eight Anacostia River ways run toward the south or west, downstream toward the Potomac.

Map of Rock Creek winding north to south through D.C., 15.6 km long, with the small streams that join it from both sides. — Rock Creek and its streams.

Flowing networks and water surfaces

Two ways OSM maps water.

OpenStreetMap separates the network of flowing water from its surface. waterway=* lines form the network: river for large rivers, stream for small ones, canal, ditch and drain for artificial channels. Each is drawn in the direction of flow, so the lines can be followed downstream like a drainage network. The surface is mapped with natural=water, refined by water=river, lake, pond or reservoir. Rivers wider than a few metres usually get both, a centre line inside a polygon of the riverbed. An older tagging, waterway=riverbank, survives in some places.

Comparing OSM's water with an official figure shows what each source measures. The Census water area comes from its own boundary and shoreline work. OSM's polygons follow whatever mappers drew, and in Washington the Potomac's surface is a multipolygon covering both banks. The two agree when both follow the same shore, and differ where piers, islands and tidal areas are drawn differently, or where small ponds are missing from one of them.

Large relations, intermittent streams and coastlines

Before mapping water.

Large water relations: multipolygons that cross an extract's border can be incomplete; see missing polygons at extract edges.
Intermittent streams: intermittent=yes marks streams that run only part of the year.
Coastlines: the sea is not natural=water; coastlines are natural=coastline lines, processed into land and water polygons at osmdata.openstreetmap.de.
Hydrology: OSM waterways are a map, not a hydrological network; check connectivity before routing flow along them.

Internal links

How to read an OpenStreetMap PBF file in Python — reading the extract
Fixing OSM lakes that show as lines — water polygons that come out as outlines
How to delineate a watershed from a pour point in Python — water from a DEM instead

FAQ

How do I get rivers from OpenStreetMap in Python?

Filter ways with the waterway key for the river lines, and natural=water with areas enabled for the water surfaces.

What is the difference between waterway=river and natural=water?

waterway=river is the centre line of the river, drawn downstream; natural=water with water=river is the polygon of its surface.

Why is my river only a thin line?

You used only the waterway lines; add the natural=water polygons for the river's width.

Do OSM waterways point downstream?

Yes, by convention each waterway way is drawn in the direction of flow.