Dissolving and Aggregating Features by Attribute

dissolve() collapses many features into one per attribute value — parcels into districts, census blocks into tracts, road segments into named routes — and the reason it surprises people is that it performs two reductions at once, one on geometry and one on every other column, with completely different rules. This guide is for anyone rolling a fine-grained layer up to a coarser boundary in GeoPandas and needing the attribute totals to survive the trip. It sits under Spatial Joins & Merging in Geospatial Data Ingestion & Processing Workflows.

Why This Approach / What Goes Wrong

A dissolve is a groupby with a geometry reduction bolted onto it. GeoPandas groups the rows by the by column, unions each group's geometries into a single geometry, aggregates the remaining columns with aggfunc, and joins the two results back together. Nothing about that pipeline is spatial except the union — the attribute half is plain pandas, and it obeys plain pandas rules.

That split is where the damage happens, because aggfunc defaults to "first". Call parcels.dissolve(by="district") and you get geometrically correct district polygons whose population column holds the population of whichever parcel happened to sort first in each group. No warning, no error, a perfectly plausible-looking number. The output has the right shape and the wrong data, which is far more dangerous than a crash. Every production dissolve should pass an explicit aggfunc, and for anything beyond one column that means a dict mapping each column to its own reduction.

The geometry half has its own failure mode. The union is a true unary union, not a bounding-box merge: GEOS nodes every input edge against every other edge in the group, then reassembles the result. Feed it a self-intersecting ring and it does not silently skip the row — it raises GEOSException: TopologyException partway through and takes the whole dissolve with it, which is why validity repair from Topology Validation & Repair belongs before the call, not after the traceback. That noding work is also why dissolve is the slowest operation in most merge pipelines: cost scales with the total vertex count inside a group, not the row count.

Prerequisites

• geopandas>=1.0 — dissolve() here delegates the union to GeoSeries.union_all(), which replaced the deprecated unary_union property
• shapely>=2.0 — supplies make_valid, set_precision and coverage_union_all
• pandas>=2.0 — the aggregation half of the operation, including dict and named reductions
• pyogrio>=0.7 — vectorized file I/O, so reading the source layer is not the bottleneck

conda install -c conda-forge "geopandas=1.0.*" "shapely=2.0.*" "pandas=2.2.*" "pyogrio=0.7.*"

Step-by-Step Implementation

1. Build a layer with a grouping key, in a projected CRS. Every area check later depends on metric units, so the parcels are placed in UTM zone 33N rather than in degrees — and never in Web Mercator, whose area distortion makes any summed hectare figure fiction. Real pipelines read the layer with gpd.read_file(...); the synthetic boxes keep this block runnable on its own.

import geopandas as gpd
from shapely.geometry import box

parcels = gpd.GeoDataFrame(
    {
        "parcel_id": [1, 2, 3, 4, 5, 6],
        "district": ["Riverside", "Riverside", "Riverside",
                     "Hillcrest", "Hillcrest", "Hillcrest"],
        "zoning": ["R1", "R1", "C2", "R2", "R2", "R2"],
        "population": [120, 90, 60, 200, 140, 80],
    },
    geometry=[
        box(0, 0, 100, 100), box(100, 0, 200, 100), box(0, 100, 100, 200),
        box(300, 0, 400, 100), box(400, 0, 500, 100), box(600, 0, 700, 100),
    ],
    crs="EPSG:32633",          # UTM 33N — metres, so .area is m²
)

2. Run the naive dissolve and read what it actually did. This is the call most people write first, and it is worth seeing its output before improving on it.

districts = parcels.dissolve(by="district")

print(districts.geometry.geom_type)
# district
# Hillcrest    MultiPolygon
# Riverside         Polygon
# Name: geometry, dtype: object

print(districts[["parcel_id", "zoning", "population"]])
#            parcel_id zoning  population
# district
# Hillcrest          4     R2         200
# Riverside          1     R1         120

Three things to notice. The by column became the index (as_index=True is the default). Riverside's three touching boxes merged into one Polygon while Hillcrest's detached parcel forced a MultiPolygon — a dissolve routinely returns a mixed geometry column, so anything downstream must tolerate both. And population reads 200 and 120, not 420 and 270, because aggfunc="first" took one arbitrary row per group. Groups are sorted alphabetically by default (sort=True), which is what makes "first" look deterministic while still being meaningless.

3. Choose aggfunc per column, and keep a count. A dict gives each column its own reduction and — importantly — drops every column you do not name, which is usually a feature: it stops stale per-parcel columns from riding along with nonsense values. Counting group members is the same mechanism: aggregate any always-populated column with "count" and rename it.

