Overlay Operations in GeoPandas: Union, Intersection, Difference Explained

Problem statement

Overlaying two polygon layers cuts them along each other's boundaries and keeps the pieces a question asks for:

• Which parcels fall inside a flood zone?
• What areas are covered by either of two planning layers?
• What land remains after removing protected areas?

This is where GeoPandas overlay operations are useful. They create a new layer by cutting and combining the boundaries from two inputs.

This is also where users often get stuck. overlay() is not the same as:

• sjoin(), which matches features based on a spatial relationship
• clip(), which trims one layer to a boundary
• direct Shapely geometry methods used on single geometries

This page focuses on choosing the correct geopandas.overlay() operation and using it correctly for polygon and multipolygon analysis.

Quick answer

Use geopandas.overlay() when you need a new layer built from two existing layers by cutting and combining their boundaries. In this page, the focus is polygon and multipolygon workflows.

• how="intersection" keeps only overlapping areas
• how="union" keeps all areas from both layers and splits polygons where they overlap
• how="difference" keeps only the parts of the first layer that are not covered by the second

In most cases, both layers should:

• use the same CRS
• contain polygon or multipolygon geometries
• have valid geometries

Choosing the right operation

Rutland parcels against flood zones: intersection, union, difference

Step 1: Load the two spatial layers

This page uses real layers for Rutland (set up under Code examples):

• parcels.gpkg: HM Land Registry's 21,253 freehold parcels
• flood_zones.gpkg: the Environment Agency's Flood Zone 3 — land with a 1% or greater annual chance of river flooding — from its Flood Map for Planning

import geopandas as gpd

parcels = gpd.read_file("data/parcels.gpkg")
flood_zones = gpd.read_file("data/flood_zones.gpkg")

print(len(parcels), "parcels:", list(parcels.columns))
print(len(flood_zones), "flood zone polygons:", list(flood_zones.columns))

21253 parcels: ['INSPIREID', 'geometry']
684 flood zone polygons: ['flood_zone', 'flood_source', 'geometry']

For the workflows in this page, both inputs should usually be polygon or multipolygon layers.

Step 2: Check CRS before overlay

Overlay analysis should be done in the same coordinate reference system.

print(parcels.crs, flood_zones.crs)

EPSG:27700 EPSG:27700

If they do not match, reproject one layer:

flood_zones = flood_zones.to_crs(parcels.crs)

If area or distance matters, use a projected CRS rather than a geographic CRS such as EPSG:4326. Both layers here are already in British National Grid (EPSG:27700), which is in metres.

Step 3: Inspect geometry types and validity

Overlay works best when both layers contain valid polygon geometries.

print(parcels.geom_type.unique().tolist(), flood_zones.geom_type.unique().tolist())
print((~parcels.geometry.is_valid).sum(), (~flood_zones.geometry.is_valid).sum())

['Polygon'] ['MultiPolygon']
0 0

Remove empty geometries if needed:

parcels = parcels[~parcels.geometry.is_empty & parcels.geometry.notnull()].copy()
flood_zones = flood_zones[~flood_zones.geometry.is_empty & flood_zones.geometry.notnull()].copy()

If invalid geometries are present, repair them before overlay:

parcels["geometry"] = parcels.geometry.make_valid()
flood_zones["geometry"] = flood_zones.geometry.make_valid()

make_valid() needs Shapely 2 and GeoPandas 0.12 or later. Both layers here are valid; the SSSI layer in Example 3 has one invalid polygon, which is where the repair earns its keep.

Step 4: Run geopandas.overlay() with how="intersection"

Use intersection when you only want areas shared by both layers.

parcel_flood_overlap = gpd.overlay(parcels, flood_zones, how="intersection")

print(len(parcel_flood_overlap), "pieces from", parcel_flood_overlap["INSPIREID"].nunique(), "parcels")
print(list(parcel_flood_overlap.columns))
print(f"{parcel_flood_overlap.area.sum() / 1e4:,.0f} ha of parcel land in Flood Zone 3")

2670 pieces from 1737 parcels
['INSPIREID', 'flood_zone', 'flood_source', 'geometry']
2,594 ha of parcel land in Flood Zone 3

This is a common workflow for finding the parts of parcels that lie inside hazard zones. Each piece carries the parcel's attributes and the flood zone's.

Step 5: Run geopandas.overlay() with how="union"

Use union when you need all coverage from both layers, split into unique pieces.

parcel_flood_union = gpd.overlay(parcels, flood_zones, how="union")

