184 lines
7.0 KiB
Python
184 lines
7.0 KiB
Python
 
|
#!/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)")
|