From f676af2fde10c7082b5783e4dbc118b2f226439b Mon Sep 17 00:00:00 2001 From: Thomas Kolb Date: Sat, 29 May 2021 15:22:29 +0200 Subject: [PATCH] Initial version: loads GeoJSON and generates SVG map --- .gitmodules | 3 + geo-countries | 1 + qsomap.py | 255 ++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 259 insertions(+) create mode 100644 .gitmodules create mode 160000 geo-countries create mode 100755 qsomap.py diff --git a/.gitmodules b/.gitmodules new file mode 100644 index 0000000..276e3eb --- /dev/null +++ b/.gitmodules @@ -0,0 +1,3 @@ +[submodule "geo-countries"] + path = geo-countries + url = https://github.com/datasets/geo-countries diff --git a/geo-countries b/geo-countries new file mode 160000 index 0000000..cd9e063 --- /dev/null +++ b/geo-countries @@ -0,0 +1 @@ +Subproject commit cd9e0635901eac20294a57ee3b3ce0684d5e3f1a diff --git a/qsomap.py b/qsomap.py new file mode 100755 index 0000000..5a94b14 --- /dev/null +++ b/qsomap.py @@ -0,0 +1,255 @@ +#!/usr/bin/env python3 + +import svgwrite +import numpy as np +import matplotlib.pyplot as pp +from matplotlib.colors import hsv_to_rgb +import json +import random + +REF_LATITUDE = 49.58666 +REF_LONGITUDE = 11.01250 +# REF_LATITUDE = 1 +# REF_LONGITUDE = 0 + + +def map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R=1): + """ Azimuthal equidistant projection. + + This function takes a point to map in latitude/longitude format as well as + a reference point which becomes the "center" of the map. + + It then projects the point into the 2D plane such that the distance from + the center is proportional to the distance on the Great Circle through the + projected and the reference point. The angle represents the azimuthal + direction of the projected point. + + Args: + lat(numpy.array): Latitudes of the point to project. + lon(numpy.array): Longitudes of the point to project. + ref_lat(float): Latitude of the reference point. + ref_lon(float): Longitude of the reference point. + R(float): Radius (scale) of the map. + + Returns: + x(numpy.array): The calculated x coordinates. + y(numpy.array): The calculated y coordinates. + + """ + dlon = lon - ref_lon + + rho_linear_norm = np.arccos(np.sin(ref_lat) * np.sin(lat) + + np.cos(ref_lat) + * np.cos(lat) * np.cos(dlon)) / np.pi + + rho = R * rho_linear_norm + + theta = np.arctan2(np.cos(lat) * np.sin(dlon), + (np.cos(ref_lat) * np.sin(lat) + - np.sin(ref_lat) * np.cos(lat) + * np.cos(dlon))) + + x = rho * np.sin(theta) + y = -rho * np.cos(theta) + + return x, y + + +def random_country_color(): + h = random.random() + s = 0.7 + v = 0.8 + r, g, b = [int(255.99*x) for x in hsv_to_rgb([h, s, v])] + return f"#{r:02x}{g:02x}{b:02x}" + + +random.seed(0) + +""" Test code +test_lat = [np.pi/2, np.pi/4, 0, -np.pi/4, -np.pi/2] +test_lon = np.arange(-np.pi, np.pi, np.pi/4) + +for lat in test_lat: + for lon in test_lon: + x, y = map_azimuthal_equidistant(np.array([lat]), np.array([lon]), np.pi/2, 0) + + print(f"{lat*180/np.pi:6.3f}, {lon*180/np.pi:6.3f} => {x[0]:6.3f}, {y[0]:6.3f}") +""" + +print("Loading Geodata…") + +with open('geo-countries/data/countries.geojson', 'r') as jfile: + geojson = json.load(jfile) + +print("Finding boundaries…") + +# key: 3-letter country identifier +# data: {full_name, numpy.array(coordinates), numpy.array(proj_coordinates)}. +# coordinates is a list of 2xN arrays, where N is the number of points. Row 0 +# contains the longitude, Row 1 the latitude. +# proj_coordinates is a list of 2xN arrays, where N is the number of points. +# Row 0 contains the projected x, Row 1 the projected y. +simplegeodata = {} + +features = geojson['features'] + +for feature in features: + name = feature['properties']['ADMIN'] + key = feature['properties']['ISO_A2'] + + # handle duplicate keys (can happen for small countries) + if key in simplegeodata.keys(): + key = name + + print(f"Preparing {key} ({name})…") + + multipoly = feature['geometry']['coordinates'] + + conv_polys = [] + + for poly in multipoly: + for subpoly in poly: + coords_list = [] # list of lists + + assert(len(subpoly[0]) == 2) + coords_list += subpoly + + # convert coordinates to numpy array and radians + coords = np.array(coords_list).T * np.pi / 180 + + conv_polys.append(coords) + + simplegeodata[key] = {"name": name, "coordinates": conv_polys} + +ref_lat = REF_LATITUDE * np.pi / 180 +ref_lon = REF_LONGITUDE * np.pi / 180 + +R = 500 + +""" +# Override data with test coordinate system +coords = [] + +N = 128 + +# constant-latitude circles +coords.append(np.array([np.linspace(-np.pi, np.pi, N), + np.ones(N) * np.pi/4])) +coords.append(np.array([np.linspace(-np.pi, np.pi, N), + np.ones(N) * 0])) +coords.append(np.array([np.linspace(-np.pi, np.pi, N), + np.ones(N) * -np.pi/4])) + +# constant-longitude half-circles +coords.append(np.array([np.ones(N) * -4*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * -3*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * -2*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * -1*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * 0*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * 1*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * 2*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) +coords.append(np.array([np.ones(N) * 3*np.pi/4, + np.linspace(-np.pi/2, np.pi/2, N)])) + +simplegeodata = {"XY": {'name': 'test', 'coordinates': coords}} +""" + +# apply azimuthal equidistant projection +for k, v in simplegeodata.items(): + proj_polys = [] + + for poly in v['coordinates']: + lat = poly[1, :] + lon = poly[0, :] + + x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) + + coords = np.array([x, y]) + + # remove any points that contain a NaN coordinate + coords = coords[:, np.any(np.invert(np.isnan(coords)), axis=0)] + + proj_polys.append(coords) + + v['proj_coordinates'] = proj_polys + +#print(simplegeodata['AW']['proj_coordinates']) +###print(simplegeodata['DE']['proj_coordinates']) + +#debug plot +# generate the SVG + +doc = svgwrite.Drawing("/tmp/test.svg", size=(2*R, 2*R)) + +doc.defs.add(doc.style(""" + .country { + stroke: black; + stroke-width: 0.01; + } +""")) + +doc.add(doc.circle(center=(R, R), r=R, fill='#ddeeff', + stroke_width=1, stroke='black')) + +for k, v in simplegeodata.items(): + print(f"Exporting {k}…") + + color = random_country_color() + + group = doc.g() + + for poly in v['proj_coordinates']: + points = poly.T + R # shift to the center of the drawing + + pgon = doc.polygon(points, **{ + 'class': 'country', + 'fill': color}) + + group.add(pgon) + + group.set_desc(title=v['name']) + doc.add(group) + +# generate equidistant circles + +d_max = 40075/2 +for distance in [500, 1000, 2000, 3000, 4000, 5000, 6000, 8000, 10000, 12000, + 14000, 16000, 18000, 20000]: + r = R * distance / d_max + doc.add(doc.circle(center=(R, R), r=r, fill='none', + stroke_width=0.1, stroke='black')) + +print(f"Saving {doc.filename}…") +doc.save(pretty=True) + +exit(0) + +# Debug Plot + +for k, v in simplegeodata.items(): + for poly in v['proj_coordinates']: + pp.plot(poly[0, :], poly[1, :]) + +pp.plot([-1, 1], [0, 0], 'k', linewidth=0.5) +pp.plot([0, 0], [-1, 1], 'k', linewidth=0.5) + +t = np.linspace(-np.pi, np.pi, 256) +ct, st = np.cos(t), np.sin(t) +pp.plot(ct, st, 'k', linewidth=0.5) + +U = 40075 +for distance in np.arange(0, U/2, 2000): + f = distance / (U/2) + pp.plot(f*ct, f*st, 'k', linewidth=0.2) + +pp.axis('equal') + +pp.show() +