How to Create a Choropleth Map in Python with GeoPandas

On this page

Problem statement

A choropleth colours each area by a value — density by district, unemployment by region — and in GeoPandas the drawing is one argument, plot(column="density"). The map's message is decided before and around that call: whether the values joined to every area, which classification turns values into colours, and what happens to the areas with no value. Each of those can change what a reader sees without changing a single number.

This page maps Germany's 401 districts and 38 regions from Eurostat and measures each decision: a boundary vintage that leaves seven districts blank, three classification schemes that tell three different stories, and a "natural breaks" legend that changes between runs.

Quick answer

The same values shown as a continuous gradient and grouped into discrete (quantile) classes.

If your polygon layer already contains the numeric column you want to map, the shortest workflow is:

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/de_districts.shp")

ax = gdf.plot(
    column="density",
    cmap="OrRd",
    legend=True,
    figsize=(10, 8),
    edgecolor="black",
    linewidth=0.3
)

ax.set_title("Population Density by District")
ax.set_axis_off()

plt.show()
print(len(gdf), "districts, density", gdf["density"].min(), "to", gdf["density"].max(), "people/km²")

401 districts, density 35.1 to 4828.6 people/km²

That is the whole drawing. The rest of this page is about the decisions that make it say the right thing.

Load polygons, check the column, plot, classify, save

The first item is the one that most often turns a map into a population map.

Load polygon data into a GeoDataFrame

GeoPandas can read shapefiles, GeoJSON, and many other GIS formats through read_file().

import geopandas as gpd

gdf = gpd.read_file("data/de_nuts3.gpkg")

print(gdf[["NUTS_ID", "NAME_LATN", "density"]].head(3))
print(gdf.geom_type.value_counts().to_dict())
print(gdf.crs)

  NUTS_ID                    NAME_LATN  density
0   DE254   Nürnberg, Kreisfreie Stadt   2854.4
1   DE255  Schwabach, Kreisfreie Stadt    997.9
2   DE256           Ansbach, Landkreis     95.1
{'MultiPolygon': 401}
EPSG:3035

Check that the geometry type is Polygon or MultiPolygon. A choropleth map is designed for area features, not points.

If your layer contains mixed geometry types, filter to polygons before plotting:

gdf = gdf[gdf.geom_type.isin(["Polygon", "MultiPolygon"])].copy()

If you are working with a shapefile:

gdf = gpd.read_file("data/de_districts.shp")

Check the attribute column you want to map

Before plotting, inspect the available fields and confirm the target column is numeric.

print(list(gdf.columns))
print(gdf["density"].dtype)
print(gdf["density"].describe().round(1).to_dict())

['NUTS_ID', 'NAME_LATN', 'density', 'gdp_per_head', 'geometry']
float64
{'count': 394.0, 'mean': 557.0, 'std': 732.8, 'min': 35.1, '25%': 117.8, '50%': 209.6, '75%': 725.1, 'max': 4828.6}

The count is 394 of 401: seven districts have no density. They are all in Thuringia, whose districts were redrawn after the 2021 NUTS boundaries this map uses, so the 2024 figures carry codes the 2021 geometry does not have. A join on a code is only as good as the two vintages agreeing.

If the values are stored as text, convert them with pandas:

import pandas as pd

gdf["density"] = pd.to_numeric(gdf["density"], errors="coerce")

Using errors="coerce" turns invalid values into NaN, which you can handle later.

Create a basic choropleth map

The core GeoPandas choropleth workflow uses plot() with the column argument.

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="density",
    ax=ax,
    legend=True
)

plt.show()

This colors each polygon according to the value in density.

Improve the map styling

A default plot works, but it is usually better to control color, borders, title, and axes.

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="density",
    ax=ax,
    cmap="YlGnBu",
    legend=True,
    edgecolor="black",
    linewidth=0.2
)

ax.set_title("Population Density by Administrative Area", fontsize=14)
ax.set_axis_off()

plt.show()

A sequential colormap such as OrRd or YlGnBu suits a value that only goes up; a diverging one needs a meaningful midpoint.

Handle missing values

Missing values often appear as blank polygons unless you define a missing-data style.

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="density",
    ax=ax,
    cmap="YlGnBu",
    legend=True,
    edgecolor="black",
    linewidth=0.2,
    missing_kwds={
        "color": "lightgrey",
        "edgecolor": "red",
        "hatch": "///",
        "label": "Missing values"
    }
)

ax.set_title("Population Density by Administrative Area")
ax.set_axis_off()

plt.show()

You can also remove records with null values before plotting:

