How to Download Administrative Boundaries in Python

Problem statement

Almost every analysis needs a boundary layer — countries, regions, municipalities, wards. It is the layer you join to, aggregate by, and clip against. And it is the layer most likely to be quietly wrong, because the standard way to get one stopped working:

import geopandas as gpd

try:
    world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
except AttributeError as err:
    print(f"AttributeError: {err}")

AttributeError: The geopandas.dataset has been deprecated and was removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.

That line appears in thousands of tutorials and notebooks. It has not worked since GeoPandas 1.0, and the replacements are not obvious.

Worse, "administrative boundary" is not one thing. A country's ADM2 might be a county, a district or a municipality, and two providers will disagree about both the level and the geometry — including whose territory is whose.

Quick answer

For a programmatic, licence-clear source with global coverage, use the geoBoundaries API:

import geopandas as gpd
import requests

HEADERS = {"User-Agent": "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"}
r = requests.get("https://www.geoboundaries.org/api/current/gbOpen/NLD/ADM2/", headers=HEADERS, timeout=90)
r.raise_for_status()
meta = r.json()

boundaries = gpd.read_file(meta["gjDownloadURL"])
print(len(boundaries), meta["boundaryLicense"])
print(boundaries.columns.tolist())

344 CC0 1.0 Universal (CC0 1.0) Public Domain Dedication
['shapeName', 'shapeISO', 'shapeID', 'shapeGroup', 'shapeType', 'geometry']

The four practical routes, in the order you should try them:

Authority or coverage, ADM levels, metadata

1. Decide whether you need authority or coverage

These pull in opposite directions.

Authority means the source has the standing to define the boundary — a national mapping agency or statistics office. Use it whenever a number will be reported against an official area, because your municipality totals must match the ones the municipality publishes.

Coverage means one schema across many countries. A cross-border study cannot stitch together twelve national portals with twelve schemas, twelve CRSs and twelve update cycles.

You rarely get both. Pick deliberately and write down which you picked.

2. Understand what ADM0, ADM1, ADM2 mean

The ADMn levels are a comparative framing, not a definition:

• ADM0 — the country
• ADM1 — first subdivision (state, province, region, landen)
• ADM2 — second subdivision (county, district, municipality)

What ADM2 is varies by country: a Dutch ADM2 is a municipality (344 of them), a US ADM2 is a county (about 3,143). Never compare "ADM2 units" across countries as if they were the same kind of thing — that is the modifiable areal unit problem wearing a uniform.

Not every country has every level:

for level in ["ADM0", "ADM1", "ADM2", "ADM3"]:
    r = requests.get(f"https://www.geoboundaries.org/api/current/gbOpen/NLD/{level}/",
                     headers=HEADERS, timeout=60)
    print(level, r.status_code)

ADM0 200
ADM1 200
ADM2 200
ADM3 404

A 404 is the answer "this level does not exist here", delivered as an HTTP error. Handle it, or a loop over twenty countries dies on the first one with fewer levels.

3. Read the metadata before the geometry

The geoBoundaries API returns a metadata document, and the download URL is one field in it. The other fields are the ones that matter later:

for key in ["boundaryName", "boundaryType", "boundaryLicense",
            "boundarySource", "sourceDataUpdateDate"]:
    print(f"{key:22} {meta[key]}")

boundaryName           Netherlands
boundaryType           ADM2
boundaryLicense        CC0 1.0 Universal (CC0 1.0) Public Domain Dedication
boundarySource         National Georegister
sourceDataUpdateDate   Thu Jan 19 07:31:04 2023

sourceDataUpdateDate is the important one. This boundary set reflects the situation in early 2023 — fine for most work, wrong if a municipality has merged since. Record it with the data, exactly as in recording provenance at fetch time.

4. Check the identifiers before you join

The whole point of a boundary layer is joining statistics to it, and joins fail on codes:

print(boundaries[["shapeName", "shapeISO", "shapeID"]].head(3).to_string(index=False))
print("duplicate names:", boundaries["shapeName"].duplicated().sum())

  shapeName shapeISO                shapeID
  Groningen          6326949B19731519413910
     Almere          6326949B15537896627734
