update: second split method

server/map/get_relation.py

 import geopandas as gpd
-from shapely.geometry import Polygon, Point, LineString, MultiLineString, MultiPolygon
+from shapely.geometry import Polygon, Point, LineString
 import matplotlib.pyplot as plt
 from shapely.ops import linemerge, unary_union, polygonize
 import matplotlib.ticker as ticker
 import numpy as np
-import pandas as pd
 
 
 def get_relation(self, body_of_water: str):
@@ -22,32 +21,22 @@ def get_relation(self, body_of_water: str):
         print("Failed to convert to polygons")
         return
 
-    ########################################################################################
     print("Performing div call")
-    # NB: any cell size under 0.02 will require patience
-    # For example, 0.007 requires up to 3 minutes to compute
-    # tiles = divide_relation(polygons, 0.02)
-
     test_line = LineString([(10.4600, 60.7451), (11.2000, 60.7451)])
-    line_pol = test_line.buffer(1e-3)
 
-    new_polygon = polygons[0].difference(line_pol)
+    new_polygons = []
+    for polygon in polygons:
+        new_polygon = cut_polygon_by_line(polygon, test_line)
+        new_polygons.extend(new_polygon)
 
-    tiles = gpd.GeoDataFrame(geometry=[new_polygon])
+    tiles = gpd.GeoDataFrame(geometry=new_polygons)
 
     fig, ax = plt.subplots()
-
     ax.set_aspect(1.5)
     ax.xaxis.set_major_locator(ticker.MultipleLocator(0.2))
     tiles.plot(ax=ax, color='blue', edgecolor='blue', alpha=0.5, label='Original Polygon')
     plt.show()
 
-    if test_line.intersects(polygons[0].boundary):
-        print("Intersects!")
-    else:
-        print("Test line does not intersect the polygon")
-
-    ########################################################################################
     # Convert GeoDataFrame to GeoJSON
     tiles_json = tiles.to_json()
 
@@ -60,47 +49,12 @@ def get_relation(self, body_of_water: str):
     self.wfile.write(tiles_json.encode('utf-8'))
 
 
-######################### Approach 1 #################################
-def divide_relation(polygons, cell_size):
-    # Calculate the combined bounding box of all polygons
-    combined_bounds = polygons[0].bounds
-    for polygon in polygons[1:]:
-        combined_bounds = (min(combined_bounds[0], polygon.bounds[0]),
-                           min(combined_bounds[1], polygon.bounds[1]),
-                           max(combined_bounds[2], polygon.bounds[2]),
-                           max(combined_bounds[3], polygon.bounds[3]))
-
-    # Create an empty list to store polygons representing tiles
-    tiles = []
-
-    # Iterate over each polygon
-    for polygon in polygons:
-        # Iterate over each cell
-        for i in np.arange(combined_bounds[1], combined_bounds[3], cell_size):
-            for j in np.arange(combined_bounds[0], combined_bounds[2], cell_size):
-                # Define the bounds of the cell
-                cell_bounds = (j, i, j + cell_size, i + cell_size)
-
-                # Create a polygon representing the cell
-                cell_polygon = Polygon([(cell_bounds[0], cell_bounds[1]),
-                                        (cell_bounds[2], cell_bounds[1]),
-                                        (cell_bounds[2], cell_bounds[3]),
-                                        (cell_bounds[0], cell_bounds[3])])
-
-                # Check if the cell polygon intersects with the current polygon
-                if polygon.intersects(cell_polygon):
-                    tiles.append(cell_polygon)
-
-    # Create a GeoDataFrame from the list of tiles
-    tiles_gdf = gpd.GeoDataFrame(geometry=tiles)
-
-    return tiles_gdf
-
-
-######################### Approach 2 #################################
-
+# From https://stackoverflow.com/questions/39338550/cut-a-polygon-with-two-lines-in-shapely
 def cut_polygon_by_line(polygon, line):
-    merged = linemerge([polygon.boundary, line])
+    if polygon.geom_type == 'Polygon':
+        polygon = [polygon]
+    merged = linemerge([poly.boundary for poly in polygon] + [line])
     borders = unary_union(merged)
     polygons = polygonize(borders)
     return list(polygons)
+