gdf_clean = gdf[gdf["density"].notna()].copy()
print(len(gdf), "->", len(gdf_clean), "districts;", sorted(gdf.loc[gdf["density"].isna(), "NAME_LATN"]))

401 -> 394 districts; ['Eisenach, Kreisfreie Stadt', 'Ilm-Kreis', 'Saalfeld-Rudolstadt', 'Schmalkalden-Meiningen', 'Sonneberg', 'Suhl, Kreisfreie Stadt', 'Wartburgkreis']

Dropping them leaves a hole in Thuringia with no explanation; missing_kwds keeps the hole and labels it.

Classify values for clearer interpretation

A continuous scale is not always the best choice. If values are unevenly distributed, classified ranges are easier to read.

GeoPandas supports classification schemes through the scheme argument. This requires the optional mapclassify package.

Install it if needed:

pip install mapclassify

Then use a classification method such as quantiles:

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="density",
    ax=ax,
    cmap="OrRd",
    legend=True,
    scheme="quantiles",
    k=5,
    edgecolor="black",
    linewidth=0.2
)

ax.set_title("Population Density by Quantile Class")
ax.set_axis_off()

plt.show()

How many districts each scheme puts in each class decides what the map says:

import mapclassify

values = gdf["density"].dropna()
for name in ("EqualInterval", "Quantiles", "FisherJenks"):
    c = getattr(mapclassify, name)(values, k=5)
    print(f"{name:14} counts {c.counts.tolist()}  upper bounds {[round(b) for b in c.bins]}")
top = mapclassify.EqualInterval(values, k=5).bins[-2]
print("darkest equal-interval class:", sorted(gdf.loc[gdf["density"] > top, "NAME_LATN"]))

EqualInterval  counts [315, 53, 19, 5, 2]  upper bounds [994, 1952, 2911, 3870, 4829]
Quantiles      counts [79, 79, 78, 79, 79]  upper bounds [107, 162, 301, 973, 4829]
FisherJenks    counts [256, 60, 52, 24, 2]  upper bounds [386, 998, 1905, 3073, 4829]
darkest equal-interval class: ['Berlin', 'München, Kreisfreie Stadt']

Equal intervals put 315 of 394 districts in the lightest class and two — Berlin and Munich — in the darkest: a map of Germany that is pale everywhere except two cities. Quantiles give every class about 79 districts, at the price of putting 301 and 973 people/km² in neighbouring classes. Fisher-Jenks finds the gaps in the data instead.

natural_breaks is k-means with random starting points, so two runs can give two legends (Example 3). What each scheme assumes is the subject of choropleth classification explained.

What the mapclassify loop above printed.

Save the choropleth map

To export the final map to an image:

from pathlib import Path
import matplotlib.pyplot as plt

Path("output").mkdir(parents=True, exist_ok=True)

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="density",
    ax=ax,
    cmap="OrRd",
    legend=True,
    edgecolor="black",
    linewidth=0.2
)

ax.set_title("Population Density by Administrative Area")
ax.set_axis_off()

plt.savefig("output/population_density_map.png", dpi=300, bbox_inches="tight")
plt.close()
print(f"{Path('output/population_density_map.png').stat().st_size / 1e3:.0f} KB")

597 KB

Use dpi=300 for reports or print output. bbox_inches="tight" helps remove unnecessary margins.

German districts and regions as choropleths

The data used on this page

Germany's districts and regions from Eurostat (© EuroGeographics for the administrative boundaries; statistics © European Union): the NUTS 2021 boundaries at 1:20 million from GISCO, level 3 (400 districts and the "Extra-Regio" code) and level 2 (38 regions), joined to three Eurostat tables — population density for 2024 (demo_r_d3dens), GDP per inhabitant for 2023 (nama_10r_3gdp) and the 2024 unemployment rate for ages 15–74 (lfst_r_lfu3rt). The district layer is also written as a shapefile, the format the first examples read.

import io
from pathlib import Path

import geopandas as gpd
import pandas as pd
import requests

def eurostat(dataset, key, year):
    r = requests.get(f"https://ec.europa.eu/eurostat/api/dissemination/sdmx/2.1/data/{dataset}/{key}"
                     f"?format=SDMX-CSV&startPeriod={year}&endPeriod={year}", timeout=300)
    r.raise_for_status()
    return pd.read_csv(io.StringIO(r.text))

def nuts(level):
    g = gpd.read_file("https://gisco-services.ec.europa.eu/distribution/v2/nuts/geojson/"
                      f"NUTS_RG_20M_2021_4326_LEVL_{level}.geojson")
    return g[g["CNTR_CODE"] == "DE"][["NUTS_ID", "NAME_LATN", "geometry"]]