Stadskanaal          6326949B23876978853895
duplicate names: 0

shapeID is a geoBoundaries identifier, stable within a release and not the code your statistics use; shapeISO here is the province code, shared by dozens of municipalities. Joining on shapeName is all that is left, and the explanation below shows how that goes against the official statistics. Where an official code exists, get the boundary from the source that carries it.

5. Reproject before measuring anything

print(boundaries.crs)
areas = boundaries.to_crs("EPSG:28992")
areas["km2"] = areas.area / 1e6
print(areas.nlargest(3, "km2")[["shapeName", "km2"]].round(1).to_string(index=False))

EPSG:4326
      shapeName   km2
Súdwest-Fryslân 907.9
   Het Hogeland 903.1
       Lelystad 765.5

Global boundary sets are always EPSG:4326. .area on that is square degrees, which is not a unit of anything. See calculating area and distance correctly.

A boundary fetcher with metadata and caching

Example 1 — a boundary fetcher with metadata and caching

import json
from pathlib import Path

import geopandas as gpd
import requests

CACHE = Path("data/boundaries")


def geoboundaries(iso3, level="ADM1", *, release="gbOpen", refresh=False):
    """Fetch one country's boundaries, caching both the geometry and its provenance."""
    CACHE.mkdir(parents=True, exist_ok=True)
    gpkg = CACHE / f"{iso3}_{level}.gpkg"
    side = CACHE / f"{iso3}_{level}.json"

    if gpkg.exists() and not refresh:
        return gpd.read_file(gpkg), json.loads(side.read_text())

    url = f"https://www.geoboundaries.org/api/current/{release}/{iso3}/{level}/"
    r = requests.get(url, headers=HEADERS, timeout=90)
    if r.status_code == 404:
        raise LookupError(f"{iso3} has no {level} in {release}")
    r.raise_for_status()
    meta = r.json()
    meta = meta[0] if isinstance(meta, list) else meta

    gdf = gpd.read_file(meta["gjDownloadURL"])
    gdf.to_file(gpkg, driver="GPKG")
    side.write_text(json.dumps({
        "iso3": iso3, "level": level,
        "licence": meta["boundaryLicense"],
        "source": meta["boundarySource"],
        "source_updated": meta["sourceDataUpdateDate"],
        "features": len(gdf),
    }, indent=2))
    return gdf, json.loads(side.read_text())


gdf, meta = geoboundaries("NLD", "ADM2")
print(f"{meta['features']} features · {meta['source']} · updated {meta['source_updated']}")

344 features · National Georegister · updated Thu Jan 19 07:31:04 2023

The sidecar is what lets you answer "where did this come from" in six months, and the cache means a rerun of the pipeline does not re-download 344 polygons.

Example 2 — several countries, with the gaps reported

import pandas as pd

COUNTRIES = ["NLD", "BEL", "DEU", "LUX", "XKX"]

rows = []
for iso3 in COUNTRIES:
    for level in ["ADM1", "ADM2"]:
        try:
            gdf, meta = geoboundaries(iso3, level)
        except LookupError as exc:
            rows.append({"iso3": iso3, "level": level, "n": None, "note": str(exc)})
            continue
        rows.append({"iso3": iso3, "level": level, "n": len(gdf),
                     "note": meta["source_updated"][-4:]})

print(pd.DataFrame(rows).to_string(index=False))

iso3 level   n note
 NLD  ADM1  12 2023
 NLD  ADM2 344 2023
 BEL  ADM1   3 2023
 BEL  ADM2  11 2023
 DEU  ADM1  16 2023
 DEU  ADM2  38 2023
 LUX  ADM1  12 2023
 LUX  ADM2 105 2023
 XKX  ADM1   7 2023
 XKX  ADM2  38 2023

