Fixing OpenStreetMap Lakes and Parks That Show as Lines Instead of Polygons
On this page
Problem statement
You extract lakes, parks or buildings from OpenStreetMap and they come out as outlines that cannot be filled, have no area and cannot be used in a spatial join. Worse, the largest features, a wide river, a big park, a lake with islands, are missing altogether. Both symptoms have the same cause: the reader was not asked to assemble areas. Without area assembly, closed ways come back as rings of line, and multipolygon relations, which is how large and complex areas are mapped, are not built at all.
This page reproduces both symptoms with the water bodies of Washington, D.C. from the OpenStreetMap extract of 1 January 2026, measures what is lost, and fixes it with pyosmium's area assembly.
Quick answer
Enable area assembly and keep only the areas:
fp = (osmium.FileProcessor(str(PBF)).with_areas()
.with_filter(osmium.filter.EntityFilter(osmium.osm.AREA))
.with_filter(osmium.filter.TagFilter(("natural", "water")))
.with_filter(osmium.filter.GeoInterfaceFilter(tags=["natural", "water", "name"])))
good = gpd.GeoDataFrame.from_features(fp, crs=4326)
print(f"water areas: {len(good):,}; {good.geom_type.value_counts().to_dict()}; {good.to_crs(32618).area.sum() / 1e6:.2f} km²")
water areas: 325; {'MultiPolygon': 325}; 23.88 km²
With area assembly, the extract yields 325 water polygons covering 23.88 km², including the parts of the extract outside the District.
Read without areas, see what went missing, why rings are not enough
1. Reproduce: read without areas
fp = (osmium.FileProcessor(str(PBF)).with_locations()
.with_filter(osmium.filter.TagFilter(("natural", "water")))
.with_filter(osmium.filter.GeoInterfaceFilter(tags=["natural", "water", "name"])))
bad = gpd.GeoDataFrame.from_features(fp, crs=4326)
print(f"natural=water without areas: {len(bad):,} features; {bad.geom_type.value_counts().to_dict()}")
print(f"area of those geometries: {bad.to_crs(32618).area.sum():.0f} m²")
natural=water without areas: 303 features; {'LineString': 303}
area of those geometries: 0 m²
Without it, the same filter returns 303 lines with a total area of zero square metres.
2. What was missing entirely
Multipolygon relations are not built without area assembly. Count the water relations in the file:
con = duckdb.connect(); con.sql("INSTALL spatial; LOAD spatial")
rel = con.sql(f"""SELECT id, tags['name'] AS name, tags['water'] AS water, len(refs) AS members FROM ST_ReadOSM('{PBF}')
WHERE kind = 'relation' AND tags['natural'] = 'water' ORDER BY members DESC""").df()
print(f"natural=water multipolygon relations: {len(rel)}; largest by members:")
print(rel.head(5).to_string(index=False))
natural=water multipolygon relations: 22; largest by members:
id name water members
1857635 NaN river 15
1858408 NaN river 13
9370920 NaN river 12
1870168 Kingman Lake reservoir 7
5991610 Washington Channel river 6
The file has 22 water multipolygon relations, the largest of them unnamed sections of the Potomac and Anacostia surfaces. Kingman Lake and the Washington Channel are among the named ones. None of them is built without area assembly.
3. Turning rings into polygons yourself is not enough
Converting each closed line to a polygon recovers simple ponds, but not the relations, and it fills islands:
from shapely.geometry import Polygon
rings = bad[bad.geometry.is_closed].geometry.apply(lambda g: Polygon(g.coords))
print(f"closed rings converted: {len(rings):,}, {gpd.GeoSeries(rings, crs=4326).to_crs(32618).area.sum() / 1e6:.2f} km², against {good.to_crs(32618).area.sum() / 1e6:.2f} km² with area assembly")
closed rings converted: 303, 1.74 km², against 23.88 km² with area assembly
Converting the 303 rings to polygons recovers only 1.74 km², less than a tenth of the 23.88 km² that area assembly finds, because the rivers are relations.
Parks, the same problem
The data used on this page
The Geofabrik OpenStreetMap extract for the District of Columbia dated 1 January 2026 (© OpenStreetMap contributors, ODbL), read with pyosmium 4 and DuckDB's spatial extension.
from pathlib import Path
import duckdb
import geopandas as gpd
import osmium
import requests
Path("data").mkdir(exist_ok=True)
PBF = Path("data/district-of-columbia-260101.osm.pbf")
if not PBF.exists():
r = requests.get("https://download.geofabrik.de/north-america/us/district-of-columbia-260101.osm.pbf", timeout=600)
r.raise_for_status()
PBF.write_bytes(r.content)
print(f"{PBF.name} ready")
district-of-columbia-260101.osm.pbf ready
The dated D.C. extract is ready.
Example 1: parks, the same problem
Parks are mapped the same way. Compare leisure=park with and without area assembly:
def read(areas):
fp = osmium.FileProcessor(str(PBF))
fp = fp.with_areas().with_filter(osmium.filter.EntityFilter(osmium.osm.NODE | osmium.osm.AREA)) if areas else fp.with_locations()
fp = fp.with_filter(osmium.filter.TagFilter(("leisure", "park"))).with_filter(osmium.filter.GeoInterfaceFilter(tags=["name"]))
return gpd.GeoDataFrame.from_features(fp, crs=4326)
for label, g in (("without areas", read(False)), ("with areas", read(True))):
print(f"{label:>13}: {len(g):,} parks; {g.geom_type.value_counts().to_dict()}; polygon area {g.to_crs(32618).area.sum() / 1e6:.2f} km²")
without areas: 730 parks; {'LineString': 691, 'Point': 39}; polygon area 0.00 km²
with areas: 775 parks; {'MultiPolygon': 736, 'Point': 39}; polygon area 24.73 km²
Parks behave the same: without areas there are 730 features, 691 of them lines, with no area at all. With areas there are 775, of which 736 are polygons covering 24.73 km².
Closed ways and multipolygon relations
OpenStreetMap stores areas in two ways. A simple area is a closed way: a ring of nodes with tags such as natural=water or leisure=park. A complex area is a multipolygon relation: a relation whose members are ways with the roles outer and inner, assembled into polygons with holes. Large features, such as a river surface hundreds of kilometres long, a lake with islands, or a park made of many parcels, are almost always relations, because a single way has a limit of 2,000 nodes and cannot have holes.
A reader that only processes ways sees the closed way as a sequence of points, so it returns a line. It never assembles relations, so multipolygons disappear, and with them the biggest features. Converting rings to polygons afterwards fixes the simple cases but cannot recover relations and ignores inner rings. The only correct fix is to let the reader assemble areas, with_areas() in pyosmium, and to keep the areas rather than the ways. Then every lake and park arrives as a polygon, holes and all.
Incomplete relations, coastlines and invalid areas
- Incomplete relations: a relation whose members run outside the extract cannot be assembled and is silently dropped; see polygons missing at extract edges.
- Coastlines: seas are not
natural=water; coastline polygons are built separately fromnatural=coastlinelines. - Invalid areas: broken rings in the data fail assembly; pyosmium skips them, so compare counts with the raw relations.
- Duplicates: with areas enabled, also drop the way copies; see features counted twice.
Internal links
- How to get rivers and water bodies from OpenStreetMap in Python — water in full
- Fixing OSM features counted twice — the other side of
with_areas() - The OpenStreetMap data model explained — ways and relations
FAQ
Why are OSM lakes lines instead of polygons?
The reader returned closed ways as lines because area assembly was off; enable it (with_areas() in pyosmium) and keep the areas.
Why are large lakes and parks missing from my OSM data?
They are multipolygon relations, which are only built when areas are assembled.
Can I just convert closed lines to polygons?
Only for simple areas; relations stay missing and inner holes are filled.
What is a multipolygon relation?
A relation of outer and inner member ways that together describe one area, possibly with holes.