Vector Overlay Operations

4.4. Vector Overlay Operations#

4.4.1. Learning objectives#

By the end of this section, you should be able to:

  • distinguish an overlay from a spatial join;

  • choose intersection, union, difference, symmetric difference, or identity from a question, and distinguish an overlay from clip;

  • explain how an overlay builds new geometry and why output rows multiply;

  • state the area identity that each operation must satisfy;

  • prepare compatible polygon layers: one projected CRS, valid geometry, and no unnoticed slivers;

  • control which geometry types an overlay keeps; and

  • validate area, identifiers, and geometry quality after an overlay.

Section goal

Use Calgary parks, community districts, and the city boundary to see how an overlay creates new geometries at shared boundaries. The source layers remain unchanged, and every derived area is measured in a projected CRS.

4.4.2. Overlay is not a spatial join#

A spatial join attaches attributes while retaining the geometry of the left feature. An overlay computes new geometry where boundaries intersect, so one input feature may become many output features. Overlays behave like set operations on areas, summarized in Table 4.16:

Table 4.16 Overlay operations, their set logic, and how to interpret each.#

Operation

Set logic

Interpretation

intersection

A AND B

Keep only the shared area of two layers

union

A OR B

Keep all areas and split at every boundary

difference

A AND NOT B

Keep the first layer after removing the second

symmetric_difference

A XOR B

Keep areas belonging to exactly one layer

identity

A, split by B

Keep all of the first layer, split where the second crosses it

Use clip when the question is simply “keep the part inside this mask”; use an overlay when the attributes and boundaries of both polygon layers matter.

4.4.3. How an overlay builds new geometry#

A spatial join leaves geometries untouched; an overlay cuts them. Given two polygon layers A and B, an overlay finds every place where a boundary of A crosses a boundary of B, splits the polygons at those crossings, and sorts the resulting pieces by where they lie: only in A, only in B, or in both. Each operation in the table above is a rule for which pieces to keep. Every kept piece inherits the attributes of the A feature and the B feature it came from, so a piece found only in A has nulls in B’s columns, and the nulls record where the piece came from.

This has three consequences for the output:

  • Rows multiply. A feature crossed by several boundaries becomes several pieces, so the output can have more rows than either input, and an identifier repeats across rows.

  • Areas change. Pieces are new geometries, so their areas are measured after the overlay and never copied from an input column.

  • Columns combine. The output carries both layers’ attributes, and clashing names get _1 and _2 suffixes. difference keeps only A’s columns, because B contributes no pieces.

Because the pieces tile the combined area, every operation obeys an exact area identity, listed in Table 4.17. The identities hold when the features within each layer do not overlap one another, and they are the cheapest way to validate an overlay.

Table 4.17 Area identity that each overlay operation satisfies.#

Operation

Pieces kept

Area identity

intersection

In both A and B

area(A ∩ B) is at most min(area A, area B)

difference

In A only

area(A − B) = area(A) − area(A ∩ B)

union

In A, in B, or in both

area(A ∪ B) = area(A) + area(B) − area(A ∩ B)

symmetric_difference

In exactly one layer

area(A △ B) = area(A) + area(B) − 2 · area(A ∩ B)

identity

Every piece of A, split by B

Total area = area(A)

clip(gdf, mask) is the lightweight cousin of intersection: it keeps the part of gdf inside the mask and returns only gdf’s columns, without the mask’s attributes. The toy layers below are two 4 by 4 squares that share a 2 by 2 corner.

Hide code cell source

 
import math

import geopandas as gpd
import matplotlib.pyplot as plt
from shapely.geometry import box

# Toy layers in metres: two 4 x 4 squares that overlap in a 2 x 2 corner.
layer_a = gpd.GeoDataFrame({'a': ['A']}, geometry=[box(0, 0, 4, 4)])
layer_b = gpd.GeoDataFrame({'b': ['B']}, geometry=[box(2, 2, 6, 6)])

