Intersection, union, and difference—the foundation of spatial analysis
2026-02-27
After the 2013 Alberta floods, the provincial government faced a spatial question that was simple to state and complex to compute: which properties were inside the floodway (where reconstruction was prohibited), which were in the flood fringe (where reconstruction was permitted with conditions), and which in each case were already developed versus vacant? The answer required intersecting four polygon layers: the floodway boundary, the flood fringe boundary, municipal parcel boundaries, and existing land-use zoning. The result of that intersection — a new set of polygons carrying attributes inherited from all four input layers — was the map that determined whether individual landowners would be bought out, required to elevate, or left undisturbed.
Polygon overlay is the operation that combines two polygon layers to produce a third whose geometry is the intersection, union, or difference of the input geometries, and whose attribute table carries columns from both inputs. It is one of the most powerful and most computationally demanding operations in vector GIS. The challenge is that real polygon datasets have vertices that are numerically close but not identical, edges that nearly coincide, and topological configurations that require careful handling to avoid slivers and gaps in the output. This model derives the line-intersection test that underlies all polygon overlay, implements the Weiler-Atherton polygon clipping algorithm, and discusses the numerical precision issues that make robust implementations difficult.
Which parcels are both inside the floodplain AND zoned for development?
Overlay operations combine two spatial datasets using set operations:
Intersection (AND): Areas in both datasets - Forest ∩ Protected areas = Protected forest - Floodplain ∩ Development zones = Flood-risk development areas - Species habitat ∩ Private land = Conservation priority areas
Union (OR): Areas in either dataset - Parks ∪ Preserves = All protected areas - Urban ∪ Agricultural = Developed land - Fire risk zones ∪ Drought zones = Combined hazard areas
Difference (NOT): Areas in first but not second - Total land Water bodies = Land area only - Proposed development Wetlands = Developable area - National forest Wilderness areas = Multiple-use forest
The mathematical question: How do we compute the intersection, union, or difference of two polygons?
Polygons as sets of points:
Polygon A = \{(x,y) : \text{point inside polygon A}\}
Boolean operations:
A \cap B = \{(x,y) : (x,y) \in A \text{ AND } (x,y) \in B\}
A \cup B = \{(x,y) : (x,y) \in A \text{ OR } (x,y) \in B\}
A \setminus B = \{(x,y) : (x,y) \in A \text{ AND } (x,y) \notin B\}
Result: New polygon(s) representing the combined area.
Special case: Clip polygon A to window B (intersection).
Examples: - Clip state boundaries to study area rectangle - Clip satellite image to county boundary - Clip parcels to watershed boundary
Challenge: Find vertices of resulting polygon.
Required pieces: 1. Vertices from A inside B 2. Vertices from B inside A 3. Edge-edge intersection points 4. Correct traversal order
Foundation: Must detect where polygon edges cross.
Segments: P_1P_2 and Q_1Q_2
Parametric form:
P(s) = P_1 + s(P_2 - P_1), \quad s \in [0,1] Q(t) = Q_1 + t(Q_2 - Q_1), \quad t \in [0,1]
Intersection: Find s, t where P(s) = Q(t):
P_1 + s(P_2 - P_1) = Q_1 + t(Q_2 - Q_1)
Solve 2×2 system:
Let \mathbf{d}_P = P_2 - P_1, \mathbf{d}_Q = Q_2 - Q_1, \mathbf{r} = P_1 - Q_1
s \mathbf{d}_P - t \mathbf{d}_Q = -\mathbf{r}
Matrix form:
\begin{pmatrix} d_{Px} & -d_{Qx} \\ d_{Py} & -d_{Qy} \end{pmatrix} \begin{pmatrix} s \\ t \end{pmatrix} = \begin{pmatrix} -r_x \\ -r_y \end{pmatrix}
Determinant:
D = d_{Px} \cdot d_{Qy} - d_{Py} \cdot d_{Qx}
If D = 0: segments are parallel (no intersection or infinite intersections)
Solution:
s = \frac{(Q_1 - P_1) \times (Q_2 - Q_1)}{(P_2 - P_1) \times (Q_2 - Q_1)}
t = \frac{(Q_1 - P_1) \times (P_2 - P_1)}{(P_2 - P_1) \times (Q_2 - Q_1)}
Where \times is 2D cross product: (x_1, y_1) \times (x_2, y_2) = x_1 y_2 - y_1 x_2
Intersection exists if: $0 s $ AND $0 t $
Intersection point:
I = P_1 + s(P_2 - P_1)
Input: Two simple polygons A (subject) and B (clip)
Output: Intersection polygon A \cap B
Algorithm steps:
1. Find all edge intersections
For each edge of A and each edge of B: - Test for intersection - Record intersection points - Label as “entering” or “leaving”
Entering: Edge goes from outside B to inside B
Leaving: Edge goes from inside B to outside B
2. Build vertex list
For polygon A: - Original vertices - Intersection points (inserted in traversal order) - Mark each vertex: inside/outside B
For polygon B: - Original vertices
- Intersection points - Mark each vertex: inside/outside A
3. Trace boundary
Start at an intersection point marked “entering”: - Follow A edges while inside B - At “leaving” intersection, switch to following B edges - At next “entering” intersection, switch back to A - Continue until returning to start
Output: Sequence of vertices forming intersection polygon.
Union = all area covered by either polygon
Weiler-Atherton modification:
Instead of following edges “inside” clip polygon, follow edges “outside”: - Start at “leaving” intersection - Follow edges while outside the other polygon - Switch at intersections
Multiple components: Union may produce multiple disconnected polygons (if originals don’t overlap).
Difference A \setminus B = area in A but not in B
Implementation:
Follow A edges normally, but when inside B: - Switch to following B boundary in reverse direction - This “cuts out” the B region
Equivalent: A \setminus B = A \cap \overline{B} (intersection with complement)
Problem: Find intersection of: - Rectangle A: corners at (0,0), (4,0), (4,3), (0,3) - Rectangle B: corners at (2,1), (6,1), (6,4), (2,4)
Step 1: Find edge intersections
Edges of A: - Bottom: (0,0) to (4,0) - Right: (4,0) to (4,3) - Top: (4,3) to (0,3) - Left: (0,3) to (0,0)
Edges of B: - Bottom: (2,1) to (6,1) - Right: (6,1) to (6,4) - Top: (6,4) to (2,4) - Left: (2,4) to (2,1)
Test all pairs:
A right (4,0) \to (4,3) intersects B bottom (2,1) \to (6,1): - At point (4,1) ✓
A right (4,0) \to (4,3) intersects B top (6,4) \to (2,4): - No (doesn’t reach y=4)
A top (4,3) \to (0,3) intersects B left (2,4) \to (2,1): - At point (2,3) ✓
Intersections found: (4,1) and (2,3)
Step 2: Classify intersections
At (4,1): - Edge (4,0) \to (4,3) goes from outside B to inside B → Entering
At (2,3): - Edge (4,3) \to (0,3) goes from inside B to outside B → Leaving
Step 3: Build intersection polygon
Start at entering intersection (4,1): - Follow A edge to (4,3) - Continue to leaving intersection (2,3) - Switch to B edges - Follow B bottom from (2,1) to entering intersection (4,1)
Intersection polygon: (4,1) \to (4,3) \to (2,3) \to (2,1) \to (4,1)
Result: Rectangle with corners (2,1), (4,1), (4,3), (2,3)
Area: (4-2) \times (3-1) = 2 \times 2 = 4 square units ✓
Below is an interactive polygon overlay demonstration.
<label>
Operation:
<select id="overlay-op">
<option value="intersection" selected>Intersection (A ∩ B)</option>
<option value="union">Union (A ∪ B)</option>
<option value="difference">Difference (A \ B)</option>
<option value="sym-difference">Symmetric Difference (A ⊕ B)</option>
</select>
</label>
<label>
Show inputs:
<input type="checkbox" id="show-inputs" checked>
</label>
<label>
Show intersections:
<input type="checkbox" id="show-intersections" checked>
</label>
<div class="overlay-info">
<p><strong>Polygon A area:</strong> <span id="area-a">--</span></p>
<p><strong>Polygon B area:</strong> <span id="area-b">--</span></p>
<p><strong>Result area:</strong> <span id="area-result">--</span></p>
<p><strong>Intersections:</strong> <span id="intersection-count">--</span></p>
</div><canvas id="overlay-canvas" width="600" height="500" style="border: 1px solid #ddd;"></canvas>Try this: - Intersection: Purple area where polygons overlap (A AND B) - Union: Yellow area covered by either polygon (A OR B) - Difference: Red area in A but not in B (A NOT B) - Show intersections: Red dots where edges cross - Hide inputs: See only the result polygon - Notice: Area(A ∩ B) + Area(A B) + Area(B A) = Area(A ∪ B)
Key insight: Overlay operations are set operations on continuous space—the foundation of “What areas meet multiple criteria?”
Problem: Find suitable locations for solar farm: - Must be: flat land (slope < 5°) - Must be: within 2km of transmission lines - Must NOT be: protected areas or wetlands
Solution using overlay:
1. suitable_slope = reclassify(DEM, slope < 5°)
2. near_power = buffer(transmission_lines, 2000m)
3. buildable = suitable_slope ∩ near_power
4. final = buildable \ (protected ∪ wetlands)Each step is a polygon operation.
Question: How much residential land is in the 100-year floodplain?
residential = parcels WHERE zone = 'R1'
at_risk = residential ∩ floodplain_100yr
area_at_risk = sum(area(at_risk))Result: Quantifies flood exposure for land use planning.
Measure connectivity: Find contiguous forest patches.
forest = land_cover WHERE type = 'forest'
continuous = union(forest) # Merge touching polygons
patches = count_components(continuous)Fragmentation index: More patches = more fragmented = lower wildlife connectivity.
Touching polygons (shared edge, no overlap):
Intersection = line segment (zero area) or empty
Expected: Empty polygon
Reality: May return degenerate polygon (all vertices collinear)
Solution: Check result area > threshold, discard if zero/tiny.
Touching at single vertex:
Similarly produces degenerate result.
Floating-point arithmetic causes issues:
Intersection point calculated as (5.0000000001, 3.0)
Expected: (5.0, 3.0)Result: Duplicate vertices, tiny slivers, self-intersecting polygons
Solution: - Round coordinates to tolerance (e.g., 6 decimal places) - Snap vertices within threshold - Use robust predicates (exact arithmetic for orientation tests)
Very thin polygons from near-parallel edges:
Edge of A at y = 10.00001
Edge of B at y = 10.00000Intersection: Polygon with width 0.00001 units
Problem: Inflates vertex count, causes numerical instability, visually ugly
Solution: Remove slivers below area threshold.
Union of non-overlapping polygons:
Result should be two separate polygons, not one.
Simple implementation: Returns single polygon with disconnected components (invalid).
Proper handling: Return MultiPolygon (collection of polygons).
Real datasets have multiple polygons:
All protected areas = collection of many parks/preserves.
MultiPolygon: Set of disjoint polygons treated as single feature.
Overlay operations:
\text{MultiPoly}_1 \cap \text{MultiPoly}_2 = \bigcup_{i,j} (P_{1i} \cap P_{2j})
Algorithm: 1. For each polygon in MultiPolygon1: 2. For each polygon in MultiPolygon2: 3. Compute pairwise intersection 4. Collect all non-empty results
Union similar: Collect all polygons, merge overlapping ones.
Complexity: O(nm) where n, m are polygon counts.
Optimization: Spatial index (R-tree) to skip non-overlapping pairs.
Intersection (AND):
A \cap B = \{x : x \in A \text{ AND } x \in B\}
Properties: - Commutative: A \cap B = B \cap A - Associative: (A \cap B) \cap C = A \cap (B \cap C) - Identity: A \cap U = A (where U is universal set)
Union (OR):
A \cup B = \{x : x \in A \text{ OR } x \in B\}
Properties: - Commutative: A \cup B = B \cup A - Associative: (A \cup B) \cup C = A \cup (B \cup C) - Identity: A \cup \emptyset = A
Difference (NOT):
A \setminus B = \{x : x \in A \text{ AND } x \notin B\}
NOT commutative: A \setminus B \neq B \setminus A (generally)
Symmetric Difference (XOR):
A \triangle B = (A \setminus B) \cup (B \setminus A)
Equivalent: A \triangle B = (A \cup B) \setminus (A \cap B)
\overline{A \cup B} = \overline{A} \cap \overline{B}
\overline{A \cap B} = \overline{A} \cup \overline{B}
Application: A \setminus B = A \cap \overline{B}
Spatial meaning: Difference is intersection with complement.