diff --git a/qsomap.py b/qsomap.py index 7cc4104..927bd83 100755 --- a/qsomap.py +++ b/qsomap.py @@ -65,6 +65,78 @@ def random_country_color(): return f"#{r:02x}{g:02x}{b:02x}" +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 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 render(ref_lat, ref_lon, output_stream): random.seed(0) @@ -127,6 +199,12 @@ def render(ref_lat, ref_lon, output_stream): ref_lat = ref_lat * np.pi / 180 ref_lon = ref_lon * np.pi / 180 + antipodal_lat = -ref_lat + antipodal_lon = ref_lon + np.pi + + if antipodal_lon > np.pi: + antipodal_lon -= 2*np.pi + R = 500 """ @@ -183,7 +261,6 @@ def render(ref_lat, ref_lon, output_stream): v['proj_coordinates'] = proj_polys - # generate the SVG doc = svgwrite.Drawing("/tmp/test.svg", size=(2*R, 2*R)) @@ -224,14 +301,26 @@ def render(ref_lat, ref_lon, output_stream): group = doc.g() - for poly in v['proj_coordinates']: + for i in range(len(v['proj_coordinates'])): + poly = v['proj_coordinates'][i] points = poly.T + R # shift to the center of the drawing - pgon = doc.polygon(points, **{ - 'class': 'country', - 'fill': color}) + # check if the antipodal point is inside this polygon. If so, it + # needs to be "inverted", i.e. subtracted from the surrounding map + # circle. - group.add(pgon) + if is_point_in_polygon((antipodal_lon, antipodal_lat), + v['coordinates'][i].T): + print("!!! Found polygon containing the antipodal point!") + 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)