operations = ['intersection', 'union', 'difference', 'symmetric_difference', 'identity']
results = {how: gpd.overlay(layer_a, layer_b, how=how) for how in operations}
results['clip'] = gpd.clip(layer_a, layer_b)

print(f'{"operation":22}{"rows":>5}{"area":>7}  attribute columns')
for name, result in results.items():
    attributes = [column for column in result.columns if column != 'geometry']
    print(f'{name:22}{len(result):>5}{result.area.sum():>7.1f}  {attributes}')

area_a, area_b = layer_a.area.sum(), layer_b.area.sum()
shared = results['intersection'].area.sum()
identities = {
    'difference': area_a - shared,
    'union': area_a + area_b - shared,
    'symmetric_difference': area_a + area_b - 2 * shared,
    'identity': area_a,
}
print()
for name, expected in identities.items():
    print(f'{name:22} area identity holds: {math.isclose(results[name].area.sum(), expected)}')

print()
display(results['union'].drop(columns='geometry'))

operation              rows   area  attribute columns
intersection              1    4.0  ['a', 'b']
union                     3   28.0  ['a', 'b']
difference                1   12.0  ['a']
symmetric_difference      2   24.0  ['a', 'b']
identity                  2   16.0  ['a', 'b']
clip                      1    4.0  ['a']

difference             area identity holds: True
union                  area identity holds: True
symmetric_difference   area identity holds: True
identity               area identity holds: True
a b
0 A B
1 A NaN
2 NaN B

The areas confirm the identities. Each square has area 16 and they share an area of 4, so union is 16 + 16 − 4 = 28, difference is 16 − 4 = 12, symmetric_difference is 16 + 16 − 2 × 4 = 24, and identity returns A’s own 16. In the union table, row 0 has both a and b (the shared corner), row 1 only a, and row 2 only b: the nulls record which layer each piece came from.

../_images/overlay_toy_operations.png

Fig. 4.7 The five overlay operations and clip applied to two overlapping squares; each colour is one output row and the dashed lines are the input boundaries.#

In Fig. 4.7, each colour is one output row, and the dashed lines are the input boundaries. union returns three pieces, and a null in a or b marks a piece that exists in only one layer. difference keeps A’s shape minus the shared corner and only A’s column. symmetric_difference keeps the two non-shared pieces, and identity keeps all of A split at B’s boundary, so its pieces add up to A. clip returns the same shape as intersection but without B’s attribute. Overlay is also order-sensitive: difference removes the second layer from the first, and swapping the arguments answers a different question.

4.4.3.1. Which geometry types an overlay keeps#

Two polygons that only touch share a line or a point, not an area. An overlay can therefore produce a piece of lower dimension, such as a LineString along a shared edge. By default (keep_geom_type=True) GeoPandas [Jordahl et al., 2020] drops pieces whose geometry type differs from the input, which suits area questions because those pieces have zero area. With keep_geom_type=False they are kept, so no piece is dropped silently, and the row count rises accordingly.

Hide code cell source

 
import warnings

# A neighbour that shares only the edge x = 4 with layer A.
neighbour = gpd.GeoDataFrame({'b': ['C']}, geometry=[box(4, 0, 8, 4)])

for keep in [True, False]:
    with warnings.catch_warnings():
        warnings.simplefilter('ignore')
        touching = gpd.overlay(layer_a, neighbour, how='intersection', keep_geom_type=keep)
    print(f'keep_geom_type={keep!s:5} rows={len(touching)}  types={touching.geom_type.tolist()}')

keep_geom_type=True  rows=0  types=[]
keep_geom_type=False rows=1  types=['LineString']

4.4.4. Load the layers#

The next cell loads the three polygon layers [Open Calgary: The City of Calgary Open Data Portal, 2026]: a deterministic 300-feature sample of parks, the community districts, and the city boundary. It also creates metric copies, marked with the suffix _m, in a local UTM CRS [Snyder, 1987].

