🔧 Fix Urban form & morphology Census & demographics Urban Morphology · step 18 Notebook By · Updated

Density Depends on the Boundary

On this page

Problem statement

The population density you computed for Den Haag does not match the published one. The building coverage of a district is half what a planning report says. Two cities rank in the opposite order from last year's study. The numerator — people, footprint area, floor area — is usually the same in both calculations. The denominator is not: the land area the count is divided by depends on the boundary, and boundaries include water, parks, streets and farmland in very different proportions.

Density is a ratio with an argument hidden in its denominator. Fixing a mismatch means finding which land each figure divided by.

Quick answer

List the candidate denominators and compute the density on each:

import io, zipfile, requests, numpy as np, rasterio, rasterio.mask, geopandas as gpd, osmnx as ox

name = "GHS_POP_E2020_GLOBE_R2023A_4326_30ss"
zipfile.ZipFile(io.BytesIO(requests.get(
    "https://jeodpp.jrc.ec.europa.eu/ftp/jrc-opendata/GHSL/GHS_POP_GLOBE_R2023A/"
    f"{name}/V1-0/tiles/{name}_V1_0_R4_C19.zip", timeout=600).content)).extractall("ghsl_nl")
src = rasterio.open(f"ghsl_nl/{name}_V1_0_R4_C19.tif")

ox.settings.http_user_agent = "spatialworkflow-examples/1.0"
ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"   # mirror; overpass-api.de refused
ox.settings.overpass_rate_limit = False
gem = ox.geocode_to_gdf("Den Haag, Netherlands").to_crs(4326)
people = np.clip(rasterio.mask.mask(src, gem.geometry, crop=True, nodata=0)[0], 0, None).sum()

land = gpd.read_file("https://naciscdn.org/naturalearth/10m/physical/ne_10m_land.zip").to_crs(4326)
water = ox.features_from_polygon(gem.geometry.iloc[0], {"natural": "water"})
water = water[water.geom_type.isin(["Polygon", "MultiPolygon"])]
parks = ox.features_from_polygon(gem.geometry.iloc[0], {"leisure": "park", "natural": "wood"})
parks = parks[parks.geom_type.isin(["Polygon", "MultiPolygon"])]

km2 = lambda g: gpd.GeoSeries([g], crs=4326).to_crs(28992).area.iloc[0] / 1e6
poly = gem.geometry.iloc[0]
on_land = poly.intersection(land.union_all())
no_water = on_land.difference(water.union_all())
no_green = no_water.difference(parks.union_all())
print(f"population (GHSL 2020) inside the municipal polygon: {people:,.0f}")
for label, g in (("municipal polygon", poly), ("… on land", on_land),
                 ("… minus inland water", no_water), ("… minus parks and woods", no_green)):
    print(f"{label:<26} {km2(g):6.1f} km² → {people / km2(g):7,.0f} people/km²")
population (GHSL 2020) inside the municipal polygon: 550,960
municipal polygon            98.1 km² →   5,615 people/km²
… on land                    81.0 km² →   6,806 people/km²
… minus inland water         76.3 km² →   7,223 people/km²
… minus parks and woods      69.9 km² →   7,888 people/km²

The same people, four denominators, and densities from 5,615 to 7,888 per km² — a 40% spread, far more than any year-to-year change. None is wrong; each is a different definition of "the land".

Write down both parts, check water, decide gross or net

1. Write down both parts of the published figure

Find the numerator and, above all, the area the source divided by. Statistical offices usually state it: total area, land area, or net residential area.

2. Check for water

Coastal and riverside municipalities include water. Land area and total area can differ by a large fraction.

3. Decide gross or net

Gross density divides by all land; net density by residential land only, excluding parks, roads and non-residential uses. Planning documents often use net; statistical yearbooks gross.

4. Check the boundary itself

A geocoded polygon may be the municipality, a city district, or a place outline. Print its name and area before trusting it.

5. Recompute with the stated denominator

Once both figures use the same land, differences that remain are data differences — census against GHSL, dates, and so on.

Density mismatches with the likely denominator difference and the fix for each. a different denominator, not a different citymismatchcause → fixyours lower on a coastpolygon includes sea → clip tolandplan figure much highernet density → remove parks,streetswildly different areawrong geocode → print name and km²same area, still offdifferent data → compare sourcesDen Haag: 5,615 to 7,888 people/km² for the same people spatialworkflow.io
Most density mismatches are a different denominator, not a different city.

Coverage on three land units

The data used on this page

Building footprints with LiDAR heights from the 3D BAG (the Dutch building register with heights from AHN), fetched from its WFS for central Den Haag — 4.285–4.330° E, 52.065–52.095° N — exactly as in How to download building footprints from OSM and Overture, and saved as denhaag_bag.parquet in the Dutch national grid (EPSG:28992). The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import zipfile

import geopandas as gpd
import requests
from pyproj import Transformer

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-density-depends-on-boundary-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

