Simplifying Coastlines with Visvalingam-Whyatt in Python

Visvalingam-Whyatt simplifies a coastline by ranking every vertex by the area of the triangle it forms with its two immediate neighbours and removing the smallest-area vertices first; because a vertex’s importance is measured by the area it encloses rather than its distance from a chord, the algorithm preserves a coastline’s visual character — bays, spits, and headlands — far better than Douglas-Peucker at aggressive vertex reduction. A heap-based implementation runs in O(n log n) and lets you stop on either an effective-area threshold or a target vertex budget.

Core Algorithm and Workflow

The insight behind Visvalingam-Whyatt is that a vertex matters in proportion to how much shape it contributes. Take any interior vertex B with neighbours A and C: the triangle ABC has an effective area equal to 0.5 * abs((Ax - Cx) * (By - Ay) - (Ax - Bx) * (Cy - Ay)). A vertex sitting almost on the straight line A → C encloses a near-zero triangle and can be dropped with negligible visual loss; a vertex at the tip of a headland encloses a large triangle and must survive. The algorithm removes vertices in ascending order of effective area, recomputing the two affected neighbours after each removal.

The diagram below shows the effective-area triangle for a single vertex and the order in which the low-area vertices are removed.

Visvalingam-Whyatt effective area and removal order A jagged coastline polyline. One interior vertex B is joined to its neighbours A and C by a shaded triangle showing its effective area. Three low-area vertices along the line are numbered 1, 2 and 3 to show they are removed first, smallest area first. A B (keep: large area) C effective area = area(ABC) remove 2 remove 1 remove 3 Hollow vertices enclose small triangles and are dropped first; solid B survives.

Numbered vertices remove 1, remove 2, and remove 3 sit nearly on the line through their neighbours, so their triangles are tiny and they leave first. Vertex B — the headland tip — encloses a large triangle and stays until the very end. This is exactly the behaviour that separates the method from Douglas-Peucker, whose comparison against a chord tolerance would discard an entire shallow bay in one pass. The parent overview of generalization and simplification covers where each algorithm is the right default.

Production-Ready Python Implementation

The implementation below uses a min-heap keyed on effective area and a doubly linked ring so that removing a vertex and re-linking its neighbours is O(1), with only the two touched neighbours re-pushed onto the heap. Lazy deletion — a removed flag plus a per-vertex version stamp — avoids the cost of an arbitrary heap delete. It targets shapely 2.0.x and numpy 1.26.x; see the Shapely geometry docs for the constructors used here.

import heapq
from shapely.geometry import Polygon, MultiPolygon, LinearRing


def _triangle_area(a, b, c) -> float:
    """Twice-signed area halved: the effective area of vertex b given neighbours a, c."""
    return 0.5 * abs(
        (a[0] - c[0]) * (b[1] - a[1]) - (a[0] - b[0]) * (c[1] - a[1])
    )


def simplify_ring_vw(coords, area_threshold=None, target_vertices=None):
    """
    Visvalingam-Whyatt simplification of a single ring or line.

    coords          : list of (x, y) tuples in a PROJECTED CRS (metres).
    area_threshold  : stop once the smallest effective area exceeds this
                      value (in CRS units squared). Coarser = larger value.
    target_vertices : alternatively, stop when this many vertices remain.

    Exactly one of area_threshold / target_vertices should be supplied.
    Returns the surviving coordinates in original order.
    """
    is_closed = coords[0] == coords[-1]
    pts = coords[:-1] if is_closed else list(coords)
    n = len(pts)
    if n <= 3:
        return list(coords)

    # Doubly linked ring via prev/next index arrays.
    prev = [(i - 1) % n for i in range(n)]
    nxt = [(i + 1) % n for i in range(n)]
    alive = [True] * n
    version = [0] * n

    def eff_area(i) -> float:
        # Endpoints of an open line are pinned so they are never removed.
        if not is_closed and (i == 0 or i == n - 1):
            return float("inf")
        return _triangle_area(pts[prev[i]], pts[i], pts[nxt[i]])

    heap = [(eff_area(i), i, version[i]) for i in range(n)]
    heapq.heapify(heap)

    # Minimum survivors: 3 for a closed ring, 2 for an open line.
    floor_count = 3 if is_closed else 2
    remaining = n

    while heap and remaining > floor_count:
        area, i, ver = heapq.heappop(heap)
        if not alive[i] or ver != version[i]:
            continue  # stale entry from lazy deletion
        if target_vertices is not None and remaining <= target_vertices:
            break
        if area_threshold is not None and area > area_threshold:
            break

        # Splice vertex i out and re-link its neighbours.
        p, q = prev[i], nxt[i]
        nxt[p], prev[q] = q, p
        alive[i] = False
        remaining -= 1

        # Only the two neighbours change; re-push them with a fresh version.
        for nb in (p, q):
            if alive[nb]:
                version[nb] += 1
                heapq.heappush(heap, (eff_area(nb), nb, version[nb]))

    # Walk the surviving ring in order.
    start = next(i for i in range(n) if alive[i])
    out, cur = [], start
    while True:
        out.append(pts[cur])
        cur = nxt[cur]
        if cur == start:
            break
    if is_closed:
        out.append(out[0])
    return out