Path("data").mkdir(exist_ok=True)
density = eurostat("demo_r_d3dens", "A.PER_KM2.", 2024).set_index("geo")["OBS_VALUE"]
gdp = eurostat("nama_10r_3gdp", "A.EUR_HAB.", 2023).set_index("geo")["OBS_VALUE"]
lfs = eurostat("lfst_r_lfu3rt", "all", 2024)
unemployment = lfs[(lfs["sex"] == "T") & (lfs["age"] == "Y15-74") & (lfs["isced11"] == "TOTAL")].set_index("geo")["OBS_VALUE"]

nuts3 = nuts(3).to_crs(3035)
nuts3["density"] = nuts3["NUTS_ID"].map(density)
nuts3["gdp_per_head"] = nuts3["NUTS_ID"].map(gdp)
nuts3.to_file("data/de_nuts3.gpkg")
nuts3.to_file("data/de_districts.shp")
nuts2 = nuts(2).to_crs(3035)
nuts2["unemployment_rate"] = nuts2["NUTS_ID"].map(unemployment)
nuts2.to_file("data/de_nuts2.gpkg")
print(len(nuts3), "districts,", nuts3["density"].notna().sum(), "with a density;",
      len(nuts2), "regions,", nuts2["unemployment_rate"].notna().sum(), "with a rate")

401 districts, 394 with a density; 38 regions, 37 with a rate

Example 1: Basic choropleth from a shapefile

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/de_districts.shp")

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="density",
    ax=ax,
    cmap="OrRd",
    legend=True
)

ax.set_title("Population Density by District")
ax.set_axis_off()

plt.show()
print(list(gdf.columns)[:4], len(gdf))

['NUTS_ID', 'NAME_LATN', 'density', 'gdp_per_he'] 401

The shapefile has quietly renamed a column: gdp_per_head is eleven characters and the format allows ten, so it comes back as gdp_per_he. Code that plots column="gdp_per_head" from the shapefile fails with a KeyError; the GeoPackage keeps the name, which is why the next examples read that.

Example 2: Choropleth from GeoJSON with custom styling

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/de_nuts3.gpkg")

fig, ax = plt.subplots(figsize=(11, 8))

gdf.plot(
    column="gdp_per_head",
    ax=ax,
    cmap="YlGnBu",
    legend=True,
    edgecolor="white",
    linewidth=0.5,
    missing_kwds={"color": "lightgrey", "label": "No data"}
)

ax.set_title("GDP per Inhabitant by District, 2023 (EUR)", fontsize=14)
ax.set_axis_off()

plt.show()
print(gdf["gdp_per_head"].isna().sum(), "districts drawn as No data;",
      gdf.loc[gdf["gdp_per_head"].idxmax(), "NAME_LATN"], "highest at", gdf["gdp_per_head"].max())

7 districts drawn as No data; Wolfsburg, Kreisfreie Stadt highest at 185100.0

Example 3: Choropleth with classified value ranges

This example requires mapclassify to be installed.

import geopandas as gpd
import matplotlib.pyplot as plt

gdf = gpd.read_file("data/de_nuts2.gpkg")

fig, ax = plt.subplots(figsize=(10, 8))

gdf.plot(
    column="unemployment_rate",
    ax=ax,
    cmap="PuRd",
    legend=True,
    scheme="fisher_jenks",
    k=5,
    edgecolor="black",
    linewidth=0.2
)

ax.set_title("Unemployment Rate by Region, 2024 (%)")
ax.set_axis_off()

plt.show()

import mapclassify
rates = gdf["unemployment_rate"].dropna()
import numpy as np
print("missing:", gdf.loc[gdf["unemployment_rate"].isna(), "NAME_LATN"].tolist())
print("fisher_jenks:           ", [round(float(b), 1) for b in mapclassify.FisherJenks(rates, k=5).bins])
for seed in (1, 2, 3):
    np.random.seed(seed)
    print(f"natural_breaks (seed {seed}):", [round(float(b), 1) for b in mapclassify.NaturalBreaks(rates, k=5).bins])

missing: ['Trier']
fisher_jenks:            [2.8, 3.3, 3.7, 4.3, 5.6]
natural_breaks (seed 1): [2.8, 3.3, 3.7, 4.3, 5.6]
natural_breaks (seed 2): [2.8, 3.3, 3.7, 4.3, 5.6]
natural_breaks (seed 3): [2.9, 3.4, 3.7, 4.3, 5.6]