Hide code cell source

 
from pathlib import Path

import geopandas as gpd
import pandas as pd

DATA_DIR = Path('../data') if Path('../data').exists() else Path('data')

parks = gpd.read_file(DATA_DIR / 'Parks_Sites_20261019.geojson')
parks = (
    parks.loc[parks.geometry.notna()]
    .sort_values('asset_cd', kind='stable')
    .head(300)
    .copy()
)

communities_raw = pd.read_csv(DATA_DIR / 'Community_District_Boundaries_20240620.csv')
communities = gpd.GeoDataFrame(
    communities_raw,
    geometry=gpd.GeoSeries.from_wkt(communities_raw['MULTIPOLYGON']),
    crs='EPSG:4326',
)

city_boundary = gpd.read_file(DATA_DIR / 'City Boundary_20240620.gpkg')
metric_crs = parks.estimate_utm_crs()
parks_m = parks.to_crs(metric_crs)
communities_m = communities.to_crs(metric_crs)
city_m = city_boundary.to_crs(metric_crs)

4.4.5. Prepare compatible inputs#

An overlay is only as reliable as its inputs. Check four conditions before running one:

  • One CRS. GeoPandas warns on a CRS mismatch but still computes, and the result is meaningless. Assert that the CRSs are equal instead of relying on the warning.

  • A projected CRS. Pieces are measured after the overlay, so areas need metres, as in S02.

  • Valid geometry. An invalid ring can make the overlay fail or return a misleading result, as the bow-tie in S01 showed. Repair the feature instead of dropping it.

  • No unnoticed slivers. Layers digitised separately rarely agree exactly on a shared boundary. An offset of a millimetre turns a shared edge into a thin overlap, and the overlay returns it as a real polygon of almost zero area. Slivers inflate row counts and are not real features, so detect them by area or thinness and document the threshold you use.

The toy layer below should share an edge with layer A but is offset by one millimetre.

Hide code cell source

 
# Should share the edge x = 4 with layer A, but is off by 0.001 m.
offset_neighbour = gpd.GeoDataFrame({'d': ['D']}, geometry=[box(3.999, 0, 8, 4)])

sliver = gpd.overlay(layer_a, offset_neighbour, how='intersection')
minx, miny, maxx, maxy = sliver.total_bounds
print(f'intersection rows: {len(sliver)}')
print(f'sliver area: {sliver.area.sum():.4f} m², width: {maxx - minx:.3f} m, height: {maxy - miny:.1f} m')
print(f'valid: {bool(sliver.geometry.is_valid.all())}  (nothing flags it as a problem)')

intersection rows: 1
sliver area: 0.0040 m², width: 0.001 m, height: 4.0 m
valid: True  (nothing flags it as a problem)

4.4.5.1. Prepare the Calgary layers#

A full parks layer can be expensive to overlay against every district while learning, so this example uses a deterministic 300-feature sample. The sample is for understanding the operation, not for estimating Calgary-wide park area. One district polygon has a self-intersecting ring, which would make overlay results unreliable, so the repair step below fixes it before any overlay runs. Section S06 covers this repair pattern in detail.

Hide code cell source

 
for name, layer in {
    'parks': parks_m,
    'communities': communities_m,
    'city boundary': city_m,
}.items():
    print(name, layer.crs, layer.geometry.geom_type.value_counts().to_dict())
    print('  invalid:', int((~layer.geometry.is_valid).sum()))

assert parks_m.crs == communities_m.crs == city_m.crs
assert parks_m.geometry.notna().all()
assert communities_m.geometry.notna().all()

# One district ring self-intersects. Repair it explicitly rather than dropping the
# feature, and confirm the layer is fully valid before overlaying.
from shapely.validation import make_valid

communities_m['geometry'] = communities_m.geometry.map(make_valid)
assert communities_m.geometry.is_valid.all()

parks EPSG:32611 {'MultiPolygon': 300}
  invalid: 0
