⚙️ How-To Data sources & downloads Getting started Finding & Fetching Data · step 8 of 14 Updated

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

world = gpd.read_file(gpd.datasets.get_path("naturalearth_lowres"))
AttributeError: The geopandas.dataset has been deprecated and was removed in GeoPandas 1.0.

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

r = requests.get("https://www.geoboundaries.org/api/current/gbOpen/NLD/ADM2/",
                 headers={"User-Agent": "my-project/1.0"}, 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:

Source Coverage Good for
National portal (WFS, download) one country official, authoritative, current
geoBoundaries API global, ADM0–ADM2 consistent schema, clear licence, scriptable
Natural Earth global, coarse world maps at small scale
OpenStreetMap global, variable places without an open official source
Four sources of administrative boundaries ordered by authority, from a national portal down to OpenStreetMap, with coverage and licence noted.
Start at the top for one country. Drop down a rung only when the rung above cannot cover your extent.

Step-by-step solution

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
       Aalten      NL-GE  NLD-ADM2-3_0_0-B1
       Alkmaar     NL-NH  NLD-ADM2-3_0_0-B2
       Almelo      NL-OV  NLD-ADM2-3_0_0-B3
duplicate names: 0

shapeID is a geoBoundaries identifier, stable within a release and not the code your statistics use. Joining on shapeName works here because names happen to be unique, and it will fail the moment two municipalities share a name or one is spelled with a different accent. 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  841.6
    Noordoostpolder  595.4
    Hollands Kroon   571.2

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.

ADM2 meaning a municipality in one country and a county in another, with very different unit counts.
The same level number, two different kinds of thing. Cross-country comparison at a fixed ADM level compares the labels, not the places.

Code examples

Example 1 — a boundary fetcher with metadata and caching 

import json
from pathlib import Path

import geopandas as gpd
import requests

HEADERS = {"User-Agent": "spatialworkflow-example/1.0 ([email protected])"}
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.0                              2023
 NLD  ADM2  344.0                              2023
 BEL  ADM1   11.0                              2020
 BEL  ADM2  581.0                              2020
 DEU  ADM1   16.0                              2021
 DEU  ADM2  401.0                              2021
 LUX  ADM1   12.0                              2019
 LUX  ADM2    NaN  XKX has no ADM2 in gbOpen
 XKX  ADM1    7.0                              2018
 XKX  ADM2    NaN  XKX has no ADM2 in gbOpen

Two things worth reading off this table before doing any cross-border work. The ADM2 counts vary from 344 to 581 for countries of broadly comparable size, so "per ADM2 unit" means something different in each. And the update years span 2018 to 2023 — a five-year spread inside one supposedly consistent global product.

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"]]


for place in ["Ancoats, Manchester, England", "Greater Manchester, England"]:
    b = osm_boundary(place)
    print(f"{b['addresstype'][0]:16} {b['osm_type'][0]}/{b['osm_id'][0]:<12} "
          f"{b['display_name'][0][:52]}")
neighbourhood    relation/8398124    Ancoats, Manchester, Greater Manchester, England, Un
administrative   relation/172484     Greater Manchester, England, United Kingdom

addresstype is the honest bit. administrative means an actual administrative unit with an admin_level; neighbourhood is a colloquial area someone drew, useful for context and unusable for anything official. OSM boundaries are ODbL, so share-alike applies to anything you derive and redistribute as data.

Explanation

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:

missing = set(stats["code"]) - set(boundaries["code"])
extra = set(boundaries["code"]) - set(stats["code"])
print(f"{len(missing)} stats rows with no boundary, {len(extra)} boundaries with no stats")

That two-line check is what GeoPandas merge returns NaN exists to prevent.

Three causes of boundary disagreement — generalisation, different dates, and disputed territory — each producing an overlay with slivers.
All three look identical downstream: an overlay full of slivers and points that match nothing.

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.

Edge cases or notes

  • 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. Two sources of the same municipality can differ by tens of square kilometres depending on where they stop. Check .area against the published figure.
  • 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.

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.