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.

Effective area: the quantity Visvalingam–Whyatt actually ranks A coastline polyline with seven vertices. Around the third vertex, the triangle formed with its two neighbours is shaded and labelled with an effective area of 210 square metres. Around the fifth vertex, a much flatter triangle is labelled 18 square metres. A note explains that the flattest triangle is removed first, its two neighbours have their areas recomputed, and the heap is updated, so removal order is by shape contribution rather than by distance from a chord. A = 210 m² a genuine headland — kept A = 18 m² survey noise — removed first effective area = area of the triangle with the two current neighbours After a removal the two surviving neighbours get new triangles and new areas, so the heap must be updated rather than rebuilt — that update is the whole cost difference against a naive re-sort.
Ranking by area rather than by perpendicular distance is why this algorithm degrades gracefully: at any tolerance the vertices that survive are the ones carrying shape, not the ones that happen to be far from a chord.

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.
Retention and fidelity as the threshold rises Two curves of retained vertices against threshold for the same coastline. Visvalingam–Whyatt retains fewer vertices at every threshold while holding a higher shape-fidelity score, plotted as a dashed line. Douglas–Peucker keeps more vertices but concentrates them along straight runs, so its fidelity score falls faster once the threshold passes the survey-noise band. 0% 50% 100% threshold, increasing → Visvalingam–Whyatt retention Douglas–Peucker retention shape fidelity, VW shape fidelity, DP The useful reading is the vertical gap between each retention curve and its fidelity curve: it is how much recognisable coastline each algorithm buys per vertex it keeps, and it is where the two clearly part.
Both algorithms remove vertices; they disagree about which ones are expendable. On coastlines that difference shows up as whether bays survive the same file-size budget.

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