decePubClient/tools/globe/prepare-countries.py

#!/usr/bin/env python3
# Builds wwwroot/vendor/globe/countries.json, the globe's country borders and name positions, from world-atlas 2.0.2's
# countries-50m.json (Natural Earth 1:50m, public domain; packaged under ISC by Mike Bostock), downloaded once from
#   https://registry.npmjs.org/world-atlas/-/world-atlas-2.0.2.tgz
# usage: tools/globe/prepare-countries.py <countries-50m.json> <iso_3166-1.json> <output json>     (not part of the build)
# iso_3166-1.json is the iso-codes package's (/usr/share/iso-codes/json/iso_3166-1.json); it maps world-atlas's numeric
# ids to the two-letter codes the browser names countries by (Intl.DisplayNames), so names follow the reader's language.
#
# Output: {"borders": [[lng, lat, lng, lat, ...], ...], "countries": [{"code", "name", "lat", "lng", "area"}, ...]}
# - borders are the arcs two countries share, land borders only: coastlines are in the texture already. Coordinates are
#   rounded to 0.01° (about a kilometre).
# - each country's name sits in its largest polygon, at the point nearest the polygon's centroid that is still at least
#   half as far from the edges as the pole of inaccessibility (the point farthest from them, as Mapbox's polylabel finds
#   it): inside the country even for a crescent, and in the middle of a long one (central Italy, not the Po valley).
# - area is the country's, size the square root of its largest polygon's (where the name is), both in degrees scaled by
#   the cosine of the latitude; the globe shows a name once its polygon is large enough on screen to hold it.
import heapq
import json
import math
import sys


def decode(topology):
    scale_x, scale_y = topology["transform"]["scale"]
    translate_x, translate_y = topology["transform"]["translate"]
    arcs = []
    for arc in topology["arcs"]:
        x = y = 0
        points = []
        for dx, dy in arc:
            x += dx
            y += dy
            points.append((x * scale_x + translate_x, y * scale_y + translate_y))
        arcs.append(points)
    return arcs


def ring(arcs, indexes):
    points = []
    for index in indexes:
        arc = arcs[index] if index >= 0 else list(reversed(arcs[~index]))
        points.extend(arc if not points else arc[1:])
    return points


def polygons(geometry):
    if geometry["type"] == "Polygon":
        return [geometry["arcs"]]
    if geometry["type"] == "MultiPolygon":
        return geometry["arcs"]
    return []


def planar(points, latitude):
    # an equirectangular plane squeezed by the cosine of the polygon's latitude, so distances are roughly true there
    k = math.cos(math.radians(latitude))
    return [(x * k, y) for x, y in points]


def area(points):
    return abs(sum(x1 * y2 - x2 * y1 for (x1, y1), (x2, y2) in zip(points, points[1:] + points[:1]))) / 2


def distance_to_rings(x, y, rings):
    inside = False
    best = math.inf
    for points in rings:
        for (ax, ay), (bx, by) in zip(points, points[1:] + points[:1]):
            if (ay > y) != (by > y) and x < (bx - ax) * (y - ay) / (by - ay) + ax:
                inside = not inside
            dx, dy = bx - ax, by - ay
            t = 0 if dx == dy == 0 else max(0, min(1, ((x - ax) * dx + (y - ay) * dy) / (dx * dx + dy * dy)))
            best = min(best, (x - ax - t * dx) ** 2 + (y - ay - t * dy) ** 2)
    return (1 if inside else -1) * math.sqrt(best)


def polylabel(rings, precision=0.05):
    xs = [x for x, _ in rings[0]]
    ys = [y for _, y in rings[0]]
    min_x, min_y, max_x, max_y = min(xs), min(ys), max(xs), max(ys)
    size = min(max_x - min_x, max_y - min_y)
    if size == 0:
        return min_x, min_y
    half = size / 2
    cells = []

    def push(x, y, h):
        d = distance_to_rings(x, y, rings)
        heapq.heappush(cells, (-(d + h * math.sqrt(2)), d, x, y, h))

    x = min_x
    while x < max_x:
        y = min_y
        while y < max_y:
            push(x + half, y + half, half)
            y += size
        x += size
    best_d = distance_to_rings((min_x + max_x) / 2, (min_y + max_y) / 2, rings)
    best = ((min_x + max_x) / 2, (min_y + max_y) / 2)
    while cells:
        bound, d, x, y, h = heapq.heappop(cells)
        if d > best_d:
            best_d, best = d, (x, y)
        if -bound - best_d <= precision:
            continue
        h /= 2
        for sx in (-1, 1):
            for sy in (-1, 1):
                push(x + sx * h, y + sy * h, h)
    return best, best_d


def centroid(points):
    twice = 0
    x = y = 0
    for (x1, y1), (x2, y2) in zip(points, points[1:] + points[:1]):
        cross = x1 * y2 - x2 * y1
        twice += cross
        x += (x1 + x2) * cross
        y += (y1 + y2) * cross
    if twice == 0:
        return points[0]
    return x / (3 * twice), y / (3 * twice)


def label_point(rings):
    (best_x, best_y), best_d = polylabel(rings)
    if best_d <= 0:
        return best_x, best_y
    cx, cy = centroid(rings[0])
    xs = [x for x, _ in rings[0]]
    ys = [y for _, y in rings[0]]
    step = max(max(xs) - min(xs), max(ys) - min(ys)) / 80
    chosen, chosen_distance = (best_x, best_y), math.hypot(best_x - cx, best_y - cy)
    x = min(xs)
    while x <= max(xs):
        y = min(ys)
        while y <= max(ys):
            to_centroid = math.hypot(x - cx, y - cy)
            if to_centroid < chosen_distance and distance_to_rings(x, y, rings) >= best_d / 2:
                chosen, chosen_distance = (x, y), to_centroid
            y += step
        x += step
    return chosen


topology = json.load(open(sys.argv[1]))
codes = {entry["numeric"]: entry["alpha_2"] for entry in json.load(open(sys.argv[2]))["3166-1"]}
arcs = decode(topology)
geometries = topology["objects"]["countries"]["geometries"]

uses = [0] * len(arcs)
for geometry in geometries:
    for polygon in polygons(geometry):
        for indexes in polygon:
            for index in indexes:
                uses[index if index >= 0 else ~index] += 1

borders = [[round(c, 2) for point in arc for c in point] for arc, used in zip(arcs, uses) if used >= 2]

countries = []
for geometry in geometries:
    shapes = []
    for polygon in polygons(geometry):
        outer = ring(arcs, polygon[0])
        latitude = sum(y for _, y in outer) / len(outer)
        shapes.append((area(planar(outer, latitude)), latitude, [ring(arcs, r) for r in polygon]))
    if not shapes:
        continue
    largest, latitude, rings = max(shapes, key=lambda s: s[0])
    k = math.cos(math.radians(latitude))
    x, y = label_point([planar(r, latitude) for r in rings])
    countries.append({
        "code": codes.get(str(geometry.get("id", "")).zfill(3)),
        "name": geometry.get("properties", {}).get("name"),
        "lat": round(y, 2),
        "lng": round(x / k, 2),
        "area": round(sum(s[0] for s in shapes), 1),
        "size": round(math.sqrt(largest), 2)
    })

json.dump({"borders": borders, "countries": sorted(countries, key=lambda c: -c["area"])}, open(sys.argv[3], "w"), separators=(",", ":"))
print(f"{len(borders)} border lines ({sum(len(b) // 2 for b in borders)} points), {len(countries)} countries "
      f"({sum(1 for c in countries if c['code'] is None)} without an ISO code)")