How to Split Multipart Geometries into Single Parts in GeoPandas #

Problem statement #

A "parcel" layer has 4,200 rows but 5,900 polygons. One row holds an island group, another holds three disconnected fragments of the same estate, and a third is a MultiLineString that is really four separate river reaches. Every per-feature calculation you run — area, centroid, nearest neighbour, count-per-zone — silently treats those bundles as one thing.

>>> gdf.geom_type.value_counts()
Polygon         2951
MultiPolygon    1249
Name: count, dtype: int64

>>> gdf.geometry.apply(lambda g: len(g.geoms) if g.geom_type.startswith("Multi") else 1).sum()
5903

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 produced GeometryCollection rows
• 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:

1. give every row a stable id before exploding
2. explode(index_parts=False, ignore_index=True) for a flat result
3. explode(index_parts=True) when you need to know which part is which
4. filter the parts you want (drop slivers, keep polygons only)
5. 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")
gdf["feature_id"] = range(len(gdf))

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()}")

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 #

Step-by-step solution #

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()))

If 40 rows out of 4,200 are multipart, you may 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 one

Use 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])                 # [(0, 0), (0, 1), (1, 0), ...]

# 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"})

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

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 0.4 m² shard left by a bad digitising session.

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 = metric.loc[~slivers].to_crs(polys.crs)

Always measure area in a projected CRS. In EPSG:4326 the numbers are square degrees and the threshold means nothing.

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

Copying a total to every part and then summing is the single most common error introduced by exploding — and it inflates results silently.

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

dissolve is the inverse of explode for geometry, and it lets you choose the aggregation for each attribute rather than guessing.

Code examples #

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}")

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)

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)

# 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(detailed.groupby(["group_name", "island_no"]).size().head())

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

print(reaches.groupby("river_id")["length_m"].agg(["count", "sum"]).head())

The per-reach lengths sum to the original multipart length, which is a quick check that nothing was lost.

Explanation #

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.

Edge cases or notes #

• ignore_index=True discards 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_empty after 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_index reflects storage order, not size or importance. Sort by area if you need "the largest part".
• Round-tripping is not byte-identical: explode then dissolve produces an equivalent multipart geometry, but node order and ring orientation may differ.
• 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.