Support for command line arguments
This commit is contained in:
parent
af794b0598
commit
a766274d79
429
qsomap.py
429
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)
|
||||
|
|
|
|||
Loading…
Reference in a new issue