Creating proximity zones around geographic features
2026-02-27
Alberta’s Water Act specifies a 30-metre setback from any permanent or intermittent watercourse for most development and land disturbance activities. BC’s Riparian Areas Regulation sets a 30-metre Streamside Protection and Enhancement Area around fish-bearing streams. These legal instruments work, at the level of individual permits, as polygon intersection problems: does the proposed development boundary overlap the 30-metre buffer polygon of any regulated watercourse? The buffer is not metaphorical — it is a precisely computed geometric object, and its computational properties determine whether the regulatory test is tractable for the thousands of permit applications processed annually.
Buffer operations generate proximity zones: all locations within a specified distance of a geographic feature. The mathematics underlying them is Euclidean distance — the Pythagorean theorem applied cell by cell across a grid, or computed analytically as an offset to a polygon boundary. For point features, the buffer is a circle. For line features, it is a stadium shape (rectangle with semicircular ends). For polygon features, it is an inward or outward offset of the boundary. The computational challenge is that real boundary lines are not straight — they are sequences of vertices connected by line segments, and the offset of a polyline must handle inside corners (which close) and outside corners (which either round or miter). This model derives the distance geometry, builds both raster distance fields and vector offset operations, and shows how they are used in environmental setback analysis and service area mapping.
How far is this location from the nearest protected stream?
Buffer zones appear everywhere in spatial analysis: - Environmental: 100m riparian buffer along streams (no development) - Public safety: 1km evacuation zone around chemical plants - Urban planning: 400m walkable distance to transit stops - Epidemiology: 2km contact tracing radius from infection case - Retail: 5km delivery service area around store
The mathematical question: Given a geographic feature (point, line, or polygon), create a new polygon representing all locations within distance d.
Variations: - Point buffer → circle - Line buffer → “stadium” shape (rounded rectangle) - Polygon buffer → offset boundary (larger or smaller)
Distance from point (x_1, y_1) to point (x_2, y_2):
d = \sqrt{(x_2 - x_1)^2 + (y_2 - y_1)^2}
This is the Pythagorean theorem in coordinate form.
Buffer definition: All points within distance d from a feature.
\text{Buffer} = \{(x, y) : \text{distance}(x, y, \text{feature}) \leq d\}
1. Point buffer:
Distance from point to center equals d → circle of radius d
x^2 + y^2 = d^2
2. Line buffer:
All points within distance d from any point on the line.
For horizontal line segment from (x_1, y) to (x_2, y): - Points above/below: parallel lines at y \pm d - Points at ends: semicircles of radius d - Result: “stadium” or “racetrack” shape
3. Polygon buffer:
Outward buffer (positive): Expand boundary by distance d
Inward buffer (negative): Shrink boundary by distance d
Each edge moves perpendicular to itself, circular arcs at vertices.
Center at (c_x, c_y), radius r:
\{(x, y) : (x - c_x)^2 + (y - c_y)^2 \leq r^2\}
Parametric form (for drawing):
x = c_x + r\cos\theta y = c_y + r\sin\theta
Where \theta \in [0, 2\pi].
Discretization (for polygon representation):
for i = 0 to n-1:
theta = 2π * i / n
x = cx + r * cos(theta)
y = cy + r * sin(theta)
add (x, y) to buffer polygonTypically n = 32 or $64$ points for smooth circle.
For segment from A = (x_1, y_1) to B = (x_2, y_2) with buffer distance d:
Step 1: Calculate perpendicular direction
Direction vector: \mathbf{v} = (x_2 - x_1, y_2 - y_1)
Length: L = \sqrt{(x_2-x_1)^2 + (y_2-y_1)^2}
Unit vector: \hat{v} = \mathbf{v}/L = \left(\frac{x_2-x_1}{L}, \frac{y_2-y_1}{L}\right)
Perpendicular (rotate 90°): \hat{n} = \left(-\frac{y_2-y_1}{L}, \frac{x_2-x_1}{L}\right)
Step 2: Offset parallel edges
Left edge: - A' = A + d\hat{n} - B' = B + d\hat{n}
Right edge: - A'' = A - d\hat{n} - B'' = B - d\hat{n}
Step 3: Add semicircular caps
At A: semicircle from A'' to A' (radius d, center A)
At B: semicircle from B' to B'' (radius d, center B)
Buffer polygon: A' \to B' \to [arc] \to B'' \to A'' \to [arc] \to A'
Minkowski sum approach:
Buffer of polygon P by distance d = Minkowski sum of P and disk of radius d.
Intuitive: Slide a disk of radius d around the polygon boundary.
Implementation strategy:
For each edge: 1. Offset edge outward by d (perpendicular direction) 2. At each vertex, add circular arc connecting adjacent offset edges 3. Handle intersections (offset edges may cross)
Challenges: - Convex vertices: Circular arc expands outward - Concave vertices: May cause self-intersection → requires clipping - Negative buffers: May eliminate small features entirely
Robust implementation: Use computational geometry libraries (GEOS, Clipper).
Problem: Create a 50-meter buffer around a line segment from A = (0, 0) to B = (300, 400) (meters).
Step 1: Calculate length
L = \sqrt{300^2 + 400^2} = \sqrt{90000 + 160000} = \sqrt{250000} = 500 \text{ m}
Step 2: Unit direction vector
\hat{v} = \left(\frac{300}{500}, \frac{400}{500}\right) = (0.6, 0.8)
Step 3: Perpendicular (rotate 90° CCW)
\hat{n} = (-0.8, 0.6)
Step 4: Offset edges by d = 50 m
Left edge (outward):
A' = (0, 0) + 50(-0.8, 0.6) = (-40, 30) B' = (300, 400) + 50(-0.8, 0.6) = (260, 430)
Right edge (inward):
A'' = (0, 0) - 50(-0.8, 0.6) = (40, -30) B'' = (300, 400) - 50(-0.8, 0.6) = (340, 370)
Step 5: Circular caps
At A: Semicircle from (40, -30) to (-40, 30), radius 50, center (0, 0)
At B: Semicircle from (260, 430) to (340, 370), radius 50, center (300, 400)
Buffer polygon vertices (approximate with 8 points per semicircle):
A'(-40, 30) \to B'(260, 430) \to [8 arc points] \to B''(340, 370) \to A''(40, -30) \to [8 arc points] \to A'
Total: ~20 vertices defining the buffer polygon.
Below is an interactive buffer generator.
<label>
Feature type:
<select id="feature-type">
<option value="point">Point</option>
<option value="line" selected>Line</option>
<option value="polygon">Polygon</option>
</select>
</label>
<label>
Buffer distance (pixels):
<input type="range" id="buffer-distance" min="10" max="100" step="5" value="50">
<span id="buffer-value">50</span>
</label>
<label>
Show construction:
<input type="checkbox" id="show-construction" checked>
</label>
<label>
Buffer style:
<select id="buffer-style">
<option value="filled" selected>Filled</option>
<option value="outline">Outline only</option>
</select>
</label><canvas id="buffer-canvas" width="600" height="500" style="border: 1px solid #ddd;"></canvas><p><strong>Buffer area:</strong> <span id="buffer-area">--</span> px²</p>
<p><strong>Buffer perimeter:</strong> <span id="buffer-perimeter">--</span> px</p>Try this: - Point buffer: Perfect circle around point - Line buffer: “Stadium” shape with semicircular caps - Polygon buffer: Offset boundary with rounded corners - Show construction: See the parallel edges and perpendicular offsets - Adjust distance: Watch buffer grow/shrink - Notice: Line buffer = rectangle + two semicircles (area formula confirms)
Key insight: Buffer is the Minkowski sum of the feature and a disk—equivalent to sweeping a circle along the boundary.
Distance field = raster where each cell contains distance to nearest feature.
Algorithm (Euclidean Distance Transform):
for each cell (x, y):
min_dist = infinity
for each feature point (fx, fy):
dist = sqrt((x - fx)² + (y - fy)²)
if dist < min_dist:
min_dist = dist
distance_field[x, y] = min_distOptimization: Use sweep algorithms (Danielsson, Chamfer) for O(n) instead of O(n \times m).
Application: Iso-distance contours from distance field = buffer boundaries at different distances.
Union of buffers:
Buffer around rivers = union of all individual stream buffers.
Algorithm: 1. Buffer each feature independently 2. Compute polygon union (overlay operation) 3. Result: Single merged buffer polygon
Example: 100m buffers around all streams in watershed → riparian protection zone.
Erosion: Shrink polygon by distance d.
Use cases: - Create “core areas” (remove edge effects) - Model sea-level rise (coastline moves inland) - Setbacks (building must be 10m inside property line)
Challenge: Polygon may disappear entirely if d exceeds half the minimum width.
Example: 50m buffer on 80m-wide forest strip → 80m - 2(50m) = -20m → vanishes!
Multiple features close together:
Buffers overlap → merged area may not be simple sum of individual buffers.
Example: Two points 60m apart, each with 50m buffer: - Individual areas: $2 (50^2) = 15,708$ m² - Merged area: < 15,708 m² (overlap region counted once)
Solution: Use union operation, calculate area from merged polygon.
Concave polygon + large buffer:
Offset edges may cross each other.
Example: U-shaped polygon with buffer larger than gap width → buffer fills the gap.
Handling: 1. Detect self-intersections 2. Remove invalid portions (keep outermost boundary) 3. Use robust geometry libraries (GEOS, Shapely)
Very small buffers (< 1 pixel):
May produce degenerate polygons or disappear entirely in raster representation.
Very large buffers:
May exceed coordinate system bounds or cause numerical overflow.
Solution: Validate buffer distance is reasonable for coordinate system and scale.
Buffer may change polygon topology:
Original: Single polygon
After negative buffer: May split into multiple polygons (if narrow connections < 2d)
Example: Dumbbell-shaped polygon with narrow neck. Inward buffer may sever the neck → two separate polygons.
Problem: Euclidean distance is inaccurate over large areas (Earth is curved).
Geodetic buffer: Uses great-circle distance on ellipsoid (WGS84).
Difference:
At equator, 1° longitude ≈ 111 km
At 60°N, 1° longitude ≈ 55 km
Euclidean buffer at 60°N: Would be elliptical (stretched E-W)
Geodetic buffer: Remains circular in ground distance
Implementation: Project to equal-area projection (e.g., Albers), buffer, project back.
When to use: - Features spanning > 1° latitude/longitude - Polar regions (extreme distortion) - Global or continental analysis
When Euclidean is fine: - Local analysis (< 100 km²) - Projected coordinates (UTM, State Plane)
To rotate vector (x, y) by 90° counter-clockwise:
(x, y) \to (-y, x)
Why: Rotation matrix for 90°:
\begin{pmatrix} \cos 90° & -\sin 90° \\ \sin 90° & \cos 90° \end{pmatrix} = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}
\begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} -y \\ x \end{pmatrix}
Given segment from A to B:
Direction vector: \mathbf{v} = B - A
Two perpendiculars (90° rotations): - Left: \mathbf{n}_L = (-v_y, v_x) - Right: \mathbf{n}_R = (v_y, -v_x)
Unit perpendicular (length 1):
\hat{n} = \frac{1}{\|\mathbf{v}\|}(-v_y, v_x)
Offset point by distance d:
P_{\text{offset}} = P + d\hat{n}
This moves point P perpendicular to line direction by distance d.