From a766274d79278b04b99ddfe4a5a4e416ea99a3d5 Mon Sep 17 00:00:00 2001 From: Thomas Kolb Date: Sun, 30 May 2021 22:20:40 +0200 Subject: [PATCH] Support for command line arguments --- qsomap.py | 429 ++++++++++++++++++++++++++++-------------------------- 1 file changed, 224 insertions(+), 205 deletions(-) diff --git a/qsomap.py b/qsomap.py index de05838..7cc4104 100755 --- a/qsomap.py +++ b/qsomap.py @@ -1,16 +1,18 @@ #!/usr/bin/env python3 +import sys import svgwrite import numpy as np import matplotlib.pyplot as pp from matplotlib.colors import hsv_to_rgb import json import random +import argparse REF_LATITUDE = 49.58666 REF_LONGITUDE = 11.01250 -# REF_LATITUDE = 1 -# REF_LONGITUDE = 0 +# REF_LATITUDE = -30 +# REF_LONGITUDE = 90 def map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R=1): @@ -63,270 +65,287 @@ def random_country_color(): return f"#{r:02x}{g:02x}{b:02x}" -random.seed(0) +def render(ref_lat, ref_lon, output_stream): + 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) + """ 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) + 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(f"{lat*180/np.pi:6.3f}, {lon*180/np.pi:6.3f} => {x[0]:6.3f}, {y[0]:6.3f}", file=sys.stderr) + """ -print("Loading Geodata…") + print("Loading Geodata…", file=sys.stderr) -with open('geo-countries/data/countries.geojson', 'r') as jfile: - geojson = json.load(jfile) + with open('geo-countries/data/countries.geojson', 'r') as jfile: + geojson = json.load(jfile) -print("Finding boundaries…") + print("Finding boundaries…", file=sys.stderr) -# 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 = {} + # 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'] + features = geojson['features'] -for feature in features: - name = feature['properties']['ADMIN'] - key = feature['properties']['ISO_A2'] + 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 + # handle duplicate keys (can happen for small countries) + if key in simplegeodata.keys(): + key = name - print(f"Preparing {key} ({name})…") + print(f"Preparing {key} ({name})…", file=sys.stderr) - multipoly = feature['geometry']['coordinates'] + multipoly = feature['geometry']['coordinates'] - conv_polys = [] + conv_polys = [] - for poly in multipoly: - for subpoly in poly: - coords_list = [] # list of lists + for poly in multipoly: + for subpoly in poly: + coords_list = [] # list of lists - assert(len(subpoly[0]) == 2) - coords_list += subpoly + assert(len(subpoly[0]) == 2) + coords_list += subpoly - # convert coordinates to numpy array and radians - coords = np.array(coords_list).T * np.pi / 180 + # convert coordinates to numpy array and radians + coords = np.array(coords_list).T * np.pi / 180 - conv_polys.append(coords) + conv_polys.append(coords) - simplegeodata[key] = {"name": name, "coordinates": conv_polys} + simplegeodata[key] = {"name": name, "coordinates": conv_polys} -ref_lat = REF_LATITUDE * np.pi / 180 -ref_lon = REF_LONGITUDE * np.pi / 180 + ref_lat = ref_lat * np.pi / 180 + ref_lon = ref_lon * np.pi / 180 -R = 500 + R = 500 -""" -# Override data with test coordinate system -coords = [] + """ + # Override data with test coordinate system + coords = [] -N = 128 + 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-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)])) + # 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}} -""" + simplegeodata = {"XY": {'name': 'test', 'coordinates': coords}} + """ -# apply azimuthal equidistant projection -for k, v in simplegeodata.items(): - proj_polys = [] + # apply azimuthal equidistant projection + for k, v in simplegeodata.items(): + proj_polys = [] - for poly in v['coordinates']: - lat = poly[1, :] - lon = poly[0, :] + for poly in v['coordinates']: + lat = poly[1, :] + lon = poly[0, :] - x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) + x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) - coords = np.array([x, y]) + coords = np.array([x, y]) - # remove any points that contain a NaN coordinate - coords = coords[:, np.any(np.invert(np.isnan(coords)), axis=0)] + # remove any points that contain a NaN coordinate + coords = coords[:, np.any(np.invert(np.isnan(coords)), axis=0)] - proj_polys.append(coords) + proj_polys.append(coords) - v['proj_coordinates'] = proj_polys + v['proj_coordinates'] = proj_polys -# generate the SVG + # generate the SVG -doc = svgwrite.Drawing("/tmp/test.svg", size=(2*R, 2*R)) + doc = svgwrite.Drawing("/tmp/test.svg", size=(2*R, 2*R)) -doc.defs.add(doc.style(""" - .country { - stroke: black; - stroke-width: 0.01px; - } + doc.defs.add(doc.style(""" + .country { + stroke: black; + stroke-width: 0.01px; + } - .dist_circle_label, .azimuth_line_label { - font-size: 3px; - font: sans-serif; - text-anchor: middle; - } + .dist_circle_label, .azimuth_line_label { + font-size: 3px; + font: sans-serif; + text-anchor: middle; + } - .dist_circle, .azimuth_line { - fill: none; - stroke: black; - stroke-width: 0.1px; - } + .dist_circle, .azimuth_line { + fill: none; + stroke: black; + stroke-width: 0.1px; + } - .maidenhead_line { - fill: none; - stroke: red; - stroke-width: 0.1px; - opacity: 0.5; - } -""")) + .maidenhead_line { + fill: none; + stroke: red; + stroke-width: 0.1px; + opacity: 0.5; + } + """)) -doc.add(doc.circle(center=(R, R), r=R, fill='#ddeeff', - stroke_width=1, stroke='black')) + 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}…") + for k, v in simplegeodata.items(): + print(f"Exporting {k}…", file=sys.stderr) - color = random_country_color() + 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, + **{'class': 'dist_circle'})) + + doc.add(doc.text(f"{distance} km", (R, R-r+5), + **{'class': 'dist_circle_label'})) + + # generate azimuth lines in 30° steps + + for azimuth in np.arange(0, np.pi, np.pi/6): + start_x = R + R * np.cos(azimuth-np.pi/2) + start_y = R + R * np.sin(azimuth-np.pi/2) + end_x = R - R * np.cos(azimuth-np.pi/2) + end_y = R - R * np.sin(azimuth-np.pi/2) + + doc.add(doc.line((start_x, start_y), (end_x, end_y), + **{'class': 'azimuth_line'})) + + azimuth_deg = int(np.round(azimuth * 180 / np.pi)) + textpos = (2*R - 10, R - 2) + + txt = doc.text(f"{azimuth_deg:d} °", textpos, + **{'class': 'azimuth_line_label'}) + txt.rotate(azimuth_deg - 90, center=(R, R)) + doc.add(txt) + + txt = doc.text(f"{azimuth_deg+180:d} °", textpos, + **{'class': 'azimuth_line_label'}) + txt.rotate(azimuth_deg - 90 + 180, center=(R, R)) + doc.add(txt) + + # generate Maidenhead locator grid (first two letters only) group = doc.g() - for poly in v['proj_coordinates']: - points = poly.T + R # shift to the center of the drawing + N = 18 # subdivisions of Earth + resolution = 4096 - pgon = doc.polygon(points, **{ - 'class': 'country', - 'fill': color}) + for x in range(0, N): + lon = x * 2 * np.pi / N + lat = np.linspace(-np.pi/2, np.pi/2, resolution) - group.add(pgon) + x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) + points = np.array([x, y]).T + R + + group.add(doc.polyline(points, **{'class': 'maidenhead_line'})) + + for y in range(0, N): + lon = np.linspace(-np.pi, np.pi, resolution) + lat = y * np.pi / N - np.pi/2 + + x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) + points = np.array([x, y]).T + R + + group.add(doc.polyline(points, **{'class': 'maidenhead_line'})) - group.set_desc(title=v['name']) doc.add(group) -# generate equidistant circles + """ + for x in range(0, 26): + for y in range(0, 26): + sectorname = chr(ord('A')+x) + chr(ord('A')+y) + """ -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, - **{'class': 'dist_circle'})) + print(f"Writing output…", file=sys.stderr) + doc.write(output_stream, pretty=True) - doc.add(doc.text(f"{distance} km", (R, R-r+5), - **{'class': 'dist_circle_label'})) + return -# generate azimuth lines in 30° steps + # Debug Plot -for azimuth in np.arange(0, np.pi, np.pi/6): - start_x = R + R * np.cos(azimuth-np.pi/2) - start_y = R + R * np.sin(azimuth-np.pi/2) - end_x = R - R * np.cos(azimuth-np.pi/2) - end_y = R - R * np.sin(azimuth-np.pi/2) + for k, v in simplegeodata.items(): + for poly in v['proj_coordinates']: + pp.plot(poly[0, :], poly[1, :]) - doc.add(doc.line((start_x, start_y), (end_x, end_y), - **{'class': 'azimuth_line'})) + pp.plot([-1, 1], [0, 0], 'k', linewidth=0.5) + pp.plot([0, 0], [-1, 1], 'k', linewidth=0.5) - azimuth_deg = int(np.round(azimuth * 180 / np.pi)) - textpos = (2*R - 10, R - 2) + t = np.linspace(-np.pi, np.pi, 256) + ct, st = np.cos(t), np.sin(t) + pp.plot(ct, st, 'k', linewidth=0.5) - txt = doc.text(f"{azimuth_deg:d} °", textpos, - **{'class': 'azimuth_line_label'}) - txt.rotate(azimuth_deg - 90, center=(R, R)) - doc.add(txt) + 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) - txt = doc.text(f"{azimuth_deg+180:d} °", textpos, - **{'class': 'azimuth_line_label'}) - txt.rotate(azimuth_deg - 90 + 180, center=(R, R)) - doc.add(txt) + pp.axis('equal') -# generate Maidenhead locator grid (first two letters only) - -group = doc.g() - -N = 18 # subdivisions of Earth -resolution = 4096 - -for x in range(0, N): - lon = x * 2 * np.pi / N - lat = np.linspace(-np.pi/2, np.pi/2, resolution) - - x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) - points = np.array([x, y]).T + R - - group.add(doc.polyline(points, **{'class': 'maidenhead_line'})) - -for y in range(0, N): - lon = np.linspace(-np.pi, np.pi, resolution) - lat = y * np.pi / N - np.pi/2 - - x, y = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon, R) - points = np.array([x, y]).T + R - - group.add(doc.polyline(points, **{'class': 'maidenhead_line'})) - -doc.add(group) - -""" -for x in range(0, 26): - for y in range(0, 26): - sectorname = chr(ord('A')+x) + chr(ord('A')+y) -""" + pp.show() -print(f"Saving {doc.filename}…") -doc.save(pretty=True) +if __name__ == "__main__": + parser = argparse.ArgumentParser( + description="Render an azimuthal equidistant map of the world " + + "centered on the given point") -exit(0) + parser.add_argument(metavar='ref-lat', type=float, dest='ref_lat', + help='Reference Latitude') + parser.add_argument(metavar='ref-lon', type=float, dest='ref_lon', + help='Reference Longitude') + parser.add_argument('-o', '--output-file', type=argparse.FileType('w'), + help='The output SVG file (default: print to stdout)', + default=sys.stdout) -# 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() + args = parser.parse_args() + render(args.ref_lat, args.ref_lon, args.output_file)