halkyon/apps/learn/tool/maps/outlines.py
Tizian 「ティツィアーン」 649499be23
chore: Halkyon Learn moves into apps/learn
Its server, client, app and curator, and the tools that write its
assets (maps, Anki fixtures, icons), as they are; paths are fixed in
the commits that follow.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Signed-off-by: Tizian 「ティツィアーン」 <raisondetredev@pm.me>
2026-10-22 18:11:02 +02:00

102 lines
3.7 KiB
Python

"""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}