from math import * import numpy as np import sys import json from tqdm import tqdm, trange 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): ''' :param ray: numpy.array(3) Ray cast from [0, 0, 0] :return: A tuple of [bool, Optional[numpy.array(3)]]: The first component is True if the ray intersects the triangle, False otherwise. The second component marks the intersection of the triangle's plane (even if outside the triangle), or None if the ray does not intersect the plane at all. ''' # 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 False, None s = -self.v1 u = np.dot(s, ray_cross_e2) / dot s_cross_e1 = np.cross(s, self.e1) v = np.dot(ray, s_cross_e1) / dot t = np.dot(self.e2, s_cross_e1) / dot if (u < 0 and abs(u) > EPS) or (u > 1 and abs(u - 1) > EPS): return False, ray * t if (v < 0 and abs(v) > EPS) or (u + v > 1 and abs(u + v - 1) > EPS): return False, ray * t if t <= EPS: return False, None return True, 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]) class Face: def __init__(self, vs): self.vertices = vs self.tris = [] for i in range(1, len(vs)-1): self.tris.append(Triangle(vs[0], vs[i], vs[i+1])) self.plane_coords = [self.tris[0].plane_coords(v) for v in self.vertices] 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, ref=None): if ref is None: ref = tri p0 = convert_to_cartesian(*poly[0]) p0b, p0i = tri.intersect(p0) mapped = [] new = True for i in range(1, len(poly)): p1 = convert_to_cartesian(*poly[i]) p1b, p1i = tri.intersect(p1) if p0b or p1b and p0i is not None and p1i is not None: if new: mapped.append([ref.plane_coords(p0i), ref.plane_coords(p1i)]) new = False else: mapped[-1].append(ref.plane_coords(p1i)) else: new = True p0 = p1 p0i = p1i p0b = p1b return mapped if __name__ == '__main__': fs = [Face([a, b, c, d]) for a, b, c, d in [ (np.array([1,1,1]), np.array([-1,1,1]), np.array([-1,-1,1]), np.array([1,-1,1])), # top (np.array([1,1,-1]), np.array([-1,1,-1]), np.array([-1,-1,-1]), np.array([1,-1,-1])), # bottom (np.array([-1,-1,-1]), np.array([-1,1,-1]), np.array([-1,1,1]), np.array([-1,-1,1])), # left (np.array([1,-1,-1]), np.array([1,1,-1]), np.array([1,1,1]), np.array([1,-1,1])), # right (np.array([-1,-1,-1]), np.array([1,-1,-1]), np.array([1,-1,1]), np.array([-1,-1,1])), # front (np.array([-1,1,-1]), np.array([1,1,-1]), np.array([1,1,1]), np.array([-1,1,1])), # back ]] with open('countries.geojson', 'r') as f: j = json.load(f) for i, face in tqdm(list(enumerate(fs)), desc='Faces '): countries = {} for t in face.tris: for f in tqdm(j['features'], desc='Features', total=len(face.tris)*len(j['features'])): cc = f['properties']['ADMIN'] if f['geometry']['type'] == 'MultiPolygon': for poly in f['geometry']['coordinates']: countries.setdefault(cc, []).extend(map_poly(t, poly[0], face.tris[0])) elif f['geometry']['type'] == 'Polygon': countries.setdefault(cc, []).extend(map_poly(t, f['geometry']['coordinates'][0], face.tris[0])) minx = min(v[0] for v in face.plane_coords) miny = min(v[1] for v in face.plane_coords) maxx = max(v[0] for v in face.plane_coords) maxy = max(v[1] for v in face.plane_coords) w = maxx - minx h = maxy - miny with open(f'face{i}.svg', 'w') as svg: svg.write(f'\n') svg.write(f' \n \n \n \n') for c, l in countries.items(): if l: svg.write(f' \n') for ls in l: svg.write(f' \n') svg.write(' \n') svg.write(' \n\n')