From 60f59379e5388b1d3c86c5e699f91d77bf7fb7db Mon Sep 17 00:00:00 2001 From: Thomas Kolb Date: Thu, 10 Jun 2021 21:32:49 +0200 Subject: [PATCH] Refactoring: move processing steps to separate functions --- qsomap.py | 322 ++++++++++++++++++++++++++++-------------------------- 1 file changed, 170 insertions(+), 152 deletions(-) diff --git a/qsomap.py b/qsomap.py index a2ce520..de1c955 100755 --- a/qsomap.py +++ b/qsomap.py @@ -137,16 +137,7 @@ def svg_make_inverse_country_path(doc, map_radius, polygon, **kwargs): return doc.path(commands, **kwargs) -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) - +def simplify_geojson(geojson): # key: 3-letter country identifier # data: {full_name, # numpy.array(coordinates), @@ -187,15 +178,178 @@ def render(ref_lat, ref_lon, output_stream): simplegeodata[key] = {"name": name, "coordinates": conv_polys} - ref_lat = ref_lat * np.pi / 180 - ref_lon = ref_lon * np.pi / 180 + 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 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"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 """ @@ -233,24 +387,7 @@ def render(ref_lat, ref_lon, output_stream): 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 + map_all_polygons(simplegeodata, ref_lat, ref_lon, R) # 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', stroke_width=1, stroke='black')) - 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 + 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) - """ + svg_add_countries(doc, simplegeodata, ref_lat, ref_lon, 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)