Refactoring: move processing steps to separate functions
This commit is contained in:
parent
67caf16b14
commit
60f59379e5
322
qsomap.py
322
qsomap.py
|
|
@ -137,16 +137,7 @@ def svg_make_inverse_country_path(doc, map_radius, polygon, **kwargs):
|
||||||
return doc.path(commands, **kwargs)
|
return doc.path(commands, **kwargs)
|
||||||
|
|
||||||
|
|
||||||
def render(ref_lat, ref_lon, output_stream):
|
def simplify_geojson(geojson):
|
||||||
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)
|
|
||||||
|
|
||||||
# key: 3-letter country identifier
|
# key: 3-letter country identifier
|
||||||
# data: {full_name,
|
# data: {full_name,
|
||||||
# numpy.array(coordinates),
|
# numpy.array(coordinates),
|
||||||
|
|
@ -187,15 +178,178 @@ def render(ref_lat, ref_lon, output_stream):
|
||||||
|
|
||||||
simplegeodata[key] = {"name": name, "coordinates": conv_polys}
|
simplegeodata[key] = {"name": name, "coordinates": conv_polys}
|
||||||
|
|
||||||
ref_lat = ref_lat * np.pi / 180
|
return simplegeodata
|
||||||
ref_lon = ref_lon * np.pi / 180
|
|
||||||
|
|
||||||
|
|
||||||
|
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 svg_add_countries(doc, simplegeodata, ref_lat, ref_lon, map_radius):
|
||||||
antipodal_lat = -ref_lat
|
antipodal_lat = -ref_lat
|
||||||
antipodal_lon = ref_lon + np.pi
|
antipodal_lon = ref_lon + np.pi
|
||||||
|
|
||||||
if antipodal_lon > np.pi:
|
if antipodal_lon > np.pi:
|
||||||
antipodal_lon -= 2*np.pi
|
antipodal_lon -= 2*np.pi
|
||||||
|
|
||||||
|
for k, v in simplegeodata.items():
|
||||||
|
print(f"Exporting {k}…", file=sys.stderr)
|
||||||
|
|
||||||
|
color = random_country_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_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
|
R = 500
|
||||||
|
|
||||||
"""
|
"""
|
||||||
|
|
@ -233,24 +387,7 @@ def render(ref_lat, ref_lon, output_stream):
|
||||||
simplegeodata = {"XY": {'name': 'test', 'coordinates': coords}}
|
simplegeodata = {"XY": {'name': 'test', 'coordinates': coords}}
|
||||||
"""
|
"""
|
||||||
|
|
||||||
# apply azimuthal equidistant projection
|
map_all_polygons(simplegeodata, ref_lat, ref_lon, R)
|
||||||
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
|
|
||||||
|
|
||||||
# generate the SVG
|
# generate the SVG
|
||||||
|
|
||||||
|
|
@ -294,128 +431,9 @@ def render(ref_lat, ref_lon, output_stream):
|
||||||
doc.add(doc.circle(center=(R, R), r=R, fill='#ddeeff',
|
doc.add(doc.circle(center=(R, R), r=R, fill='#ddeeff',
|
||||||
stroke_width=1, stroke='black'))
|
stroke_width=1, stroke='black'))
|
||||||
|
|
||||||
for k, v in simplegeodata.items():
|
svg_add_countries(doc, simplegeodata, ref_lat, ref_lon, R)
|
||||||
print(f"Exporting {k}…", file=sys.stderr)
|
svg_add_maidenhead_grid(doc, ref_lat, ref_lon, R)
|
||||||
|
svg_add_distance_azimuth_lines(doc, ref_lat, ref_lon, R)
|
||||||
color = random_country_color()
|
|
||||||
|
|
||||||
group = doc.g()
|
|
||||||
|
|
||||||
for i in range(len(v['proj_coordinates'])):
|
|
||||||
poly = v['proj_coordinates'][i]
|
|
||||||
points = poly.T + R # 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, R, 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)
|
|
||||||
|
|
||||||
# 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'}))
|
|
||||||
|
|
||||||
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, R)
|
|
||||||
|
|
||||||
font_size = 10
|
|
||||||
if y == 0 or y == N-1:
|
|
||||||
font_size = 3
|
|
||||||
|
|
||||||
group.add(doc.text(sectorname, (tx + R, ty + R),
|
|
||||||
**{'class': 'maidenhead_label',
|
|
||||||
'font-size': font_size}))
|
|
||||||
|
|
||||||
doc.add(group) # Maidenhead grid
|
|
||||||
|
|
||||||
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 = R * distance / d_max
|
|
||||||
group.add(doc.circle(center=(R, R), r=r,
|
|
||||||
**{'class': 'dist_circle'}))
|
|
||||||
|
|
||||||
group.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)
|
|
||||||
|
|
||||||
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*R - 10, R - 2)
|
|
||||||
|
|
||||||
txt = doc.text(f"{azimuth_deg:d} °", textpos,
|
|
||||||
**{'class': 'azimuth_line_label'})
|
|
||||||
txt.rotate(azimuth_deg - 90, center=(R, R))
|
|
||||||
group.add(txt)
|
|
||||||
|
|
||||||
txt = doc.text(f"{azimuth_deg+180:d} °", textpos,
|
|
||||||
**{'class': 'azimuth_line_label'})
|
|
||||||
txt.rotate(azimuth_deg - 90 + 180, center=(R, R))
|
|
||||||
group.add(txt)
|
|
||||||
|
|
||||||
doc.add(group) # Circles, azimuth lines and labels
|
|
||||||
|
|
||||||
"""
|
|
||||||
for x in range(0, 26):
|
|
||||||
for y in range(0, 26):
|
|
||||||
sectorname = chr(ord('A')+x) + chr(ord('A')+y)
|
|
||||||
"""
|
|
||||||
|
|
||||||
print("Writing output…", file=sys.stderr)
|
print("Writing output…", file=sys.stderr)
|
||||||
doc.write(output_stream, pretty=True)
|
doc.write(output_stream, pretty=True)
|
||||||
|
|
|
||||||
Loading…
Reference in a new issue