print(len(parcel_flood_union), "pieces")
print("parcel only:", parcel_flood_union["flood_zone"].isna().sum(),
      "| flood zone only:", parcel_flood_union["INSPIREID"].isna().sum(),
      "| both:", (parcel_flood_union["flood_zone"].notna() & parcel_flood_union["INSPIREID"].notna()).sum())

24407 pieces
parcel only: 21135 | flood zone only: 602 | both: 2670

This keeps:

• overlap areas
• parcel-only areas
• flood-zone-only areas — the parts of the flood zone no registered parcel covers

The output is often split into many smaller polygons because GeoPandas cuts the layer wherever boundaries intersect.

Step 6: Run geopandas.overlay() with how="difference"

Use difference when you want the first layer minus the second layer.

parcels_outside_flood = gpd.overlay(parcels, flood_zones, how="difference")
flood_outside_parcels = gpd.overlay(flood_zones, parcels, how="difference")

print(f"parcels minus flood zone: {len(parcels_outside_flood):,} features, "
      f"{parcels_outside_flood.area.sum() / 1e4:,.0f} ha")
print(f"flood zone minus parcels: {len(flood_outside_parcels):,} features, "
      f"{flood_outside_parcels.area.sum() / 1e4:,.0f} ha")

parcels minus flood zone: 21,135 features, 40,304 ha
flood zone minus parcels: 602 features, 2,661 ha

Order matters here:

• overlay(parcels, flood_zones, how="difference") means parcels minus flood zones
• overlay(flood_zones, parcels, how="difference") means flood zones minus parcels

These are different results, with different attributes: each keeps only the first layer's columns.

Step 7: Review the output attributes

Overlay outputs combine columns from both inputs.

For intersection:

• output features contain attributes from both parcels and flood zones

For union:

• some polygons come only from one layer
• fields from the other layer will be NaN

print(parcel_flood_union[["INSPIREID", "flood_zone", "flood_source"]].sample(6, random_state=1).to_string())

        INSPIREID flood_zone flood_source
1522   32214191.0        FZ3        river
5065   57714461.0        NaN          NaN
19836  61371413.0        NaN          NaN
14078  31997954.0        NaN          NaN
16469  31566599.0        NaN          NaN
18709  32075655.0        NaN          NaN

Look at INSPIREID too: an integer column in the parcels, it comes out of the union as floats (32214191.0), because pandas' default integer type cannot hold the NaN of the flood-zone-only pieces. Cast it back with .astype("Int64") before joining on it.

If both inputs have a field with the same name, GeoPandas suffixes them to avoid collisions. Overlaying the flood map's Zone 2 on its Zone 3 — two layers with the same schema — shows it:

zone2 = gpd.read_file("data/flood_zones_all.gpkg").query("flood_zone == 'FZ2'")
print(list(gpd.overlay(zone2, flood_zones, how="intersection").columns))

['flood_zone_1', 'flood_source_1', 'flood_zone_2', 'flood_source_2', 'geometry']

Step 8: Save the result to a new file

GeoPackage is usually the safest format for modern workflows.

from pathlib import Path
Path("output").mkdir(exist_ok=True)

parcel_flood_overlap.to_file("output/parcel_flood_overlap.gpkg", layer="intersection", driver="GPKG")
parcel_flood_union.to_file("output/parcel_flood_union.gpkg", layer="union", driver="GPKG")
parcels_outside_flood.to_file("output/parcels_outside_flood.gpkg", layer="difference", driver="GPKG")

You can also export to Shapefile or GeoJSON, but GeoPackage handles field names and schema more reliably.

Overlaying parcels, flood zones and protected sites

The data used on this page

HM Land Registry's INSPIRE parcels for Rutland, the Environment Agency's Flood Map for Planning (Flood Zones 2 and 3, rivers and sea) from its WFS, and Natural England's Sites of Special Scientific Interest from its open ArcGIS service — the last two requested for the box around Rutland, all in British National Grid.

import io
import zipfile
from pathlib import Path

import geopandas as gpd
import pyogrio
import requests

pyogrio.set_gdal_config_options({"GML_DOWNLOAD_SCHEMA": "NO"})
Path("data").mkdir(exist_ok=True)
BOX = (479000, 292000, 507000, 319000)

if not Path("data/parcels.gpkg").exists():
    reply = requests.Session().get("https://use-land-property-data.service.gov.uk/datasets/inspire/download/"
                                   "Rutland.zip", timeout=600)
    reply.raise_for_status()
    zipfile.ZipFile(io.BytesIO(reply.content)).extract("Land_Registry_Cadastral_Parcels.gml", "downloads")
    (gpd.read_file("downloads/Land_Registry_Cadastral_Parcels.gml")[["INSPIREID", "geometry"]]
        .to_file("data/parcels.gpkg"))

