qsomap/qsomap.py

547 lines
17 KiB
Python
Executable file

#!/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
LABEL_MIN_FONT_SIZE = 2
LABEL_MAX_FONT_SIZE = 40
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 is_point_in_polygon(point, polygon):
# Idea: draw an infinite line from the test point along the x axis to the
# right. Then check how many polygon edges this line intersects. If the
# number is even, the point is outside the polygon.
edges = [] # list of lists, containing two points each
for i in range(len(polygon)-1):
edges.append([polygon[i], polygon[i+1]])
# the closing edge
edges.append([polygon[-1], polygon[0]])
num_intersects = 0
test_x, test_y = point
for edge in edges:
start_x = edge[0][0]
start_y = edge[0][1]
end_x = edge[1][0]
end_y = edge[1][1]
# quick exclusion tests
if start_x < test_x and end_x < test_x:
continue # edge is completely left of the test point
if start_y < test_y and end_y < test_y:
continue # edge is completely below the test point
if start_y > test_y and end_y > test_y:
continue # edge is completely above the test point
# calculate the x coordinate where the edge intersects the whole
# horizontal line
intersect_x = start_x + (end_x - start_x) \
* (test_y - start_y) \
/ (end_y - start_y)
if intersect_x > test_x:
# we found an intersection!
num_intersects += 1
if num_intersects % 2 == 0:
return False # even number of intersects -> outside polygon
else:
return True # odd number of intersects -> inside polygon
def polygon_centroid(polygon):
""" Calculate the centroid (center of mass) of the given 2D polygon.
See https://en.wikipedia.org/wiki/Centroid#Of_a_polygon for the source
formulae.
Args:
polygon: 2xN numpy array with x and y coordinates of N points.
Returns:
(x, y): The coordinates of the center of mass.
"""
# Separate the coordinates
x = polygon[0, :]
y = polygon[1, :]
# First, calculate the area of the polygon
A = 0.5 * np.sum(x[:-1] * y[1:] - x[1:] * y[:-1])
# Calculate x and y of the centroid
C_x = np.sum((x[:-1] + x[1:])
* (x[:-1] * y[1:] - x[1:] * y[:-1])) / (6 * A)
C_y = np.sum((y[:-1] + y[1:])
* (x[:-1] * y[1:] - x[1:] * y[:-1])) / (6 * A)
return C_x, C_y
def svg_make_inverse_country_path(doc, map_radius, polygon, **kwargs):
# build a closed circle path covering the whole map
commands = [f"M 0, {map_radius}",
f"a {map_radius},{map_radius} 0 1,0 {map_radius*2},0",
f"a {map_radius},{map_radius} 0 1,0 {-map_radius*2},0",
"z"]
# "subtract" the country polygon
commands.append(f"M {polygon[0][0]} {polygon[0][1]}")
# add lines for each polygon point
for point in polygon[1:]:
commands.append(f"L {point[0]} {point[1]}")
# ensure straight closing line
commands.append(f"L {polygon[0][0]} {polygon[0][1]}")
# close the inner path
commands.append("z")
return doc.path(commands, **kwargs)
def simplify_geojson(geojson):
# 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})…", file=sys.stderr)
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}
return simplegeodata
def map_all_polygons(simplegeodata, ref_lat, ref_lon, map_radius):
# 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,
map_radius)
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
def hsv_to_svgstr(h, s, v):
r, g, b = [int(255.99*x) for x in hsv_to_rgb([h, s, v])]
return f"#{r:02x}{g:02x}{b:02x}"
def assign_country_colors(simplegeodata):
for k, v in simplegeodata.items():
hue = random.random()
v['polygon_color'] = hsv_to_svgstr(hue, 0.7, 0.8)
v['label_color'] = hsv_to_svgstr(hue, 0.5, 0.4)
def svg_add_countries(doc, simplegeodata, ref_lat, ref_lon, map_radius):
antipodal_lat = -ref_lat
antipodal_lon = ref_lon + np.pi
if antipodal_lon > np.pi:
antipodal_lon -= 2*np.pi
for k, v in simplegeodata.items():
print(f"Mapping {k}", file=sys.stderr)
color = v['polygon_color']
group = doc.g()
for i in range(len(v['proj_coordinates'])):
poly = v['proj_coordinates'][i]
points = poly.T + map_radius # shift to the center of the drawing
# check if the antipodal point is inside this polygon. If so, it
# needs to be "inverted", i.e. subtracted from the surrounding map
# circle.
if is_point_in_polygon((antipodal_lon, antipodal_lat),
v['coordinates'][i].T):
print("!!! Found polygon containing the antipodal point!",
file=sys.stderr)
obj = svg_make_inverse_country_path(doc, map_radius,
np.flipud(points),
**{'class': 'country',
'fill': color})
else:
obj = doc.polygon(points, **{
'class': 'country',
'fill': color})
group.add(obj)
group.set_desc(title=v['name'])
doc.add(group)
def svg_add_maidenhead_grid(doc, ref_lat, ref_lon, map_radius):
# 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, map_radius)
points = np.array([x, y]).T + map_radius
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,
map_radius)
points = np.array([x, y]).T + map_radius
group.add(doc.polyline(points, **{'class': 'maidenhead_line'}))
for x in range(0, N):
for y in range(0, N):
sectorname = chr(ord('A') + (x + N//2) % N) \
+ chr(ord('A') + y)
lon = (x + 0.5) * 2 * np.pi / N
lat = (y + 0.5) * np.pi / N - np.pi/2
tx, ty = map_azimuthal_equidistant(lat, lon, ref_lat, ref_lon,
map_radius)
font_size = 10
if y == 0 or y == N-1:
font_size = 3
group.add(doc.text(sectorname, (tx + map_radius, ty + map_radius),
**{'class': 'maidenhead_label',
'font-size': font_size}))
doc.add(group) # Maidenhead grid
def svg_add_country_names(doc, simplegeodata, map_radius):
group = doc.g()
for k, v in simplegeodata.items():
print(f"Labeling {k} ", end='')
for poly in v['proj_coordinates']:
x = poly[0, :]
y = poly[1, :]
w = np.max(x) - np.min(x)
h = np.max(y) - np.min(y)
# align text at the center of mass
center_x, center_y = polygon_centroid(poly)
center_x += map_radius
center_y += map_radius
kwargs = {
'class': 'country_label',
'fill': v['label_color']
}
label = v['name']
# rotate text by 90 degrees if polygon is higher than wide
if h > w:
font_size = h / len(label)
rotate = True
else:
font_size = w / len(label)
rotate = False
if font_size < LABEL_MIN_FONT_SIZE:
# try again with 2-letter country code
label = k
font_size = min(LABEL_MIN_FONT_SIZE*2,
(h if rotate else w) / len(label))
if font_size < LABEL_MIN_FONT_SIZE:
print('.', end='', flush=True)
continue # too small
else:
print('!', end='', flush=True)
else:
print('#', end='', flush=True)
font_size = min(LABEL_MAX_FONT_SIZE, font_size)
kwargs['font-size'] = font_size
txt = doc.text(label, (center_x, center_y),
**kwargs)
if rotate:
txt.rotate(90, center=(center_x, center_y))
group.add(txt)
print()
doc.add(group)
def svg_add_distance_azimuth_lines(doc, ref_lat, ref_lon, map_radius):
group = doc.g()
# 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 = map_radius * distance / d_max
group.add(doc.circle(center=(map_radius, map_radius), r=r,
**{'class': 'dist_circle'}))
group.add(doc.text(f"{distance} km", (map_radius, map_radius-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 = map_radius + map_radius * np.cos(azimuth-np.pi/2)
start_y = map_radius + map_radius * np.sin(azimuth-np.pi/2)
end_x = map_radius - map_radius * np.cos(azimuth-np.pi/2)
end_y = map_radius - map_radius * np.sin(azimuth-np.pi/2)
group.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*map_radius - 10, map_radius - 2)
txt = doc.text(f"{azimuth_deg:d} °", textpos,
**{'class': 'azimuth_line_label'})
txt.rotate(azimuth_deg - 90, center=(map_radius, map_radius))
group.add(txt)
txt = doc.text(f"{azimuth_deg+180:d} °", textpos,
**{'class': 'azimuth_line_label'})
txt.rotate(azimuth_deg - 90 + 180, center=(map_radius, map_radius))
group.add(txt)
doc.add(group) # Circles, azimuth lines and labels
def render(ref_lat, ref_lon, output_stream):
random.seed(0)
print("Loading Geodata…", file=sys.stderr)
with open('geo-countries/data/countries.geojson', 'r') as jfile:
geojson = json.load(jfile)
print("Finding boundaries…", file=sys.stderr)
simplegeodata = simplify_geojson(geojson)
ref_lat = ref_lat * np.pi / 180
ref_lon = ref_lon * np.pi / 180
R = 500
map_all_polygons(simplegeodata, ref_lat, ref_lon, R)
assign_country_colors(simplegeodata)
# 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.01px;
}
.dist_circle_label, .azimuth_line_label {
font-size: 3px;
font-family: sans-serif;
text-anchor: middle;
}
.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_label, .country_label {
font-family: sans-serif;
dominant-baseline: middle;
text-anchor: middle;
}
.maidenhead_label {
fill: red;
opacity: 0.25;
}
"""))
doc.add(doc.circle(center=(R, R), r=R, fill='#ddeeff',
stroke_width=1, stroke='black'))
svg_add_countries(doc, simplegeodata, ref_lat, ref_lon, R)
svg_add_country_names(doc, simplegeodata, R)
svg_add_maidenhead_grid(doc, ref_lat, ref_lon, R)
svg_add_distance_azimuth_lines(doc, ref_lat, ref_lon, R)
print("Writing output…", file=sys.stderr)
doc.write(output_stream, pretty=True)
return
# 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()
if __name__ == "__main__":
parser = argparse.ArgumentParser(
description="Render an azimuthal equidistant map of the world " +
"centered on the given point")
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)
args = parser.parse_args()
render(args.ref_lat, args.ref_lon, args.output_file)