communities EPSG:32611 {'MultiPolygon': 313}
  invalid: 1
city boundary EPSG:32611 {'MultiPolygon': 1}
  invalid: 0

As in S01, every layer is stored as MultiPolygon. The invalid ring is in the community layer (1 of 313 districts) and the park sample has none; the repair keeps all 313 districts instead of dropping the faulty one.

4.4.6. Intersection: shared park and district areas#

A park can cross more than one community district, so the intersection output may contain several rows for one park. Keep source identifiers and calculate the area of each resulting piece rather than treating every output row as a new park.

Hide code cell source

 
park_fields = parks_m[['site_name', 'asset_cd', 'geometry']]
district_fields = communities_m[['COMM_CODE', 'NAME', 'geometry']]

park_district_parts = gpd.overlay(
    park_fields,
    district_fields,
    how='intersection',
    keep_geom_type=False,
)
park_district_parts['shared_area_ha'] = (
    park_district_parts.geometry.area / 10_000
)

print('intersection rows:', len(park_district_parts))
print(park_district_parts[['site_name', 'NAME', 'shared_area_ha']].head())
assert (park_district_parts['shared_area_ha'] >= 0).all()
assert park_district_parts.geometry.notna().all()

intersection rows: 333
      site_name          NAME  shared_area_ha
0  01B GREENWAY           01B        1.759191
1  01B GREENWAY  SCENIC ACRES        0.706585
2  01B GREENWAY           01B       11.205214
3           NaN           01K        7.241670
4           NaN           02C        0.697641

Check the result against the inputs. The shared areas cannot exceed the sampled park area, and a park identifier now repeats once per piece.

Hide code cell source

 
park_area_ha = park_fields.geometry.area.sum() / 10_000
shared_ha = park_district_parts['shared_area_ha'].sum()
print(f'sampled park area: {park_area_ha:,.1f} ha')
print(f'area inside a district: {shared_ha:,.1f} ha ({shared_ha / park_area_ha:.1%})')

pieces = park_district_parts.groupby('asset_cd').size()
print('park identifiers:', park_fields['asset_cd'].nunique(), '-> pieces:', len(park_district_parts))
print('parks split into several pieces:', int((pieces > 1).sum()))

assert shared_ha <= park_area_ha * (1 + 1e-9)

sampled park area: 1,158.9 ha
area inside a district: 1,095.3 ha (94.5%)
park identifiers: 300 -> pieces: 333
parks split into several pieces: 27

The 300 sampled parks become 333 pieces, and 27 parks are split across districts (at least 60 of the 333 pieces come from those 27 parks, so some parks are cut into three or more pieces). About 94.5% of the sampled park area lies inside a district, leaving roughly 64 ha, the 5.5% that falls in no district. Two of the rows listed earlier have no site_name, so park names can be missing and asset_cd is the identifier to keep.

4.4.7. Difference and union#

difference(parks, city_boundary) keeps the sampled park pieces outside the city boundary. Reversing the arguments answers a different question. union creates a combined partition and carries attributes from both layers; it is useful for partitioning space, but it usually creates more rows than either input.

Hide code cell source

 
city_mask = city_m[['geometry']]
outside_city = gpd.overlay(
    park_fields, city_mask, how='difference', keep_geom_type=False
)
inside_city = gpd.overlay(
    park_fields, city_mask, how='intersection', keep_geom_type=False
)

print('park pieces inside city:', len(inside_city))
print('park pieces outside city:', len(outside_city))

# Use small subsets to keep the union demonstration quick.
union_example = gpd.overlay(
    communities_m.head(10)[['COMM_CODE', 'geometry']],
    parks_m.head(25)[['asset_cd', 'geometry']],
    how='union',
    keep_geom_type=False,
)
print('union rows:', len(union_example))

input_area = park_fields.geometry.area.sum()
outside_area = outside_city.geometry.area.sum()
# A difference cannot add area to the first input.
assert outside_area <= input_area