natural_breaks is mapclassify's NaturalBreaks, a k-means classifier with random starting points, and the third seed moved two class boundaries. A map made twice from the same data should not change its legend, so this example uses fisher_jenks, which finds the optimal breaks deterministically.

Example 4: Export the final map image

from pathlib import Path
import geopandas as gpd
import matplotlib.pyplot as plt

Path("output").mkdir(parents=True, exist_ok=True)

gdf = gpd.read_file("data/de_nuts3.gpkg")

fig, ax = plt.subplots(figsize=(12, 9))

gdf.plot(
    column="gdp_per_head",
    ax=ax,
    cmap="viridis",
    legend=True,
    edgecolor="black",
    linewidth=0.3
)

ax.set_title("GDP per Inhabitant by District, 2023")
ax.set_axis_off()

plt.savefig("output/gdp_per_head.png", dpi=300, bbox_inches="tight")
plt.close()
print(sorted(p.name for p in Path("output").glob("*.png")))

['gdp_per_head.png', 'population_density_map.png']

Three decisions between a table and a map

Skewed data mapped with equal intervals: 315 districts in the lightest class, Berlin and Munich alone in the darkest.

A choropleth makes three mappings, and each failed differently here. Area to value is a join, and a join is only as good as the two vintages agreeing: the 2024 density figures use Thuringia's redrawn district codes, the 2021 boundaries do not, and seven districts came through with no value. Value to class is the classification, and it is an analytical choice rather than a styling one: equal intervals put 315 of 394 districts in the lightest class and only Berlin and Munich in the darkest, quantiles force about 79 into each, and Fisher-Jenks follows the gaps in the data. Class to colour is the legend, and it should not depend on chance — which is why the classified example uses Fisher-Jenks rather than the randomly started natural_breaks.

The values themselves should be rates or densities. A district is coloured by its area on the page, so a raw total — people, not people per km² — mostly shows which districts are large.

Text columns, vintages and point layers

A numeric column stored as text plots as categories; convert it with pd.to_numeric(..., errors="coerce") first.
Missing values disappear by default. missing_kwds draws them in a neutral colour with a legend entry; dropping them leaves an unexplained hole.
Check boundary and statistics vintages before joining: matching codes from different years is how Thuringia lost seven districts.
Points are not areas. For point data use proportional symbols or a density surface, not a choropleth.

Internal links

If you need more background, read What Is a Choropleth Map in GIS.

If you need to prepare output data for other tools, see How to Export GeoJSON in Python with GeoPandas.

If your layer draws in the wrong place or looks distorted, check Why GeoPandas Plot Is Not Showing the Correct Map.

FAQ

How do I create a choropleth map in GeoPandas?

Load a polygon dataset with gpd.read_file(), then call gdf.plot(column="your_field", legend=True). Add cmap, edgecolor, and figsize for better output.

Does GeoPandas work with shapefiles and GeoJSON for choropleth maps?

Yes. GeoPandas supports both shapefiles and GeoJSON through read_file(), so you can use the same choropleth workflow with either format.

What column type should I use for a choropleth map?

Use a numeric column such as integer or float. If your values are strings, convert them to numeric before plotting.

How do I handle missing values in a GeoPandas choropleth?

Use missing_kwds to style polygons with null values, or remove null rows with gdf[gdf["field"].notna()].

How do I join attribute values from a CSV before making a choropleth?

Read the CSV with pandas and merge() it onto the polygon GeoDataFrame using a shared key column, keeping the GeoDataFrame on the left so the geometry is preserved: gdf = gdf.merge(df, on="region_id", how="left"). Then plot the merged numeric column as usual.

How do I add a north arrow, scale bar, or basemap to a GeoPandas choropleth?

GeoPandas does not draw those itself; use the contextily package for a basemap and matplotlib for annotations. Reproject your layer to EPSG:3857 first with gdf.to_crs(epsg=3857), then call contextily.add_basemap(ax) after plotting.

How do I control the color classification breakpoints manually?

Pass an explicit list to classification_kwds with the user_defined scheme: gdf.plot(column="density", scheme="user_defined", classification_kwds={"bins": [10, 50, 100, 500]}). This requires mapclassify and gives you fixed, reproducible class edges instead of data-driven ones.

How do I move or resize the choropleth legend?

The continuous color bar is controlled with legend_kwds, for example legend_kwds={"shrink": 0.5, "label": "Density"}. For classified maps the legend is a standard matplotlib legend, so you can reposition it afterward with ax.get_legend().set_bbox_to_anchor((1, 0.5)).