rollup = parcels.dissolve(
    by="district",
    aggfunc={
        "population": "sum",     # additive quantity — a real total
        "parcel_id": "count",    # non-null count = features per group
        "zoning": "first",       # representative label, chosen deliberately
    },
    as_index=False,              # keep `district` as a column, not the index
).rename(columns={"parcel_id": "parcel_count"})

print(rollup[["district", "parcel_count", "population"]])
#     district  parcel_count  population
# 0  Hillcrest             3         420
# 1  Riverside             3         270

Which reduction is safe depends on the column's meaning, not its dtype. "sum" is correct for counts and populations and wrong for densities and rates. "mean" on a parcel-level rate produces an unweighted average that quietly over-weights tiny parcels — compute the numerator and denominator separately, sum both, and divide afterwards. And a reduction that looks illegal on text often is not: pandas concatenates object columns under "sum" rather than raising, so a stray "zoning": "sum" yields "R2R2R2" instead of an error.

4. Repair geometry before the union, not after the traceback. union_all() nodes every edge in the group against every other edge, and a single self-intersecting ring aborts the entire call. Repair is vectorized, so the guard costs one pass over the column.

from shapely import set_precision

# make_valid() fixes self-intersections and bowties; it can return a
# GeometryCollection, so keep only the polygonal parts before unioning.
parcels["geometry"] = parcels.geometry.make_valid()
parcels = parcels[parcels.geometry.geom_type.isin(["Polygon", "MultiPolygon"])]

# Snap coordinates onto a 1 mm grid so near-coincident shared borders
# collapse exactly, instead of leaving hairline slivers in the union.
parcels["geometry"] = parcels.geometry.apply(lambda g: set_precision(g, 0.001))

clean_districts = parcels.dissolve(by="district", aggfunc={"population": "sum"})

5. Roll joined values up to a boundary — with the right geometry. The natural sequel to a left spatial join is "now total those readings per district", and this is where dissolve is most often misused. After a join, the active geometry column is still the left layer's — so dissolving the joined frame unions the sensor points into a MultiPoint blob, not the district polygon you were picturing.

from shapely.geometry import Point

sensors = gpd.GeoDataFrame(
    {"sensor_id": [1, 2, 3, 4], "reading": [18.4, 22.1, 9.7, 15.2]},
    geometry=[Point(50, 50), Point(150, 50), Point(350, 50), Point(650, 50)],
    crs="EPSG:32633",
)

boundaries = parcels.dissolve(
    by="district", aggfunc={"population": "sum"}, as_index=False
)
joined = gpd.sjoin(
    sensors, boundaries[["district", "geometry"]], how="inner", predicate="within"
)

# Misuse: dissolve unions the sensor POINTS, discarding the boundary geometry
print(joined.dissolve(by="district").geom_type.tolist())
# ['MultiPoint', 'MultiPoint']

# Correct: reduce the values with pandas, then attach them to the boundaries
per_district = joined.groupby("district", as_index=False)["reading"].mean()
boundaries = boundaries.merge(per_district, on="district", how="left")

print(boundaries[["district", "population", "reading"]])
#     district  population  reading
# 0  Hillcrest         420    12.45
# 1  Riverside         270    20.25

Use dissolve() when you want the merged footprint of the joined features; use a plain groupby plus merge when you want their values carried onto geometry that already exists.

6. Make it fast before the layer gets big. Dissolve slows down for two separate reasons, and the fix differs. Too many vertices inside a group makes each union expensive, because GEOS must node every edge against every other edge; too many groups makes the per-group loop dominate, since the union runs once per key. Reduce vertices with set_precision or .simplify(tolerance) before dissolving, and drop every column you are not aggregating so the pandas half moves less data. When the polygons form a clean coverage — no overlaps, no gaps, shared borders identical vertex for vertex, as administrative boundaries usually are — shapely.coverage_union_all skips the general noding step entirely and is dramatically faster than a unary union on the same input.

import shapely

def dissolve_coverage(gdf: gpd.GeoDataFrame, by: str) -> gpd.GeoDataFrame:
    """Union each group as a coverage. Requires non-overlapping, gap-free input."""
    keys, geoms = [], []
    for key, group in gdf.groupby(by, sort=True):
        geoms.append(shapely.coverage_union_all(group.geometry.to_numpy()))
        keys.append(key)
    return gpd.GeoDataFrame({by: keys}, geometry=geoms, crs=gdf.crs)