Read the ADM2 column before doing any cross-border work. It is 344 municipalities in the Netherlands, 11 units in Belgium and 38 in Germany — the same level number cut thirty times finer on one side of a border than the other — so "per ADM2 unit" means something different in each country. The update column reads 2023 for every country, which tells you when the source file was last touched, not that the boundaries in it are current; the explanation shows a Dutch merger it does not include.

Example 3 — boundaries from OpenStreetMap when nothing official is open

import osmnx as ox
from osmnx._errors import InsufficientResponseError

def osm_boundary(place):
    area = ox.geocode_to_gdf(place)
    if area.geom_type[0] not in ("Polygon", "MultiPolygon"):
        raise ValueError(f"{place!r} has no boundary polygon")
    return area[["geometry", "display_name", "osm_type", "osm_id", "addresstype", "class", "type"]]


for place in ["Ancoats, Manchester, England", "Greater Manchester, England"]:
    b = osm_boundary(place)
    print(f"{b['addresstype'][0]:15} {b['class'][0]}={b['type'][0]:15} {b['osm_type'][0]}/{b['osm_id'][0]:<10} "
          f"{b['display_name'][0][:40]}")

suburb          place=neighbourhood   way/532671584  Ancoats, Manchester, Greater Manchester,
state_district  boundary=administrative  relation/88084      Greater Manchester, England, United King

class and type are the honest bit. boundary=administrative means an actual administrative unit with an admin_level; anything else — Ancoats comes back as a place mapped as a single way — is an area someone drew, useful for context and unusable for anything official. Nominatim's answer also changes whenever the map is edited, so record the osm_type/osm_id you actually used. OSM boundaries are ODbL, so share-alike applies to anything you derive and redistribute as data.

Why gpd.datasets was removed

Why gpd.datasets was removed

The bundled Natural Earth extracts were a convenience that became a liability: they were several releases out of date, they encoded contested borders that GeoPandas had no business adjudicating, and people shipped analyses based on them without realising they were a teaching dataset.

Removing them was correct, and it means every tutorial written before GeoPandas 1.0 has a broken first line. For a world map, download Natural Earth from its own site, or install the geodatasets package which fetches and caches it explicitly.

Why two providers disagree about the same boundary

Three separate reasons, and they compound:

• Generalisation. A 1:10 million world layer and a national 1:10,000 layer describe the same coastline with a thousandfold difference in vertex count. Neither is wrong.
• Date. Municipalities merge. Any two sources sampled at different times differ in the places that changed.
• Politics. Disputed territories are drawn differently by different publishers, and some products offer several versions of the same country for exactly this reason.

None of these produces an error. They produce two layers that fail to line up, slivers in the overlay, and points that fall outside every polygon.

Why the join key is the hard part

The boundary geometry is rarely the problem. The identifier is. Statistics are published against official codes, and global boundary products carry their own synthetic ids instead.

Where possible, take boundaries from the same organisation that publishes your statistics — then the codes match by construction. Where that is not possible, join on the code you do share and check both directions:

Here is that check for the Dutch municipalities: CBS's 2023 statistics, from its own WFS on PDOK, against the geoBoundaries layer — which shares no code with them, so names are the only key:

import io

r = requests.get("https://service.pdok.nl/cbs/wijkenbuurten/2023/wfs/v1_0", params={
    "service": "WFS", "version": "2.0.0", "request": "GetFeature", "typeNames": "wijkenbuurten:gemeenten",
    "outputFormat": "application/json", "count": 1000}, headers=HEADERS, timeout=120)
cbs = gpd.read_file(io.BytesIO(r.content))
stats = cbs[cbs["water"] == "NEE"]                     # the layer also holds separate water polygons

missing = set(stats["gemeentenaam"]) - set(boundaries["shapeName"])
extra = set(boundaries["shapeName"]) - set(stats["gemeentenaam"])
print(f"{len(stats)} CBS municipalities, {len(boundaries)} geoBoundaries polygons")
print(f"{len(missing)} stats rows with no boundary: {sorted(missing)}")
print(f"{len(extra)} boundaries with no stats:  {sorted(extra)}")