if not os.path.exists("denhaag_bag.parquet"):
    W, S, E, N = 4.285, 52.065, 4.330, 52.095
    to_rd = Transformer.from_crs(4326, 28992, always_xy=True)
    (x0, y0), (x1, y1) = to_rd.transform(W, S), to_rd.transform(E, N)
    feats, start = [], 0
    while True:
        fc = requests.get("https://data.3dbag.nl/api/BAG3D/wfs", params={
            "service": "WFS", "version": "2.0.0", "request": "GetFeature", "typeNames": "BAG3D:lod12",
            "count": 5000, "startIndex": start, "outputFormat": "application/json",
            "bbox": f"{x0},{y0},{x1},{y1},EPSG:28992"}, timeout=300).json()
        feats += fc["features"]
        if len(fc["features"]) < 5000:
            break
        start += 5000
    bag = gpd.GeoDataFrame.from_features(feats, crs=28992)
    bag["height"] = bag.b3_h_70p - bag.b3_h_maaiveld
    bag[["identificatie", "b3_h_70p", "b3_h_maaiveld", "height", "b3_opp_grond", "geometry"]] \
        .to_parquet("denhaag_bag.parquet")
print(len(gpd.read_parquet("denhaag_bag.parquet")), "buildings")
23929 buildings

Example 1 — building coverage on three land units

b = gpd.read_parquet("denhaag_bag.parquet")                   # central Den Haag footprints
extent = gpd.GeoSeries([b.union_all().convex_hull], crs=28992)
bbox_area = (b.total_bounds[2] - b.total_bounds[0]) * (b.total_bounds[3] - b.total_bounds[1])
blocks_land = extent.iloc[0].difference(
    ox.graph_to_gdfs(ox.graph_from_point((52.0799, 4.3113), dist=2000, network_type="drive"),
                     nodes=False).to_crs(28992).buffer(6).union_all())
for label, area in (("bounding box", bbox_area), ("convex hull of buildings", extent.area.iloc[0]),
                    ("hull minus 12 m street corridors", blocks_land.area)):
    print(f"{label:<34} coverage {b.area.sum() / area:5.1%}")
bounding box                       coverage 28.8%
convex hull of buildings           coverage 29.9%
hull minus 12 m street corridors   coverage 37.4%

Removing the streets alone — a 12 m corridor along each street centreline — raises coverage from 30% to 37%. Plans that quote "net" coverage typically exclude streets; maps built on grid squares include them.

Example 2 — the polygon you actually got

for q in ("Den Haag, Netherlands", "Centrum, Den Haag, Netherlands"):
    g = ox.geocode_to_gdf(q)
    print(f"{q:<32} → {g.display_name.iloc[0][:60]:<60} {g.to_crs(28992).area.iloc[0] / 1e6:6.1f} km²")
Den Haag, Netherlands            → The Hague, South Holland, Netherlands                          98.1 km²
Centrum, Den Haag, Netherlands   → Centrum, The Hague, South Holland, Netherlands                  7.8 km²

Two queries, two very different polygons. A density computed on whichever one came back is a density of that polygon.

Example 3 — which denominator a source used

reported = {"total area incl. water": km2(poly), "land area": km2(on_land)}
for label, area in reported.items():
    print(f"if the source divided by {label} ({area:.1f} km²): {people / area:,.0f} people/km²")
print("compare with the source's own density figure to see which it used")
if the source divided by total area incl. water (98.1 km²): 5,615 people/km²
if the source divided by land area (81.0 km²): 6,806 people/km²
compare with the source's own density figure to see which it used

Reverse-engineering the denominator from a published density and a known population is often the fastest way to find the mismatch.

Population density of Den Haag computed on four different land areas, GHS-POP 2020. the same 550,960 people, people per km²polygon (98.1 km²)5,615on land (81.0)6,806− inland water (76.3)7,223− parks, woods (69.9)7,888040008000A 40% spread from the choice of denominator alone spatialworkflow.io
One population, four land areas, four densities.

Why the denominator varies so much

Why the denominator varies so much

Cities contain large amounts of land that no one lives on: water, parks, rail yards, industrial estates, beaches, roads. Whether those are in the denominator is a definitional choice, and in a coastal city with dunes and harbours they can be a large share of the total.

Why gross and net both have uses

Gross density describes land consumption — how much land the city uses per person. Net density describes the living environment — how closely people live where they live. Policy debates often compare one with the other by accident.

Why the boundary is part of the data

A polygon from a geocoder is an answer to a search query. It may include the sea to the territorial limit, exclude an enclave, or represent a different administrative level than intended. Treat it as data to be checked, not as a given.

Building coverage of central Den Haag on three land denominators. the same footprints, three coveragesland countedcoveragebounding box28.8%convex hull of the buildings29.9%hull minus 12 m street corridors37.4%Net coverage excludes streets; grid-based coverage includes them spatialworkflow.io
Taking out the streets adds seven points of coverage without changing a building.

Coastlines, population sources and street corridors

  • Coastline datasets differ in resolution; a coarse land mask misallocates the coastal strip.
  • Population sources (census, register, GHSL) differ by a few percent; do not mix them in one comparison.
  • Grid-based densities avoid boundary choices inside the grid but inherit them at the edge.
  • Street corridor widths are approximations when only centrelines are available.
  • Floor area density (FAR) has the same denominator problem; see the coverage guide.

FAQ

Why doesn't my density match the official figure?

The official figure probably divides by a different area — land only, or net residential land.

What is the difference between gross and net density?

Gross divides by all land in the unit; net by residential land only.

Should water be excluded?

For land-use and living-density questions, yes. Say so in the figure's definition.

How do I know what polygon the geocoder returned?

Print its name and area; compare with the administrative level you meant.

Is a grid-based density free of boundary problems?

Inside the grid, largely; at the edge of the study area and in cells that include water, no.

Does the same problem affect building coverage?

Yes — coverage on blocks, on grids and on whole districts differs because of streets and open space.