116 lines
4 KiB
Python
116 lines
4 KiB
Python
|
||
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
|
||
# Sort triangle vertices so that so that the normal points "outward"
|
||
if np.dot(self.v1, np.cross(self.v2 - self.v1, self.v3 - self.v1)) < 0:
|
||
self.v3, self.v2 = self.v2, self.v3
|
||
# Edges of the origin vertex and their L2 norms
|
||
self.e1 = self.v2 - self.v1
|
||
self.e2 = self.v3 - self.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)
|
||
self.n = self.n / np.linalg.norm(self.n)
|
||
# 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>')
|