park pieces inside city: 299
park pieces outside city: 2
union rows: 35

Intersection and difference against the same mask split each park into complementary pieces, so their areas should add up to the input area.

Hide code cell source

 
inside_ha = inside_city.geometry.area.sum() / 10_000
outside_ha = outside_area / 10_000
print(f'inside + outside: {inside_ha + outside_ha:,.1f} ha, input: {input_area / 10_000:,.1f} ha')
print(f'relative difference: {(inside_ha + outside_ha) * 10_000 / input_area - 1:.1e}')

inside + outside: 1,158.9 ha, input: 1,158.9 ha
relative difference: -2.6e-15

The two results are complementary to rounding error (a relative difference of −2.6 × 10⁻¹⁵), so the area identity for splitting by a mask holds to machine precision. Only 2 of the 301 pieces lie outside the city, so almost all the sampled parks are inside it.

4.4.8. Symmetric difference and identity#

symmetric_difference keeps the areas that belong to exactly one layer — useful for questions such as “which parts of these two planning zones do not overlap?”. identity keeps the entire first layer but splits it wherever the second layer’s boundaries cross, attaching the second layer’s attributes to the pieces inside it. Identity answers “give me all of layer A, tagged by layer B”. The example below uses a small district subset, so most sampled parks fall outside those districts; with the full district layer, the pieces would carry district attributes wherever they overlap.

Hide code cell source

 
sym_diff = gpd.overlay(
    communities_m.head(2)[['COMM_CODE', 'geometry']],
    city_mask,
    how='symmetric_difference',
    keep_geom_type=False,
)
print('symmetric difference rows:', len(sym_diff))

identity_example = gpd.overlay(
    park_fields.head(25),
    communities_m.head(10)[['COMM_CODE', 'geometry']],
    how='identity',
    keep_geom_type=False,
)
print('identity rows:', len(identity_example))
print('identity keeps park rows with and without a district:',
      int(identity_example['COMM_CODE'].isna().sum()), 'without,',
      int(identity_example['COMM_CODE'].notna().sum()), 'with')

symmetric difference rows: 1
identity rows: 25
identity keeps park rows with and without a district: 25 without, 0 with
../_images/overlay_park_pieces.png

Fig. 4.8 Sampled Calgary park pieces after intersection with community districts, colored by area in hectares; district boundaries in grey. Data: City of Calgary Open Data.#

In Fig. 4.8, each shaded polygon is a piece of a sampled park cut at a district boundary. Almost all pieces are small (dark purple, under about 50 ha), while one large piece near the eastern edge, close to 275 ha on the colour scale, stretches the colours for everything else. The pieces cluster in a few parts of the city instead of spreading evenly, because the sample is the first 300 parks by asset code and not a random draw, which is another reason not to estimate Calgary-wide park area from it.

4.4.9. Choosing the right operation#

  • Only the shared area matters → intersection

  • A complete partition of both layers → union

  • Everything in A except B → difference

  • Areas unique to either layer → symmetric_difference

  • All of A, tagged by B → identity

  • Only the extent of a mask matters, not its attributes → clip

Overlays assume valid input polygons and a shared projected CRS. Invalid rings, slivers, and degree-based area calculations are the three most common ways an overlay produces a plausible-looking but wrong result.

Slivers add a fourth cause: near-coincident boundaries create thin pieces that look like real features. When the layers should share a boundary, check the smallest piece areas before counting rows.

4.4.10. Summary#

Overlay is the right tool when shared or excluded area matters. Use one projected CRS, repair invalid polygons, preserve source identifiers, expect feature splitting, and watch for slivers and zero-area pieces. Validate with the area identities, row counts, geometry validity, and null attributes, which record where each piece came from. Sections S01–S04 now form a complete vector workflow: identify the data type, measure or transform it, join information, and combine areas. Continue with S05 to apply the full workflow to two Calgary questions.