From a9f987c07964e55a75e09049e45b654342eae16e Mon Sep 17 00:00:00 2001
From: Thomas Kolb <cfr34k-git@tkolb.de>
Date: Tue, 1 Jun 2021 22:58:33 +0200
Subject: [PATCH] Fix rendering with land at the antipodal point
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

If the antipodal point (the point opposite of the map center) is on land,
the corresponding land mass was stretched across the whole map, breaking
the rendering. This was fixed by a couple of changes:

1. A check is done for each polygon whether it contains the antipodal
   point.
2. If it does not, render it as a simple polygon as it was previously.
3. If it does contain the antipodal point, create a path that contains:
   a) A circle covering the whole map.
   b) The polygon with the points in reverse order.

This leads to filled circle with the polygon “cut out”, which is the
correct display in this case.
---
 qsomap.py | 101 ++++++++++++++++++++++++++++++++++++++++++++++++++----
 1 file changed, 95 insertions(+), 6 deletions(-)

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)