How to Dissolve Polygons in Python (GeoPandas)
On this page
Problem statement
A dissolve merges polygons that share a value — districts into regions, parcels into zones, land-use patches into classes — and summarises their attributes on the way. In GeoPandas it is one call, dissolve(by=...), and its defaults make three decisions you may not want. Rows whose key is null disappear. Every column you did not name an aggregation for keeps the value of whichever row came first. And the merged shape of a group is whatever the union gives: often several parts, sometimes with holes, never "one polygon per region" by guarantee.
This page dissolves Germany's 400 districts into its 16 Länder, France's départements into regions with population summed, and the Isle of Wight's OpenStreetMap land use by class, and measures each of those decisions.
Quick answer
Group by the field and name an aggregation for every column you keep:
import geopandas as gpd
gdf = gpd.read_file("data/de_districts.shp")
dissolved = gdf.dissolve(by="region").reset_index()
print(len(gdf), "districts ->", len(dissolved), "regions")400 districts -> 16 regions
This creates one merged polygon or multipart geometry for each unique region value — here Germany's 400 districts (NUTS 3) into its 16 Länder. If you want to summarize other fields at the same time, pass aggregation functions with aggfunc.
Inspect the key, dissolve, aggregate on purpose
Load polygon data into GeoPandas
Start by reading a polygon layer from a shapefile, GeoJSON, or GeoPackage.
import geopandas as gpd
gdf = gpd.read_file("data/de_districts.shp")
print(gdf.head())
print(gdf.geom_type.unique())
print(gdf.crs) nuts_id ... geometry
0 DE914 ... POLYGON ((10.75931 52.79583, 11.00878 52.49675...
1 DE916 ... POLYGON ((10.20409 52.00723, 10.28617 52.03383...
2 DE917 ... POLYGON ((10.93454 52.4718, 11.04657 52.37427,...
3 DE918 ... POLYGON ((10.06553 51.92736, 10.15371 51.79354...
4 DE91A ... POLYGON ((10.29191 52.4479, 10.42407 52.33023,...
[5 rows x 4 columns]
<ArrowStringArray>
['Polygon', 'MultiPolygon']
Length: 2, dtype: str
EPSG:4326
Confirm it is polygons, and valid:
if not all(geom_type in ["Polygon", "MultiPolygon"] for geom_type in gdf.geom_type.unique()):
raise ValueError("Input layer must contain polygon geometries.")
print(f"{len(gdf)} polygons, {(~gdf.is_valid).sum()} invalid")400 polygons, 0 invalid
Inspect the attribute used for dissolving
Before running a dissolve, inspect the field you want to group by.
print(gdf.columns)
print(gdf["region"].value_counts(dropna=False))Index(['nuts_id', 'name', 'region', 'geometry'], dtype='str')
region
Bayern 96
Nordrhein-Westfalen 53
Niedersachsen 45
Baden-Württemberg 44
Rheinland-Pfalz 36
Hessen 26
Thüringen 22
Brandenburg 18
Schleswig-Holstein 15
Sachsen-Anhalt 14
Sachsen 13
Mecklenburg-Vorpommern 8
Saarland 6
Bremen 2
Berlin 1
Hamburg 1
Name: count, dtype: int64
Every distinct spelling becomes its own group, so stray spaces or case differences split a region in two. Stripping spaces is always safe; changing case is not, because str.title() rewrites names that were already right:
import pandas as pd
print(gdf["region"].str.strip().nunique(), "regions after strip")
print(pd.Series(["Nordrhein-Westfalen", "Baden-Württemberg", "Provence-Alpes-Côte d’Azur"]).str.title().tolist())16 regions after strip
['Nordrhein-Westfalen', 'Baden-Württemberg', 'Provence-Alpes-Côte D’Azur']
The apostrophe in the French region name ends a "word" for str.title(), so d’Azur becomes D’Azur — a new spelling that would split one group into two if only some rows were cleaned. Normalise case with str.casefold() into a separate key column if you must, and keep the original for labels.
If the grouping field contains nulls, decide how to handle them before dissolving. By default dissolve drops them silently, like a pandas groupby — the polygon and its area disappear from the output. Pass dropna=False to keep them as a group of their own, or drop them deliberately:
gdf = gdf.dropna(subset=["region"])Dissolve polygons by attribute
To merge polygons with the same attribute value, use dissolve(by="column_name").
dissolved = gdf.dissolve(by="region")
print(dissolved.geom_type.value_counts())Polygon 12
MultiPolygon 4
Name: count, dtype: int64
Sixteen rows, one per Land — and four of them came out as MultiPolygon; Example 1 shows why.
Control how non-geometry fields are aggregated
If your input has other fields, GeoPandas can summarize them during the dissolve.
For example, to dissolve France's 101 départements into their 27 NUTS 2 regions, sum population, and keep one representative region_name value:
fr = gpd.read_file("data/fr_departements.geojson")
regions = fr.dissolve(
by="region_id",
aggfunc={
"population": "sum",
"region_name": "first"
}
)
print(len(fr), "->", len(regions))101 -> 27
Choose by the meaning of the column: counts and totals sum, a label identical within the group can take first, and a rate or density cannot be summed or plainly averaged — recompute it from the summed numerator and denominator.
Reset the index after dissolve
By default, the dissolve field becomes the index in the output.
print(dissolved.index[:3])Index(['Baden-Württemberg', 'Bayern', 'Berlin'], dtype='str', name='region')
This is normal, but it is usually easier to inspect and export the data after resetting the index:
dissolved = dissolved.reset_index()
print(dissolved.columns.tolist())['region', 'geometry', 'nuts_id', 'name']
Every other column survives with the default aggfunc="first": nuts_id and name now hold the code and name of whichever district came first in each Land. They look like data and describe nothing. Drop columns you did not aggregate on purpose, or name an aggregation for each.
Save the dissolved output
Write it as GeoPackage unless something downstream needs a shapefile, and see what the shapefile does to the French regions:
dissolved.to_file("output/dissolved_regions.gpkg", layer="regions", driver="GPKG")
regions.reset_index().to_file("output/fr_regions.shp")
print(gpd.read_file("output/fr_regions.shp").columns.tolist())['region_id', 'population', 'region_nam', 'geometry']
That is the shapefile's ten-character field limit: region_name comes back as region_nam, with only a warning at write time.
German districts, French départements and Wight land use dissolved
The data used on this page
Three real layers, built by the block below into data/:
- German districts — the 400 NUTS 3 regions from Eurostat's GISCO NUTS 2024 boundaries (1:20 million), with the name of each district's Land as
region, saved as a shapefile. - French départements — the 101 French NUTS 3 regions with their NUTS 2 region and their population on 1 January 2024 from Eurostat table
demo_r_pjanaggr3. - Isle of Wight land use — every OpenStreetMap polygon with a
landusetag, from the Geofabrik extract.
import io
import os
import geopandas as gpd
import pandas as pd
import requests
os.makedirs("data", exist_ok=True)
os.makedirs("output", exist_ok=True)
GISCO = ("https://gisco-services.ec.europa.eu/distribution/v2/nuts/geojson/"
"NUTS_RG_20M_2024_4326_LEVL_{}.geojson")
if not os.path.exists("data/de_districts.shp"):
nuts = {level: gpd.read_file(GISCO.format(level)) for level in (1, 2, 3)}
names = {level: nuts[level].set_index("NUTS_ID")["NAME_LATN"] for level in (1, 2)}
de = nuts[3][nuts[3]["CNTR_CODE"] == "DE"]
gpd.GeoDataFrame({"nuts_id": de["NUTS_ID"], "name": de["NAME_LATN"],
"region": de["NUTS_ID"].str[:3].map(names[1])},
geometry=de.geometry.values, crs=4326).to_file("data/de_districts.shp")
csv = requests.get("https://ec.europa.eu/eurostat/api/dissemination/sdmx/2.1/data/"
"demo_r_pjanaggr3/A.NR.T.TOTAL.", timeout=180,
params={"format": "SDMX-CSV", "startPeriod": "2024", "endPeriod": "2024"},
headers={"User-Agent": "spatialworkflow.io guide"}).text
population = pd.read_csv(io.StringIO(csv)).set_index("geo")["OBS_VALUE"]
population.to_csv("data/eurostat_population_2024.csv")
fr = nuts[3][nuts[3]["CNTR_CODE"] == "FR"]
gpd.GeoDataFrame({"nuts_id": fr["NUTS_ID"], "name": fr["NAME_LATN"],
"region_id": fr["NUTS_ID"].str[:4],
"region_name": fr["NUTS_ID"].str[:4].map(names[2]),
"population": fr["NUTS_ID"].map(population)},
geometry=fr.geometry.values, crs=4326).to_file("data/fr_departements.geojson")
if not os.path.exists("data/iow_landuse.gpkg"):
pbf = "isle-of-wight-260101.osm.pbf"
if not os.path.exists(pbf):
with requests.get("https://download.geofabrik.de/europe/united-kingdom/england/"
"isle-of-wight-260101.osm.pbf", stream=True, timeout=300) as r:
with open(pbf, "wb") as f:
for chunk in r.iter_content(1 << 20):
f.write(chunk)
lu = gpd.read_file(pbf, layer="multipolygons", columns=["osm_id", "osm_way_id", "landuse"])
lu[lu["landuse"].notna()].to_file("data/iow_landuse.gpkg")
for f in sorted(os.listdir("data")):
print(f)de_districts.cpg
de_districts.dbf
de_districts.prj
de_districts.shp
de_districts.shx
eurostat_population_2024.csv
fr_departements.geojson
iow_landuse.gpkg
Example 1: Dissolve polygons by a single attribute
This example reads the German districts and dissolves them by region, then looks at what the merge produced.
import geopandas as gpd
import shapely
gdf = gpd.read_file("data/de_districts.shp")
dissolved = gdf.dissolve(by="region").reset_index()
dissolved["parts"] = shapely.get_num_geometries(dissolved.geometry.values)
dissolved["holes"] = [sum(shapely.get_num_interior_rings(shapely.get_parts(g)))
for g in dissolved.geometry]
dissolved["km2"] = dissolved.to_crs("EPSG:3035").area / 1e6
print(dissolved.loc[(dissolved.parts > 1) | (dissolved.holes > 0),
["region", "parts", "holes", "km2"]].round(0).to_string(index=False)) region parts holes km2
Brandenburg 1 1 29465.0
Bremen 2 0 404.0
Hamburg 2 0 760.0
Niedersachsen 2 1 47478.0
Schleswig-Holstein 6 0 15456.0
This creates one output feature for each unique region, and a quarter of them are not simple polygons. Bremen is two parts because the city of Bremerhaven, 50 km downriver, belongs to it; Schleswig-Holstein's parts are its islands, North Sea and Baltic, and Niedersachsen's second part is Borkum. Brandenburg has a hole where Berlin is, and Niedersachsen one where Bremen is — a dissolve keeps an enclave as a hole rather than filling it.
Example 2: Dissolve polygons and aggregate attribute values
This example dissolves the French départements by NUTS 2 region, sums their 1 January 2024 population, and checks the sums against Eurostat's own regional figures.
import geopandas as gpd
import pandas as pd
gdf = gpd.read_file("data/fr_departements.geojson")
dissolved = gdf.dissolve(
by="region_id",
aggfunc={
"population": "sum",
"region_name": "first"
}
).reset_index()
eurostat = pd.read_csv("data/eurostat_population_2024.csv", index_col="geo")["OBS_VALUE"]
dissolved["eurostat"] = dissolved["region_id"].map(eurostat)
print(dissolved.sort_values("population", ascending=False)
[["region_id", "region_name", "population"]].head(3).to_string(index=False))
print(f"{len(dissolved)} regions; {(dissolved.population != dissolved.eurostat).sum()} "
f"differ from Eurostat's regional totals")region_id region_name population
FR10 Ile-de-France 12507542
FRK2 Rhône-Alpes 6868920
FRL0 Provence-Alpes-Côte d’Azur 5252362
27 regions; 0 differ from Eurostat's regional totals
This is a common pattern when building higher-level administrative boundaries from smaller polygons, and a sum is the safe case: the parts partition the whole, so the totals must add up exactly. They do. French NUTS 2 regions are the pre-2016 regions such as Rhône-Alpes, not today's 13 administrative régions — check which geography a code list means before trusting a name.
Example 3: Export the dissolved layer to a new file
This example dissolves OpenStreetMap land use on the Isle of Wight by class, repairs the inputs first, and writes the result to GeoJSON.
import geopandas as gpd
import pandas as pd
gdf = gpd.read_file("data/iow_landuse.gpkg").to_crs("EPSG:27700")
print(f"{len(gdf):,} polygons, {(~gdf.is_valid).sum()} invalid")
as_is = gdf.dissolve(by="landuse")
repaired = gdf.assign(geometry=gdf.make_valid()).dissolve(by="landuse")
print("invalid after dissolve:", (~as_is.is_valid).sum())
structure = gdf.assign(geometry=gdf.make_valid(method="structure")).dissolve(by="landuse")
gaps = pd.DataFrame({"linework_ha": (repaired.area - as_is.area) / 1e4,
"structure_ha": (structure.area - as_is.area) / 1e4})
print(gaps[gaps.abs().max(axis=1) > 0.01].round(2).to_string())
total = repaired.area.sum() / 1e4
union = repaired.union_all().area / 1e4
print(f"classes sum to {total:,.1f} ha; their union is {union:,.1f} ha")
r = repaired.reset_index()[["landuse", "geometry"]]
pairs = gpd.overlay(r, r.rename(columns={"landuse": "other"}), keep_geom_type=True)
pairs = pairs[pairs["landuse"] < pairs["other"]]
print((pairs.set_index(["landuse", "other"]).area / 1e4).nlargest(3).round(1).to_string())
repaired.reset_index().to_file("output/land_use_dissolved.geojson", driver="GeoJSON")4,984 polygons, 10 invalid
invalid after dissolve: 0
Empty DataFrame
Columns: [linework_ha, structure_ha]
Index: []
classes sum to 29,963.7 ha; their union is 29,330.7 ha
landuse other
meadow military 201.8
grass residential 49.6
forest grass 43.7
The dissolve does not fail on the invalid inputs: the union repairs them on the way, and the result is valid. How they were repaired is still a choice, so the example compares: on this extract both make_valid() modes, linework and structure, agree with the plain union to within 0.01 ha in every class, which is why the table of differences is empty. Nothing would have said so if they had not — the difference only shows if you compare.
The last line is a different problem. OSM land use polygons overlap — the largest overlap here is meadow inside a military area — so the classes sum to more than the land they cover. Dissolving by class merges polygons within a class; it does nothing about overlaps between classes. You can open the file directly in QGIS to see them.
What a union gives back
A dissolve is a groupby whose aggregation for the geometry column is a union, and almost everything surprising about it follows from the union. A union removes the boundaries shared within a group and keeps everything else. Pieces of a group that do not touch stay separate parts, so Bremen comes back as two (Bremerhaven is 50 km downriver) and Schleswig-Holstein as six with its islands. Territory of another group surrounded by this one stays out of it, so Brandenburg gets a hole where Berlin is. Both are correct, and both break code that expects one Polygon per row.
The union also repairs as it goes: ten invalid Isle of Wight polygons went in and none came out, without a message — and repairing them first with make_valid() gave the same class areas to 0.01 ha, so here the order made no difference; it is still worth the comparison. What a union cannot do is see across groups. Land-use classes that overlap each other still overlap after each class is dissolved, so their areas sum to 29,963.7 ha over 29,330.7 ha of ground.
If you want to merge all polygons into one feature instead of grouping by a field, use a constant field first:
import shapely
gdf = gpd.read_file("data/de_districts.shp")
gdf["group_all"] = 1
all_merged = gdf.dissolve(by="group_all").reset_index(drop=True)
outline = all_merged.geometry.iloc[0]
print(outline.geom_type, shapely.get_num_geometries(outline), "parts,",
sum(shapely.get_num_interior_rings(shapely.get_parts(outline))), "holes")MultiPolygon 8 parts, 0 holes
gdf.dissolve() with no by does the same in one step.
Holes, overlaps and the columns you did not name
- Enclaves stay holes. If you need the outer outline of a group — a region including the city it surrounds — fill the interior rings afterwards, or dissolve the enclave in with it.
- Overlaps between groups survive. Dissolving by class merges within a class only; sum class areas over a non-overlapping layer, or
overlaythe classes first. - Null keys vanish by default.
dropna=Falsekeeps them as a group of their own. - Unnamed columns default to
first. Drop what you did not aggregate on purpose. - Measure in an equal-area CRS. The Länder and Wight areas above were computed in
EPSG:3035andEPSG:27700, not in the file's degrees.
Internal links
- How to Read a Shapefile in Python with GeoPandas
- How to Export GeoJSON in Python with GeoPandas
- How to Reproject Spatial Data in Python (GeoPandas)
- Exposure counts are double-counted or missed — dissolving before a join, and why it matters
- Protected area data explained — why protected area totals need dissolving first
FAQ
How do I dissolve polygons by attribute in GeoPandas?
Use dissolve(by="field_name") on a polygon GeoDataFrame.
dissolved = gdf.dissolve(by="region").reset_index()This groups polygons by the chosen field and merges their geometries.
What is the difference between dissolve and unary union in GeoPandas?
dissolve() groups features by an attribute and can aggregate attributes. A unary union merges geometries without preserving grouped attribute structure. Use dissolve when you need one output feature per category.
Why does dissolve() make my grouping field become the index?
That is the default GeoPandas behavior. After dissolving, run:
dissolved = dissolved.reset_index()This moves the grouping field back into a normal column.
How do I keep or summarize attribute columns after dissolving polygons?
Pass aggfunc with aggregation rules:
dissolved = gdf.dissolve(
by="county",
aggfunc={"population": "sum", "name": "first"}
).reset_index()Use functions that match the meaning of each field.
How do I dissolve by more than one column at once?
Pass a list of column names to by, for example gdf.dissolve(by=["state", "land_use"]). GeoPandas groups on the unique combination of those fields and merges the geometries within each combination, producing a MultiIndex you can flatten with reset_index().
Why are tiny sliver gaps or overlaps left between polygons after dissolving?
Dissolve unions the geometries as they are, so boundaries that do not coincide exactly leave slivers — gaps become holes, overlaps vanish into the union. A topologically clean source leaves none: all 400 German districts dissolve into 8 parts (the mainland and islands) with no holes. For a messy source, repair with make_valid() and then snap to a common grid with set_precision() before dissolving; avoid buffer(0), which can silently drop parts of a polygon.
Does dissolve() keep my original CRS?
Yes. Dissolve does not reproject; the output carries the same CRS as the input. If you plan to measure the merged areas afterward, reproject to a projected CRS with to_crs() first, since area on a geographic CRS like EPSG:4326 is not in meters.
How is dissolve() related to union_all() under the hood?
For each group, dissolve() effectively calls union_all() on that group's geometries to merge them, then attaches the aggregated attributes. If you only need a single merged geometry with no grouping or attributes, calling gdf.union_all() directly is simpler.