The trade is strictness: coverage_union_all assumes the coverage is valid and returns a wrong answer rather than an error if it is not, so validate the assumption once on a sample before trusting it in a pipeline. Past roughly ten million rows, stop dissolving in a single process altogether. Because union is associative, the operation can run per partition and then once more over the partial results — the strategy Dask-GeoPandas applies automatically. The other escape hatch is to never materialize the frame at all: DuckDB spatial performs the identical grouping over GeoParquet on disk, streaming rather than holding the layer in memory, and returns only the dissolved result.

SELECT district,
       ST_Union_Agg(geom) AS geom,
       sum(population)     AS population,
       count(*)            AS parcel_count
FROM read_parquet('parcels.parquet')
GROUP BY district;

Verification

A dissolve is worth asserting because both of its failure modes are silent — a wrong attribute total and a quietly dropped group look exactly like a correct result. Row count, attribute totals and total area all have exact expected relationships to the input, so three cheap assertions catch nearly everything.

import numpy as np

# One output row per distinct key — anything less means keys were dropped
assert len(rollup) == parcels["district"].nunique(dropna=True)

# Additive attributes must survive the reduction unchanged
assert rollup["population"].sum() == parcels["population"].sum()

# Area is only comparable in a projected CRS
assert rollup.crs.is_projected, "Reproject before comparing areas"

# Union never invents area; it equals the input only when parts do not overlap
assert rollup.geometry.area.sum() <= parcels.geometry.area.sum() + 1e-6

# The union output must itself be valid, or the next overlay will raise
assert rollup.geometry.is_valid.all()

print(f"{len(parcels)} parcels -> {len(rollup)} districts, "
      f"{rollup.geometry.area.sum():,.0f} m2 retained")
# 6 parcels -> 2 districts, 60,000 m2 retained

If the area assertion fires with the dissolved total lower than the input, the parcels overlap — legitimate for zoning layers, a data error for cadastral ones. If population no longer sums correctly, a group key was NaN and dropna=True discarded it.

Edge Cases & Debugging

• GEOSException: TopologyException: found non-noded intersection. One ring in the group self-intersects. Run gdf.geometry = gdf.geometry.make_valid() before dissolving, and snap with set_precision(g, 0.001) if the message names coordinates on a shared border.
• Rows disappeared and totals dropped. dropna=True is the default, so every feature whose by value is NaN is discarded without a warning. Pass dropna=False to get a null group, or fillna("unassigned") first.
• Columns came back as a MultiIndex. Passing a list — aggfunc=["sum", "mean"] — gives every column both reductions under a two-level header. Use a dict for per-column control, or flatten with df.columns = ["_".join(c) for c in df.columns].
• Empty geometries appear for groups with no rows. The by column is a pandas Categorical and unused categories are still grouped. Pass observed=True to keep only categories that actually occur.
• Downstream code chokes on MultiPolygon. A dissolve emits whatever the union produces. Split the parts back out with dissolved.explode(index_parts=False), which restores one row per contiguous piece while repeating the aggregated attributes.

Frequently Asked Questions

Does dissolve reproject or otherwise fix my CRS? No. It unions coordinates exactly as they are, in whatever CRS the frame carries, and copies that CRS to the output. Dissolving in EPSG:4326 unions in degree space, which distorts near the poles and breaks outright across the antimeridian, where a group spanning ±180° produces a polygon smeared across the whole globe. Reproject to a suitable projected CRS first — a local UTM zone or an equal-area projection if you will sum areas afterwards — using the transformations described in Coordinate Reference System Transformations.

How is dissolve different from an overlay union? They union different things. dissolve() merges geometries within one layer, grouped by an attribute, and returns fewer rows than it received. overlay(..., how="union") cuts two different layers against each other and returns more rows than either input, one per distinct intersection region. If you are combining a zoning layer with a floodplain layer you want the second, covered in Computing Overlay Union & Difference with GeoPandas.

Can I dissolve everything into a single feature? Yes — call parcels.dissolve() with no by argument and you get one row whose geometry is the union of the entire layer. For geometry alone, parcels.geometry.union_all() returns the same shape as a bare Shapely object without building a frame, which is the cheaper option when you only need an outline or a clipping mask. The distinction between working on one geometry and on a table of them is drawn in Shapely vs GeoPandas: When to Use Each.

Why is my dissolved area smaller than the sum of the parts? Because a union counts overlapping ground once. If two parcels overlap by 200 m², the summed input area double-counts that strip and the dissolved geometry does not. This is the cheapest overlap detector you have: when dissolved.area.sum() is materially below parts.area.sum(), the source layer has duplicate or overlapping features that should be resolved before anything is published.