103 lines
3.4 KiB
Python
103 lines
3.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
|
|||
|
|
# 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, v - self.v1) / self.e1n
|
|||
|
|
n = np.cross(self.n, self.e1 / self.e1n)
|
|||
|
|
v = np.dot(n, v-self.v2)
|
|||
|
|
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 = []
|
|||
|
|
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:
|
|||
|
|
mapped.append((tri.plane_coords(p0i), tri.plane_coords(p1i)))
|
|||
|
|
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('/tmp/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'][0]:
|
|||
|
|
countries[cc].extend(map_poly(t, poly))
|
|||
|
|
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 a, b in l:
|
|||
|
|
print(f'<path fill="none" stroke="red" stroke-width="0.01" d="M {a[0]} {a[1]} L {b[0]} {b[1]}" />')
|
|||
|
|
print('</svg>')
|