update: backtrack to method 2

server/map/get_relation.py

 import geopandas as gpd
-from shapely.geometry import Polygon, Point, LineString, MultiPolygon
+from shapely.geometry import Polygon, LineString
+from shapely.ops import linemerge, unary_union, polygonize
 import matplotlib.pyplot as plt
 import matplotlib.ticker as ticker
 import random
@@ -18,44 +19,14 @@ def get_relation(self, body_of_water: str):  # NB: implement body_of_water
 
     # Divide all polygons into sections
     for polygon in polygons:
-
-        cell_size = 0.4
+        cell_size = 1
         # Divide the length and with of polygon into a grid of equally sized parts
-        grid_lines = create_grid_coords(polygon, cell_size)
-
-        vrt_lines = grid_lines[0]  # Vertical grid coordinates
-        hrz_lines = grid_lines[1]  # Horizontal grid coordinates
-
-        # Cut polygon into horizontal sections, bottom to top
-        for hrz_line in hrz_lines:
-
-            # Split shape into upper and lower section as hrz_line as divider
-            divided_poly = cut_polygon_by_points(polygon, hrz_line, cell_size)
-            horizontal_section = divided_poly[0]  # Save upper horizontal section
-
-            polygon = divided_poly[1]  # Set polygon to the remaining, un-split shape for next iteration
-
-            if not horizontal_section or not divided_poly[0]:
-                continue
-
-                # Cut each horizontal section into vertical sections, right to left
-            for vrt_line in vrt_lines:
-                if len(horizontal_section.exterior.coords) < 3:
-                    break  # Break from loop im remaining section has no coordinates
-
-                # Split the horizontal section into two vertical parts
-                vertical_parts = cut_polygon_by_points(horizontal_section, vrt_line, -0.1)
-
-                divided_map.append(vertical_parts[0])  # Append split vertical sections to final list of shapes
+        lines = create_grid(polygon, cell_size)
 
-                # Set horizontal_section to the remaining, un-split, horizontal section for next iteration
-                horizontal_section = vertical_parts[1]
-
-            polygon = divided_poly[1]  # Set polygon to the remaining, un-split shape for next iteration
-
-        divided_map.append(polygon)
-
-        break
+        lines.append(polygon.boundary)
+        lines = unary_union(lines)
+        lines = linemerge(lines)
+        divided_map.extend(list(polygonize(lines)))
 
     tiles = gpd.GeoDataFrame(geometry=divided_map)
 
@@ -85,77 +56,6 @@ def get_relation(self, body_of_water: str):  # NB: implement body_of_water
     self.wfile.write(tiles_json.encode('utf-8'))
 
 
-def cut_polygon_by_points(polygon: Polygon, divisor: float, cell_size: float):
-    # Extract polygon exterior coordinates
-    exterior_coords = list(polygon.exterior.coords)
-
-    # Initialize lists to store coordinates of new shapes after split
-    split_shape = []
-    remaining_shape = []
-
-    # Loop through points and check which side of the division line they are
-    for i in range(len(exterior_coords) - 1):
-        point1 = Point(exterior_coords[i])
-        point2 = Point(exterior_coords[i + 1])
-
-        # Check if both points are on the same side of the divisor
-        if (point1.y < divisor and point2.y < divisor) or (point1.y >= divisor and point2.y >= divisor):
-            # Both points on same side, add to appropriate shape
-            if point1.y < divisor:
-                split_shape.append(exterior_coords[i])
-            else:
-                remaining_shape.append(exterior_coords[i])
-        else:
-            # Points on different sides, calculate intersection with divisor
-            x_intersect = point1.x + (divisor - point1.y) * (point2.x - point1.x) / (point2.y - point1.y)
-
-            # Add intersection point to both shapes
-            split_shape.append((x_intersect, divisor))
-            remaining_shape.append((x_intersect, divisor))
-
-            # Determine which side each original point belongs to
-            if point1.y < divisor:
-                split_shape.append((point2.x, point2.y))
-            else:
-                remaining_shape.append((point2.x, point2.y))
-
-    # Populate newly created edges with points to facilitate vertical cutting
-    populated_split = populate_new_edge(split_shape, cell_size, divisor)
-    if populated_split is not None:
-        split_shape = populated_split
-
-    populated_rem = populate_new_edge(remaining_shape, cell_size, divisor)
-    if populated_rem is not None:
-        remaining_shape = populated_rem
-
-    # Ensure that the shapes are closed loops
-    if len(split_shape) > 0 and split_shape[0] != split_shape[-1]:
-        split_shape.append(split_shape[0])
-
-    if len(remaining_shape) > 0 and remaining_shape[0] != remaining_shape[-1]:
-        remaining_shape.append(remaining_shape[0])
-
-    return Polygon(split_shape), Polygon(remaining_shape)
-
-
-# Generate grid of equally spaced x and y coordinates where the grid size is determined by cell_size
-def create_grid_coords(polygon: Polygon, cell_size: float):
-    # Define boundaries of grid
-    min_x, min_y, max_x, max_y = polygon.bounds
-
-    global polygon_min_x
-    polygon_min_x = min_x
-
-    # Divide grid into sections of size *cell_size
-    x_step = max((max_x - min_x) / 10, cell_size)
-    y_step = max((max_y - min_y) / 10, cell_size)
-
-    x_coords = np.arange(min_x, max_x + x_step, x_step)
-    y_coords = np.arange(min_y, max_y + y_step, y_step)
-
-    return x_coords, y_coords
-
-
 # NB: only for testing
 def circle_polygon():
     circle_points = []
@@ -173,42 +73,25 @@ def circle_polygon():
     return Polygon(circle_points)
 
 
-def populate_new_edge(polygon, cell_size: float, divisor: float):
-    if polygon is not None and polygon_min_x is not None:
-
-        final_shape = []
-        left_edge = []
-        left_most_point = None
-
-        # Add all coordinates on the right side of the divisor first
-        for i in range(len(polygon)):
-            x, y = polygon[i]
-            if x > divisor:
-                final_shape.append((x, y))
-            elif x < divisor:
-                left_edge.append((x, y))
-
-        if len(final_shape) < 1 or len(left_edge) < 1:
-            return None
-
-        # Find most left point in the polygon
-        for point in left_edge:
-            if not left_most_point:
-                left_most_point = point[0]
-            elif point[0] < left_most_point:
-                left_most_point = point[0]
+def create_grid(poly: Polygon, cell_size):
+    # Retrieve bounds of the entire polygon
+    bounds = poly.bounds
 
-        # Sort final_shape from bottom to top
-        if final_shape[0][1] > final_shape[-1][1]:
-            final_shape = sorted(final_shape, key=lambda point: final_shape[1])
+    min_x, min_y, max_x, max_y = bounds
+    grid_lines = []
 
-        i = 0
-        while polygon[i][0] < left_most_point:  # While to haven't reached left edge
-            final_shape.append(polygon[i][0] + cell_size)
-            i += 1
+    # Horizontal lines
+    y = min_y
+    while y <= max_y:
+        line = LineString([(min_x, y), (max_x, y)])
+        grid_lines.append(line)
+        y += cell_size
 
-        # NB VERY IMPORTANT MUST CREATE LOGIC FOR SHAPES THAT ARE SPLIT
-        final_shape.append(left_edge)
-        final_shape.append(final_shape[0])
+    # Vertical lines
+    x = min_x
+    while x <= max_x:
+        line = LineString([(x, min_y), (x, max_y)])
+        grid_lines.append(line)
+        x += cell_size
 
-        return final_shape
+    return grid_lines