342 CBS municipalities, 344 geoBoundaries polygons
4 stats rows with no boundary: ['Bergen (L.)', 'Bergen (NH.)', 'Hengelo', 'Voorne aan Zee']
6 boundaries with no stats:  ['Bergen (L)', 'Bergen (NH)', 'Brielle', 'Hellevoetsluis', 'Hengelo (O)', 'Westvoorne']

Two different failures in one list. Bergen (L.) against Bergen (L), and Hengelo against Hengelo (O), are the same places spelled two ways. Voorne aan Zee is not: it was formed on 1 January 2023 from Brielle, Hellevoetsluis and Westvoorne, and the geoBoundaries layer still has the three old municipalities. No amount of name cleaning fixes that — the boundary layer is simply older than the statistics. That check, run in both directions, is what GeoPandas merge returns NaN exists to prevent.

Why to store boundaries as a versioned asset

A boundary layer is an input, and inputs that change silently make results irreproducible. Treat the downloaded file as a dated artefact:

data/boundaries/NLD_ADM2_2023-01-19.gpkg
data/boundaries/NLD_ADM2_latest.gpkg -> NLD_ADM2_2023-01-19.gpkg

Then last year's figures can still be regenerated against last year's boundaries, which is the difference between an analysis and a one-off.

Remote GeoJSON, releases and missing levels

• gpd.read_file() on a remote GeoJSON URL works but downloads on every call. Cache to a local GeoPackage as in Example 1.
• geoBoundaries has several releases: gbOpen (open licences only), gbHumanitarian (UN OCHA), gbAuthoritative. They are different data, not different formats.
• A missing level is a 404, not an empty result. Catch it, or a country loop dies partway through.
• Coastal boundaries may or may not include water. Súdwest-Fryslân is 907.9 km² in geoBoundaries — CBS's total including its share of the IJsselmeer — while CBS publishes 522.7 km² of land for it. Check .area against the published figure, and against the right one of the two.
• Multipart geometries are normal. Countries with islands are MultiPolygon, and code that assumes one ring per row will silently drop parts.
• Names are not keys. Accents, hyphens, transliterations and "Saint"/"St." all vary between sources.
• Simplify for display only. A simplified boundary is fine for a map and wrong for point-in-polygon, because points near the edge change which polygon they land in.

Internal links

• GIS data sources explained — authority, licence and provenance for boundary data
• How to download data from a WFS service in Python — the route for national portals
• Spatial web services explained — what a national portal is likely to offer
• How to cache downloaded GIS data so you fetch it once — versioning a boundary layer properly
• The modifiable areal unit problem explained — why the choice of level changes the answer
• GeoPandas merge returns NaN or no matches — the join failure boundaries usually cause
• How to calculate area and distance in GeoPandas correctly — after reprojecting out of EPSG:4326
• How to aggregate spatial data by region in GeoPandas — what the boundaries are for
• How to download census boundaries in Python — statistical rather than administrative zones

FAQ

Why does gpd.datasets.get_path() no longer work?

The bundled datasets were removed in GeoPandas 1.0. Download Natural Earth from its own site, use the geodatasets package, or use a boundary API as shown here.

What is the difference between ADM1 and ADM2?

They are relative levels of subdivision, not fixed categories. ADM2 is a municipality in the Netherlands and a county in the United States. Never assume they are comparable across countries.

Which source should I use for one country?

The national mapping agency or statistics office, if it publishes openly. Its codes will match the statistics you want to join, which is usually the deciding factor.

Are geoBoundaries boundaries authoritative?

They are compiled from official sources and carry the source name and update date, but they are not the official product. For reporting against official areas, use the national source.

Why don't my boundaries line up with another layer?

Different generalisation, different vintage, or a genuine political disagreement. Check the scale and date of both before assuming one is broken.

Can I use OpenStreetMap boundaries?

Yes, with two caveats: check addresstype is administrative rather than a colloquial place, and remember ODbL share-alike applies to derived data you redistribute.

How do I keep boundaries reproducible?

Store the downloaded file under a dated name, keep a latest symlink, and record the source and update date in a sidecar. Then old results can be regenerated against old boundaries.