halkyon/tool/maps/outlines.py

"""Outlines reduced for display, shared by the map scripts beside this file.

Douglas-Peucker simplification, islands too small to see dropped (a region
always keeps its largest), coordinates rounded. The result is an
approximation, not survey-accurate data.
"""

import hashlib
import math
import sys
import urllib.request


def download(url, sha256):
    """The bytes at url, checked against their recorded SHA-256."""
    print(f"Downloading {url}", file=sys.stderr)
    data = urllib.request.urlopen(url).read()
    digest = hashlib.sha256(data).hexdigest()
    if digest != sha256:
        sys.exit(f"The SHA-256 of {url} is {digest}, not {sha256}.")
    return data


def simplify(points, tolerance):
    """Douglas-Peucker on an open polyline of (x, y) pairs."""
    if len(points) < 3:
        return points
    first, last = points[0], points[-1]
    dx, dy = last[0] - first[0], last[1] - first[1]
    length = math.hypot(dx, dy)
    farthest, index = -1.0, 0
    for i in range(1, len(points) - 1):
        px, py = points[i]
        if length == 0:
            distance = math.hypot(px - first[0], py - first[1])
        else:
            distance = abs(dy * px - dx * py + last[0] * first[1]
                           - last[1] * first[0]) / length
        if distance > farthest:
            farthest, index = distance, i
    if farthest <= tolerance:
        return [first, last]
    left = simplify(points[: index + 1], tolerance)
    right = simplify(points[index:], tolerance)
    return left[:-1] + right


def simplify_ring(ring, tolerance, decimals):
    """A closed ring simplified, rounded, without repeated points; None when
    it collapses below a triangle."""
    # Split at the point farthest from the start, so a closed ring (first
    # point = last point) simplifies as two open halves.
    start = ring[0]
    far = max(range(len(ring)), key=lambda i: math.hypot(
        ring[i][0] - start[0], ring[i][1] - start[1]))
    halves = simplify(ring[: far + 1], tolerance)[:-1] + simplify(
        ring[far:], tolerance)
    rounded = []
    for x, y in halves:
        point = [round(x, decimals), round(y, decimals)]
        if not rounded or rounded[-1] != point:
            rounded.append(point)
    if rounded[0] != rounded[-1]:
        rounded.append(rounded[0])
    return rounded if len(rounded) >= 4 else None


def area(ring):
    """The absolute area of a ring, in square degrees (shoelace)."""
    return abs(sum(x1 * y2 - x2 * y1
                   for (x1, y1), (x2, y2) in zip(ring, ring[1:]))) / 2


def polygons_of(geometry):
    if geometry["type"] == "Polygon":
        return [geometry["coordinates"]]
    if geometry["type"] == "MultiPolygon":
        return geometry["coordinates"]
    raise ValueError(f"Unexpected geometry {geometry['type']}")


def reduce(polygons, tolerance, min_area, decimals):
    """polygons simplified: those smaller than min_area dropped, except the
    largest; a polygon that collapses entirely kept as its outline alone,
    unsimplified but rounded, so even the smallest island stays."""
    largest = max(range(len(polygons)), key=lambda i: area(polygons[i][0]))
    kept = []
    for i, polygon in enumerate(polygons):
        if i != largest and area(polygon[0]) < min_area:
            continue
        exterior = simplify_ring(polygon[0], tolerance, decimals)
        if exterior is None:
            if i != largest:
                continue
            exterior = simplify_ring(polygon[0], 0, decimals)
            if exterior is None:
                continue
        holes = [hole for hole in (simplify_ring(ring, tolerance, decimals)
                                   for ring in polygon[1:])
                 if hole is not None]
        kept.append([exterior, *holes])
    return {"type": "MultiPolygon", "coordinates": kept}