How to Split Multipart Geometries into Single Parts in GeoPandas
On this page
Problem statement
The UK's local authority districts are 361 rows but 5,514 polygons. One row is Na h-Eileanan Siar, the Western Isles, in 1,427 parts; others hold a mainland district and its offshore islands, or a town split by an estuary. Every per-feature calculation you run — area, centroid, nearest neighbour, count-per-zone — silently treats those bundles as one thing.
import geopandas as gpd
gdf = gpd.read_file("data/raw/parcels.gpkg")
print(gdf.geom_type.value_counts())
print(gdf.geometry.apply(lambda g: len(g.geoms) if g.geom_type.startswith("Multi") else 1).sum())MultiPolygon 361
Name: count, dtype: int64
5514
Note that the type counts say nothing useful here: the ONS service returns every district as a MultiPolygon, including the 238 that have a single part. Count parts, not types.
Splitting a multipart geometry into its parts — "exploding" it — is the standard fix. The complication is what happens to the attributes: each part inherits the whole row, so sums, joins and counts can double-count unless you keep track of where each part came from.
You need this when:
- computing area, length or centroid per real-world feature rather than per record
- a spatial join must match individual islands or fragments
- a format or a consumer requires single-part geometries
make_valid()or an overlay has producedGeometryCollectionrows- you want to find and drop tiny sliver fragments hiding inside a multipart shape
Quick answer
Use GeoDataFrame.explode(), keep an id, and decide what to do with the attributes:
- give every row a stable id before exploding
explode(index_parts=False, ignore_index=True)for a flat resultexplode(index_parts=True)when you need to know which part is which- filter the parts you want (drop slivers, keep polygons only)
- aggregate back with
dissolve(by=...)when you need one row per original feature again
import geopandas as gpd
gdf = gpd.read_file("data/raw/parcels.gpkg") # feature_id is the ONS district code
parts = gdf.explode(index_parts=False, ignore_index=True)
print(f"{len(gdf)} features → {len(parts)} parts")
print(parts.geom_type.value_counts())
per_feature = parts.groupby("feature_id").size()
print(f"parts per feature: mean {per_feature.mean():.2f}, max {per_feature.max()}")361 features → 5514 parts
Polygon 5514
Name: count, dtype: int64
parts per feature: mean 15.27, max 1427
feature_id is the whole trick. With it, every part can be traced back, aggregated back, and audited; without it, an exploded frame is a one-way transformation.
What explode actually does
Count, explode, and put it back together
Count the parts before you split
Knowing the size of the change ahead of time tells you whether exploding is worth it at all.
import geopandas as gpd
def part_count(geom) -> int:
if geom is None or geom.is_empty:
return 0
return len(geom.geoms) if hasattr(geom, "geoms") else 1
gdf["n_parts"] = gdf.geometry.map(part_count)
print(gdf["n_parts"].value_counts().sort_index().head())
print("multipart rows:", (gdf["n_parts"] > 1).sum())
print("total parts :", int(gdf["n_parts"].sum()))n_parts
1 238
2 20
3 19
4 9
5 8
Name: count, dtype: int64
multipart rows: 123
total parts : 5514
Here a third of the rows are multipart, so exploding the layer is the right call. If 40 rows out of 4,000 were, you might prefer to handle those 40 explicitly rather than reshaping the whole layer.
Add a stable id first
explode() copies attributes to every part, so without an id the parts are indistinguishable.
gdf = gdf.reset_index(drop=True)
gdf["feature_id"] = gdf.index # or an existing primary key, if you have oneUse a real key when the data has one — a parcel reference, a UPRN, a GSS code. A positional index is fine for a one-off, but it changes if the input order changes.
Explode, and choose your index
# flat result: a fresh 0..n index, no memory of parts
parts = gdf.explode(index_parts=False, ignore_index=True)
# keeps a MultiIndex of (original index, part number)
parts = gdf.explode(index_parts=True)
print(parts.index[:5].tolist())
# turn that MultiIndex into ordinary columns
parts = gdf.explode(index_parts=True).reset_index()
parts = parts.rename(columns={"level_0": "src_index", "level_1": "part_index"})
print(parts[["src_index", "part_index", "name"]].head(3).to_string(index=False))[(0, 0), (0, 1), (0, 2), (0, 3), (0, 4)]
src_index part_index name
0 0 Hartlepool
0 1 Hartlepool
0 2 Hartlepool
index_parts=True is worth the extra column whenever a reviewer might ask "which part of which feature is this?" — which is most of the time in cleaning work.
Handle GeometryCollections
explode() splits collections too, which is exactly how you flatten the output of make_valid() or an overlay. The parts can then be of mixed types, so filter after exploding.
POLYGONAL = ["Polygon", "MultiPolygon"]
parts = gdf.explode(index_parts=False, ignore_index=True)
print(parts.geom_type.value_counts())
polys = parts[parts.geom_type.isin(POLYGONAL)].copy()
dropped = parts[~parts.geom_type.isin(POLYGONAL)]
if len(dropped):
print(f"dropping {len(dropped)} non-polygonal parts")
dropped.to_file("data/out/dropped_parts.gpkg", driver="GPKG")Polygon 5514
Name: count, dtype: int64
These boundaries are clean, so there is nothing to drop; the output of make_valid() or overlay() is where collections with stray lines and points turn up.
Writing the discards out makes the decision reviewable rather than invisible.
Drop slivers, not real parts
Exploding often reveals fragments that were never meaningful: a shard left by a bad digitising session, or a skerry nobody meant to count.
metric = polys.to_crs(polys.estimate_utm_crs())
metric["area_m2"] = metric.geometry.area
MIN_AREA = 5.0 # square metres — your call, per dataset
slivers = metric["area_m2"] < MIN_AREA
print(f"{slivers.sum()} parts below {MIN_AREA} m² "
f"({metric.loc[slivers, 'area_m2'].sum():.1f} m² total)")
keep = polys.loc[~slivers].copy() # the original geometry, not a round trip
print(f"smallest kept part {metric.loc[~slivers, 'area_m2'].min():.1f} m²; "
f"parts under 1,000 m²: {(metric['area_m2'] < 1000).sum()}")51 parts below 5.0 m² (95.2 m² total)
smallest kept part 5.1 m²; parts under 1,000 m²: 1005
In these boundaries the smallest parts are rocks off the coast, kept by the generalisation, and one part has no area at all. Whether 5 m² is the right line depends on what the layer is for.
Always measure area in a projected CRS. In EPSG:4326 the numbers are square degrees and the threshold means nothing. Measure on a projected copy and filter the original with the result: converting the parts to UTM and back would move every vertex by a fraction of a millimetre, and nothing downstream would compare equal to the source again.
Recompute per-part attributes
Attributes copied from the parent are now wrong for anything additive.
metric = keep.to_crs(keep.estimate_utm_crs())
keep["part_area_m2"] = metric.geometry.area.values
feature_area = keep.groupby("feature_id")["part_area_m2"].transform("sum")
keep["part_share"] = keep["part_area_m2"] / feature_area
# an attribute like population must be apportioned, not copied
keep["population_part"] = keep["population"] * keep["part_share"]
print(f"population by district: {gdf['population'].sum():>15,.0f}")
print(f"summed over parts, copied: {keep['population'].sum():>15,.0f}")
print(f"summed over parts, apportioned: {keep['population_part'].sum():>15,.0f}")population by district: 69,483,939
summed over parts, copied: 606,787,372
summed over parts, apportioned: 69,483,939
Copying a total to every part and then summing is the single most common error introduced by exploding — and it inflates results silently: here the UK's population comes out almost nine times too large, 606.8 million instead of 69.5 million. Apportioning by area gets the total right; it still spreads each district's people evenly over its islands, which is only as true as that assumption.
Put it back together when you are done
# one row per original feature again, parts merged into a multipart geometry
rejoined = keep.dissolve(by="feature_id", aggfunc={"part_area_m2": "sum"}).reset_index()
print(f"{len(keep)} parts → {len(rejoined)} features")
original = gdf.set_index("feature_id").geometry
back = rejoined.set_index("feature_id").geometry
same = back.geom_equals(original.loc[back.index])
exact = back.geom_equals_exact(original.loc[back.index], tolerance=0)
print(f"same shape as before: {same.sum()} of {len(back)}; identical coordinates in order: {exact.sum()}")
dropped_from = set(polys.loc[slivers, "feature_id"])
print(f"districts that lost a sliver: {len(dropped_from)}")
round_trip = gdf.geometry.to_crs(gdf.estimate_utm_crs()).to_crs(gdf.crs)
print(f"after a UTM round trip, identical: {gdf.geometry.geom_equals_exact(round_trip, tolerance=0).sum()}"
f" of {len(gdf)}; within 1 µ°: {gdf.geometry.geom_equals_exact(round_trip, tolerance=1e-6).sum()}")5463 parts → 361 features
same shape as before: 334 of 361; identical coordinates in order: 0
districts that lost a sliver: 27
after a UTM round trip, identical: 0 of 361; within 1 µ°: 361
dissolve is the inverse of explode for geometry, and it lets you choose the aggregation for each attribute rather than guessing.
Exploding the UK's local authority districts
The data used on this page
parcels.gpkg: the UK's 361 local authority districts, December 2024, from the ONS Open Geography Portal (generalised to 20 m, clipped to the coastline), with the ONS code asfeature_idand the mid-2025 population estimate from Nomis.islands.gpkg: three of those districts — Orkney, Shetland and Na h-Eileanan Siar — as island groups.stations.gpkg: the airfields in those three groups from OurAirports.rivers.gpkg: the multipart river centrelines from Natural Earth's 1:10m rivers.
import io
import os
import geopandas as gpd
import pandas as pd
import requests
os.makedirs("data/raw", exist_ok=True)
os.makedirs("data/out", exist_ok=True)
ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
"Local_Authority_Districts_December_2024_Boundaries_UK_BGC/FeatureServer/0/query")
if not os.path.exists("data/raw/parcels.gpkg"):
features, offset = [], 0
while True:
page = requests.get(ONS, timeout=300, params={
"where": "1=1", "outFields": "LAD24CD,LAD24NM", "outSR": 4326, "f": "geojson",
"resultOffset": offset, "resultRecordCount": 50}).json()["features"]
features += page
offset += 50
if len(page) < 50:
break
lad = gpd.GeoDataFrame.from_features(features, crs=4326)
csv = requests.get("https://www.nomisweb.co.uk/api/v01/dataset/NM_2002_1.data.csv", timeout=120,
params={"geography": "TYPE424", "date": "2025", "gender": "0", "c_age": "200",
"measures": "20100", "select": "geography_code,obs_value"}).text
population = pd.read_csv(io.StringIO(csv)).set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
parcels = gpd.GeoDataFrame({"feature_id": lad["LAD24CD"], "name": lad["LAD24NM"],
"population": lad["LAD24CD"].map(population)},
geometry=lad.geometry, crs=4326)
parcels.to_file("data/raw/parcels.gpkg")
groups = {"S12000023": "Orkney", "S12000027": "Shetland", "S12000013": "Western Isles"}
islands = parcels[parcels["feature_id"].isin(groups)].copy()
islands["group_name"] = islands["feature_id"].map(groups)
islands.to_file("data/raw/islands.gpkg")
airports = pd.read_csv("https://davidmegginson.github.io/ourairports-data/airports.csv")
airports = airports[(airports["iso_country"] == "GB") & airports["type"].isin(
["small_airport", "medium_airport", "large_airport"])]
airports = gpd.GeoDataFrame(airports[["ident", "name"]], crs=4326, geometry=gpd.points_from_xy(
airports["longitude_deg"], airports["latitude_deg"]))
airports[airports.intersects(islands.buffer(0.02).union_all())].to_file("data/raw/stations.gpkg")
rivers = gpd.read_file("https://naciscdn.org/naturalearth/10m/physical/"
"ne_10m_rivers_lake_centerlines.zip")
rivers[rivers.geom_type == "MultiLineString"][["name", "geometry"]].to_file("data/raw/rivers.gpkg")
for name in ("parcels", "islands", "stations", "rivers"):
layer = gpd.read_file(f"data/raw/{name}.gpkg")
print(f"{name:9} {len(layer):>4} rows {layer.geom_type.value_counts().to_dict()}")parcels 361 rows {'MultiPolygon': 361}
islands 3 rows {'MultiPolygon': 3}
stations 22 rows {'Point': 22}
rivers 420 rows {'MultiLineString': 420}
Example 1: a reusable explode-with-audit helper
import geopandas as gpd
def explode_parts(gdf: gpd.GeoDataFrame, id_col: str = "feature_id",
keep_types=("Polygon", "MultiPolygon"), min_area_m2: float | None = None):
"""Explode to single parts, keeping a trace back to the parent feature."""
src = gdf.reset_index(drop=True).copy()
if id_col not in src.columns:
src[id_col] = src.index
parts = src.explode(index_parts=True).reset_index(level=1)
parts = parts.rename(columns={"level_1": "part_index"}).reset_index(drop=True)
parts = parts[parts.geometry.notna() & ~parts.geometry.is_empty]
report = {"features": len(src), "parts_raw": len(parts)}
if keep_types:
wrong = ~parts.geom_type.isin(keep_types)
report["dropped_wrong_type"] = int(wrong.sum())
parts = parts.loc[~wrong].copy()
if min_area_m2:
metric = parts.to_crs(parts.estimate_utm_crs())
parts["part_area_m2"] = metric.geometry.area.values
small = parts["part_area_m2"] < min_area_m2
report["dropped_slivers"] = int(small.sum())
parts = parts.loc[~small].copy()
report["parts_kept"] = len(parts)
return parts, report
parts, report = explode_parts(gpd.read_file("data/raw/parcels.gpkg"), min_area_m2=5)
for k, v in report.items():
print(f"{k:20} {v}")features 361
parts_raw 5514
dropped_wrong_type 0
dropped_slivers 51
parts_kept 5463
Example 2: find the features worth exploding
import geopandas as gpd
gdf = gpd.read_file("data/raw/parcels.gpkg")
gdf["n_parts"] = gdf.geometry.map(lambda g: len(g.geoms) if hasattr(g, "geoms") else 1)
worst = gdf.nlargest(10, "n_parts")[["feature_id", "name", "n_parts"]]
print(worst.to_string(index=False))
# explode only the multipart rows, leave the rest untouched
single = gdf[gdf["n_parts"] == 1]
multi = gdf[gdf["n_parts"] > 1].explode(index_parts=False)
result = gpd.GeoDataFrame(gpd.pd.concat([single, multi], ignore_index=True), crs=gdf.crs)
print(f"{len(single)} single-part rows kept, {len(multi)} parts from the rest → {len(result)} rows")feature_id name n_parts
S12000013 Na h-Eileanan Siar 1427
S12000035 Argyll and Bute 726
S12000017 Highland 724
S12000027 Shetland Islands 552
S12000023 Orkney Islands 177
E07000074 Maldon 140
E07000146 King's Lynn and West Norfolk 98
E07000148 Norwich 96
E07000147 North Norfolk 85
E06000035 Medway 64
238 single-part rows kept, 5276 parts from the rest → 5514 rows
Exploding only what needs it keeps the row count — and the diff — as small as possible.
Example 3: exploding for an accurate spatial join
import geopandas as gpd
islands = gpd.read_file("data/raw/islands.gpkg") # MultiPolygon per island group
stations = gpd.read_file("data/raw/stations.gpkg").to_crs(islands.crs) # airfields
# without exploding: a station matches the whole group
grouped = gpd.sjoin(stations, islands[["group_name", "geometry"]], predicate="within")
print("distinct groups matched:", grouped["group_name"].nunique())
# with exploding: it matches the individual island
islands["group_id"] = range(len(islands))
single_islands = islands.explode(index_parts=True).reset_index(level=1)
single_islands = single_islands.rename(columns={"level_1": "island_no"})
detailed = gpd.sjoin(stations, single_islands[["group_name", "island_no", "geometry"]],
predicate="within")
print("distinct islands matched:", detailed.groupby(["group_name", "island_no"]).ngroups)
print(detailed.groupby("group_name")["island_no"].nunique().to_string())
print(f"{len(stations) - len(detailed)} of {len(stations)} airfields fall outside every island")
missed = stations[~stations.index.isin(detailed.index)].to_crs(27700)
shore = single_islands.to_crs(27700)
for _, airfield in missed.iterrows():
print(f" {airfield['name']}: {shore.distance(airfield.geometry).min():.0f} m from the nearest island")distinct groups matched: 3
distinct islands matched: 18
group_name
Orkney 8
Shetland 8
Western Isles 2
1 of 22 airfields fall outside every island
Barra Airport: 303 m from the nearest island
Exploding turns "3 groups" into 18 separate islands with an airfield. The one airfield that matches nothing is Barra, whose runways are the tidal beach of Tràigh Mhòr — below high water, and so outside a boundary clipped to the coastline, 303 m from it. within is strict; a small buffer is kinder to points on a coast.
Example 4: split lines into their individual reaches
import geopandas as gpd
rivers = gpd.read_file("data/raw/rivers.gpkg")
rivers["river_id"] = range(len(rivers))
reaches = rivers.explode(index_parts=True).reset_index(level=1)
reaches = reaches.rename(columns={"level_1": "reach_no"})
metric = reaches.to_crs(reaches.estimate_utm_crs())
reaches["length_m"] = metric.geometry.length.values
per_river = reaches.groupby("river_id")["length_m"].agg(["count", "sum"])
whole = rivers.to_crs(rivers.estimate_utm_crs()).length
print(per_river.assign(name=rivers["name"]).sort_values("count", ascending=False).head(4).round(0))
print("largest difference from the unexploded length:",
f"{(per_river['sum'] - whole.values).abs().max():.6f} m") count sum name
river_id
346 14 967581.0 Churchill
347 13 600092.0 Churchill
262 11 2303910.0 Volga
387 11 596654.0 NaN
largest difference from the unexploded length: 0.000000 m
The per-reach lengths sum to the original multipart length, which is a quick check that nothing was lost. estimate_utm_crs() picks one zone for the whole layer, so the lengths are only as good as that zone is for each river — fine for a check, not for a gazetteer.
Single, multipart and collection geometries
Simple Features distinguishes single geometries (Polygon, LineString, Point) from their multipart counterparts (MultiPolygon, MultiLineString, MultiPoint), plus GeometryCollection, which can hold a mixture. A multipart geometry is one geometric object made of disconnected components — the right model for an archipelago, a river with braided channels, or a landholding split by a road.
Where it goes wrong is that a GeoDataFrame row is both a record and a shape, and those two things have different natural granularities. len(gdf) counts records; the number of polygons on the map may be much larger. Any analysis that assumes "one row equals one thing on the ground" is wrong for multipart data — nearest-neighbour searches use the aggregate geometry, and a count of features under-reports the physical objects.
explode() resolves the mismatch by making each part its own record. It copies the parent's attributes into every part, which is right for descriptive fields (name, class, owner) and wrong for extensive ones (population, area, count). That is why apportioning matters: three islands sharing a population of 900 do not each have 900 residents, and summing the exploded column would report 2,700.
The inverse operation is dissolve(by=id), which recombines parts into multipart geometries and lets you specify how each attribute aggregates. Keeping a parent id from the start is what makes that round trip possible — and it also makes the intermediate state auditable, which matters when the exploded layer is what you hand to someone else.
Lost parent ids and empty parts
ignore_index=Truediscards the parent index: Convenient for a flat result, but add your own id column first or the link to the parent is gone.- Empty geometries survive the explode: They come through as empty parts. Filter with
~gdf.geometry.is_emptyafter exploding. explode()on a Series vs a GeoDataFrame:GeoSeries.explode()returns geometries only; the DataFrame method carries attributes.- Pandas has its own
explode:DataFrame.explode()expands list-like columns. GeoPandas overrides it for geometry — be explicit about which you mean in mixed code. - Part order is not guaranteed to be meaningful:
part_indexreflects storage order, not size or importance. Sort by area if you need "the largest part". - Round-tripping is not byte-identical:
explodethendissolvegave back the same shape for every district that lost no sliver, and identical coordinate order for none of them. Compare withgeom_equals, notgeom_equals_exact. - File size grows: Every part repeats the parent attributes, so an exploded shapefile can be several times larger. GeoPackage and Parquet compress this much better.
Internal links
- How to Fix Invalid Geometries in Python (GeoPandas)
- How to Dissolve Polygons in Python (GeoPandas)
- How to Find and Remove Duplicate Geometries in GeoPandas
- The Python GIS Data Cleaning Checklist: From Raw Download to Analysis-Ready
- Mixed Geometry Types Error When Saving a GeoDataFrame: How to Fix It
- How to Calculate Area and Distance in GeoPandas (Correctly)
FAQ
What is the difference between explode() and dissolve()?
explode() splits multipart geometries into one row per part; dissolve() merges rows back into multipart geometries grouped by a column. They are inverses, and keeping a parent id lets you move between the two.
Does exploding change my attribute totals?
It changes them if you sum extensive columns afterwards, because every part inherits the parent's value. Apportion such columns by area or length share before aggregating, or aggregate back by the parent id.
What does index_parts do?
index_parts=True keeps a MultiIndex of (original index, part number) so each part is traceable. index_parts=False drops it, which is fine only when you have added your own id column.
How do I explode only the multipart rows?
Filter on the part count first — gdf.geometry.map(lambda g: len(g.geoms) if hasattr(g, "geoms") else 1) > 1 — explode that subset, and concatenate it back with the untouched single-part rows.
Will exploding fix an invalid geometry?
No. Exploding separates components; it does not repair self-intersections or bad rings. Run make_valid() first, then explode the result to flatten any collection it returns.
Should I explode before writing a shapefile?
Only if the consumer requires single-part features. Shapefiles store multipart geometries perfectly well, and exploding multiplies the attribute rows, which makes the file substantially bigger.
How do I find the largest part of each feature?
Compute area in a projected CRS after exploding, then use parts.sort_values("part_area_m2").groupby("feature_id").tail(1) to keep the biggest part per original feature.