if not Path("data/flood_zones_all.gpkg").exists():
    reply = requests.get("https://environment.data.gov.uk/spatialdata/flood-map-for-planning-flood-zones/wfs",
                         timeout=300, params={
                             "service": "WFS", "version": "2.0.0", "request": "GetFeature",
                             "typeNames": "dataset-04532375-a198-476e-985e-0579a0a11b47:Flood_Zones_2_3_Rivers_and_Sea",
                             "bbox": ",".join(map(str, BOX)) + ",urn:ogc:def:crs:EPSG::27700",
                             "outputFormat": "application/json", "count": 100000})
    reply.raise_for_status()
    (gpd.read_file(io.BytesIO(reply.content)).set_crs(27700, allow_override=True)
        [["flood_zone", "flood_source", "geometry"]].to_file("data/flood_zones_all.gpkg"))
all_zones = gpd.read_file("data/flood_zones_all.gpkg")
all_zones[all_zones["flood_zone"] == "FZ3"].to_file("data/flood_zones.gpkg")

if not Path("data/sssi.gpkg").exists():
    reply = requests.get("https://services.arcgis.com/JJzESW51TqeY9uat/arcgis/rest/services/SSSI_England/"
                         "FeatureServer/0/query", timeout=120, params={
                             "geometry": ",".join(map(str, BOX)), "geometryType": "esriGeometryEnvelope",
                             "inSR": 27700, "spatialRel": "esriSpatialRelIntersects",
                             "outFields": "NAME", "outSR": 27700, "f": "geojson"})
    reply.raise_for_status()
    gpd.read_file(io.BytesIO(reply.content)).to_file("data/sssi.gpkg")

print({name: len(gpd.read_file(f"data/{name}.gpkg")) for name in ["parcels", "flood_zones_all", "flood_zones", "sssi"]},
      "| zones:", all_zones["flood_zone"].value_counts().to_dict())

{'parcels': 21253, 'flood_zones_all': 2209, 'flood_zones': 684, 'sssi': 44} | zones: {'FZ2': 1525, 'FZ3': 684}

Example 1: Overlay two polygon layers with intersection

import geopandas as gpd

parcels = gpd.read_file("data/parcels.gpkg").to_crs("EPSG:27700")
flood_zones = gpd.read_file("data/flood_zones.gpkg").to_crs("EPSG:27700")

parcels["geometry"] = parcels.geometry.make_valid()
flood_zones["geometry"] = flood_zones.geometry.make_valid()

result = gpd.overlay(parcels, flood_zones, how="intersection")

print(result.shape)
print(result.columns.tolist())

(2670, 4)
['INSPIREID', 'flood_zone', 'flood_source', 'geometry']

Example 2: Use union to split and combine all coverage

The two designations that most often constrain rural land, flood risk and nature conservation, combined into one layer of unique pieces:

constraints = gpd.read_file("data/sssi.gpkg").to_crs(flood_zones.crs)
print("invalid SSSI polygons before repair:", (~constraints.is_valid).sum())
constraints["geometry"] = constraints.geometry.make_valid()

combined = gpd.overlay(flood_zones, constraints, how="union", keep_geom_type=True)
both = combined["flood_zone"].notna() & combined["NAME"].notna()
print(f"{len(combined):,} pieces; {both.sum()} are both Flood Zone 3 and SSSI, "
      f"{combined[both].area.sum() / 1e4:,.0f} ha:")
print(combined[both].dissolve("NAME").area.div(1e4).round(1).sort_values(ascending=False).head(4).to_string())

invalid SSSI polygons before repair: 1
769 pieces; 66 are both Flood Zone 3 and SSSI, 960 ha:
NAME
Rutland Water SSSI               770.8
Eye Brook Reservoir SSSI         163.6
Seaton Meadows SSSI                7.9
Empingham Marshy Meadows SSSI      4.5

Example 3: Use difference to remove excluded areas

The parcels left after taking out both constraints:

excluded = gpd.GeoDataFrame(geometry=[combined.union_all()], crs=combined.crs)
eligible_land = gpd.overlay(parcels, excluded, how="difference")
print(f"{len(parcels):,} parcels, {parcels.area.sum() / 1e4:,.0f} ha -> "
      f"{len(eligible_land):,} parcels, {eligible_land.area.sum() / 1e4:,.0f} ha outside both")

