globetiles/map.py

from math import *
import numpy as np
import sys
import json


EPS = sys.float_info.epsilon


class Triangle:

    def __init__(self, v1, v2, v3):
        # Vertices
        self.v1 = v1
        self.v2 = v2
        self.v3 = v3
        # Edges of the origin vertex and their L2 norms
        self.e1 = v2 - v1
        self.e2 = v3 - v1
        self.e1n = np.linalg.norm(self.e1)
        self.e2n = np.linalg.norm(self.e2)
        # Normal vector of the triangle's plane
        self.n = np.cross(self.e1 / self.e1n, self.e2 / self.e2n)
        # The vertex coordinates in plane coordinates
        self.v1p = self.plane_coords(self.v1)
        self.v2p = self.plane_coords(self.v2)
        self.v3p = self.plane_coords(self.v3)

    def intersect(self, ray):
        # https://en.wikipedia.org/wiki/Möller–Trumbore_intersection_algorithm
        ray_cross_e2 = np.cross(ray, self.e2)
        dot = np.dot(self.e1, ray_cross_e2)
        if abs(dot) < EPS:
            return None
        s = -self.v1
        u = np.dot(s, ray_cross_e2) / dot
        if (u < 0 and abs(u) > EPS) or (u > 1 and abs(u - 1) > EPS):
            return None
        s_cross_e1 = np.cross(s, self.e1)
        v = np.dot(ray, s_cross_e1) / dot
        if (v < 0 and abs(v) > EPS) or (u + v > 1 and abs(u + v - 1) > EPS):
            return None
        t = np.dot(self.e2, s_cross_e1) / dot
        if t <= EPS:
            return None
        return ray * t

    def plane_coords(self, v):
        u = np.dot(self.e1 / self.e1n, v - self.v1)
        n = np.cross(self.n, self.e1 / self.e1n)
        v = -np.dot(n, v - self.v1)
        return np.array([u, v])


def convert_to_cartesian(lon, lat):
    x = cos(lon*pi/180)
    y = sin(lon*pi/180)
    z = tan(lat*pi/180)
    mag = sqrt(x**2 + y**2 + z**2)
    return np.array([x/mag, y/mag, z/mag])


def map_poly(tri, poly):
    p0 = convert_to_cartesian(*poly[0])
    p0i = tri.intersect(p0)
    mapped = []
    new = True
    for i in range(1, len(poly)):
        p1 = convert_to_cartesian(*poly[i])
        p1i = tri.intersect(p1)
        if p0i is not None and p1i is not None:
            if new:
                mapped.append([tri.plane_coords(p0i), tri.plane_coords(p1i)])
                new = False
            else:
                mapped[-1].append(tri.plane_coords(p1i))
        else:
            new = True
        p0 = p1
        p0i = p1i
    return mapped


if __name__ == '__main__':
    t = Triangle(np.array([5,0,0]), np.array([0,5,0]), np.array([0,0,5]))
    countries = {}
    with open('countries.geojson', 'r') as f:
        j = json.load(f)
        for f in j['features']:
            cc = f['properties']['ADMIN']
            if f['geometry']['type'] == 'MultiPolygon':
                countries[cc] = []
                for poly in f['geometry']['coordinates']:
                    countries[cc].extend(map_poly(t, poly[0]))
            elif f['geometry']['type'] == 'Polygon':
                countries[cc] = map_poly(tri, f['geometry']['coordinates'][0])
    minx = min(t.v1p[0], t.v2p[0], t.v3p[0])
    miny = min(t.v1p[1], t.v2p[1], t.v3p[1])
    maxx = max(t.v1p[0], t.v2p[0], t.v3p[0])
    maxy = max(t.v1p[1], t.v2p[1], t.v3p[1])
    w = maxx - minx
    h = maxy - miny
    print(f'<svg xmlns="http://www.w3.org/2000/svg" viewBox="{minx} {miny} {w} {h}" width="{w}" height="{h}">')
    print(f'<path fill="none" stroke="black" stroke-width="0.01" d="M {t.v1p[0]} {t.v1p[1]} L {t.v2p[0]} {t.v2p[1]} L {t.v3p[0]} {t.v3p[1]} Z" />')
    for c, l in countries.items():
        for ls in l:
            print(f'<path fill="none" stroke="red" stroke-width="0.01" d="M {ls[0][0]} {ls[0][1]}', end='')
            for x, y in ls[1:]:
                print(f' L {x} {y}', end='')
            print('" />')
    print('</svg>')