def simplify_polygon_vw(poly: Polygon, **kw) -> Polygon | None:
    """Simplify the exterior and every hole (island lagoon) independently."""
    ext = simplify_ring_vw(list(poly.exterior.coords), **kw)
    if len(ext) < 4:                       # ring collapsed below a triangle
        return None
    holes = []
    for ring in poly.interiors:
        h = simplify_ring_vw(list(ring.coords), **kw)
        if len(h) >= 4:                    # drop holes that vanish
            holes.append(h)
    return Polygon(LinearRing(ext), holes)


def simplify_coastline(geom, min_feature_area=0.0, **kw):
    """
    Simplify a coastline Polygon or MultiPolygon (islands + mainland).
    Parts smaller than min_feature_area (CRS units squared) are dropped
    entirely rather than reduced to a degenerate sliver.
    """
    parts = geom.geoms if isinstance(geom, MultiPolygon) else [geom]
    kept = []
    for part in parts:
        if part.area < min_feature_area:
            continue
        simplified = simplify_polygon_vw(part, **kw)
        if simplified is not None and simplified.is_valid:
            kept.append(simplified)
    if not kept:
        return None
    return kept[0] if len(kept) == 1 else MultiPolygon(kept)


# ---------------------------------------------------------------------------
# Example: simplify a projected coastline to roughly a quarter of its detail
# ---------------------------------------------------------------------------
# import geopandas as gpd
# gdf = gpd.read_file("coastline.gpkg").to_crs("EPSG:3035")   # metres
# # area_threshold in m^2: features whose effective area is smaller vanish.
# gdf["geometry"] = gdf.geometry.apply(
#     lambda g: simplify_coastline(g, min_feature_area=50_000.0,
#                                  area_threshold=25_000.0)
# )
# gdf = gdf[gdf.geometry.notna()]
# gdf.to_file("coastline_simplified.gpkg", driver="GPKG")

The two stop conditions serve different needs. area_threshold is scale-driven: a vertex whose triangle is smaller than the area of one map pixel at the target scale cannot be seen, so it should go. target_vertices is budget-driven: when a downstream renderer or tile format has a hard vertex ceiling, cap the survivor count directly. Deriving a sensible area_threshold from a display scale is covered under scale mapping for web and print.

Performance Tuning and Cartographic Best Practices

  • Always project to metres before ranking areas. Effective area is computed in CRS units squared. Run the resolver on geographic coordinates (EPSG:4326) and a degree² threshold is meaningless and latitude-biased. Reproject to an equal-area system such as EPSG:3035 for Europe or a local UTM zone first, exactly as you would when auto-detecting UTM zones from bounding boxes for a regional dataset.
  • Prefer the area threshold for multi-feature layers. A single area_threshold gives a consistent minimum feature size across every island and inlet in the layer, whereas a per-feature target_vertices makes a large mainland and a tiny islet equally coarse. Reserve the vertex budget for single-feature exports bound by a format limit.
  • Cache effective areas as a weight, not a decision. The classic Visvalingam-Whyatt refinement is to store, per vertex, the maximum area at which it was removed and enforce monotonicity (area = max(area, last_removed_area)). Persisting that weight lets you re-derive any simplification level on the fly — ideal for zoom-dependent tiling without re-running the heap.
  • Guard against self-intersection on concave coasts. Area-ranked removal can fold a deep inlet across itself. Validate each output with geom.is_valid; on failure, re-run that ring with a smaller threshold or repair with geom.buffer(0) before it reaches the renderer.
  • Simplify shared borders once via topology. A coastline segment shared with an administrative polygon must be simplified identically on both sides or a visible gap appears. Feed such layers through TopoJSON so shared arcs are reduced a single time; independent per-polygon runs will diverge.

Integration and Next Steps

The simplify_coastline function returns standard Shapely geometries, so it drops into any GeoPandas pipeline and serializes to GeoPackage, GeoJSON, or PostGIS unchanged. Two integrations are common:

  • Vector tiles. Simplify once to a generous area_threshold, then hand the reduced GeoPackage to a tiler for per-zoom generalization. This keeps the tiler’s own simplification honest because the input is already clean — the workflow for generating MBTiles from GeoPackage with Tippecanoe pairs directly with this output.
  • Print export. For a fixed-scale plate, pick area_threshold from the plate’s DPI and scale so the smallest surviving triangle is about one device pixel, then export the simplified geometry as vectors through the high-resolution vector export workflow to keep file size and node count low without visible loss.

Store the per-vertex removal weight described above alongside your source geometry and a single simplification pass can serve every zoom level, turning coastline generalization from a batch step into a lookup.

Back to Generalization and Simplification