21,253 parcels, 42,899 ha -> 21,106 parcels, 38,905 ha outside both

Dissolving the constraints into one shape first keeps the difference simple: one subtraction per parcel rather than one per parcel per constraint.

Example 4: Export overlay output

eligible_land.to_file(
    "output/eligible_land.gpkg",
    layer="eligible_land",
    driver="GPKG"
)

check = gpd.read_file("output/eligible_land.gpkg", layer="eligible_land")
print(check.crs)
print(check.geom_type.unique().tolist())
print(check.columns.tolist())

EPSG:27700
['MultiPolygon']
['INSPIREID', 'geometry']

Overlay compared with clip and join, and input order

What overlay does compared to clip and spatial join

Overlay creates new geometries by cutting and combining polygon boundaries from both layers.

Use clip() when you want to trim one layer to a mask boundary, but you do not need attributes from both layers.

Use sjoin() when you want to attach attributes based on spatial relationships like intersects or within, without cutting the geometries.

Use overlay() when the geometry itself must change as part of the analysis.

When to use intersection

Use intersection when only shared areas matter.

Example:

• parcel area inside a flood zone
• habitat overlap with zoning districts

This is often the right choice for regulated-area analysis.

When to use union

Use union when you need a complete partition of space from both layers.

Example:

• combining zoning and environmental policy layers
• preparing polygons for area summaries by all combinations of attributes

This creates the most detailed output and often the most features.

When to use difference

Use difference when removing one layer from another.

Example:

• buildable land minus wetlands
• service area minus excluded zones

The first input defines what you are keeping.

Why input order matters

Input order matters most for difference, because A - B is not the same as B - A.

In practice, define the first layer as the one you want to keep or analyze, and the second as the layer to compare against or subtract.

Polygon scope, invalid shapes and union gaps

Polygon scope for this page

This page focuses on polygon and multipolygon overlays. If your inputs include mixed geometry types such as points or lines, handle those workflows separately.

Invalid geometries

Self-intersections, duplicate rings, and broken polygons can cause overlay failures or incorrect output. Validate and repair geometries first.

If .make_valid() is not available in your environment, use a geometry repair method supported by your GeoPandas and Shapely versions before running overlay.

CRS mismatch

If the CRS does not match, results may be wrong or the operation may fail. Always align CRS before overlay. If area is important, use a projected CRS.

Performance on large datasets

Overlay can be slow and memory-intensive for large polygon layers. If possible:

• filter to the study area first
• keep only necessary columns
• process in smaller subsets

Missing values in union output

With how="union", some output polygons exist only in one input layer. Those features will have null values for fields from the other layer. This is expected.

Internal links

• Python for GIS: What It Is and When to Use It
• GeoPandas Basics: Working with Spatial Data in Python
• Coordinate Reference Systems (CRS) Explained for Python GIS

FAQ

What is the difference between overlay() and sjoin() in GeoPandas?

overlay() creates new geometries by cutting and combining boundaries. sjoin() keeps the original geometries and only transfers attributes based on spatial relationships.

When should I use intersection instead of clip?

Use intersection when you want attributes from both layers and need the true overlap polygons. Use clip when you only want to trim one layer to a mask boundary.

Why does difference return unexpected results?

The most common reasons are:

• the inputs are in different CRS
• geometries are invalid
• input order is reversed

Remember that difference keeps only the parts of the first layer outside the second.

Do both layers need the same CRS for GeoPandas overlay operations?

Yes. In practice, both layers should be in the same CRS before running overlay. Reproject one layer with .to_crs() if needed.

What does how="symmetric_difference" do in gpd.overlay()?

It keeps the parts of each layer that do not overlap the other, dropping the shared area entirely. It is effectively the union minus the intersection, useful for finding where two layers disagree.

How do I merge overlapping polygons within a single layer instead of between two layers?

overlay() works between two layers, so to dissolve overlaps inside one layer use gdf.union_all() for a single combined geometry, or dissolve() to merge by an attribute. Overlay is the wrong tool when there is only one input layer.

Why is my overlay output empty even though the layers clearly overlap on a map?

The most common cause is mismatched CRS, where the two layers are in different coordinate systems so their geometries never actually intersect numerically. Confirm both with print(gdf.crs) and reproject one with to_crs() before running overlay().

How do I compute the area of each parcel that falls inside a flood zone?

Run gpd.overlay(parcels, flood_zones, how="intersection"), then in a projected CRS add a column with result["area_m2"] = result.geometry.area. Make sure both inputs are in a metre-based projected CRS first, otherwise .area returns meaningless square-degree values.