diff --git a/server/map/__pycache__/get_relation.cpython-311.pyc b/server/map/__pycache__/get_relation.cpython-311.pyc
index b60cf86f63fea40cec5c55d8ee98c439fd98b3da..68c0cc6a2c8de1d9bad2c4a677526c0144087a78 100644
Binary files a/server/map/__pycache__/get_relation.cpython-311.pyc and b/server/map/__pycache__/get_relation.cpython-311.pyc differ
diff --git a/server/map/get_relation.py b/server/map/get_relation.py
index 4d32b06fa0c21e602b43acceefa5b50cb37b4e8c..81f7ce5f73a028319021e576f52c5c58fa75227e 100644
--- a/server/map/get_relation.py
+++ b/server/map/get_relation.py
@@ -96,46 +96,42 @@ def cut_polygon_by_points(polygon: Polygon, divisor: float, cell_size: float):
     remaining_shape = []
 
     # Loop through points and check which side of the division line they are
-    if cell_size > 0:  # Horizontal split
-        for point in exterior_coords:
-            point = Point(point)  # Convert coordinates to Shapely Point object
-            if point.y < divisor:  # Check if point is over or below divisor
-                split_shape.append(point)  # Append to appropriate shape
+    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(point)
-
-        # Populate newly created edges with points to facilitate vertical cutting
-        populated_shapes = populate_new_edge(split_shape, remaining_shape, cell_size, divisor)
-        if populated_shapes is not None:
-            split_shape = populated_shapes[0]
-            remaining_shape = populated_shapes[1]
-
-    else:  # Vertical split
-        for point in exterior_coords:
-            point = Point(point)  # Convert coordinates to Shapely Point object
-            if point.x < divisor:
-                split_shape.append(point)
+                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(point)
-
-    # Check if the split_shape has enough coordinates to create a polygon
-    if len(split_shape) < 3:
-        # print("Not enough coordinates to create valid polygon: Split shape: ", len(split_shape))
-        split_shape = None
-    else:
-        split_shape.append(split_shape[0])  # Append first coord to create closed loop
-
-    # Check if the remaining_shape has enough coordinates to create a polygon
-    if len(remaining_shape) < 3:
-        # print("Not enough coordinates to create valid polygon: Remaining shape: ", len(remaining_shape))
-        remaining_shape = None
-    else:
-        remaining_shape.append(remaining_shape[0])  # Append first coord to create closed loop
-
-    # Return split polygons as Shapely Polygon objects
+                remaining_shape.append((point2.x, point2.y))
+
+    # 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
@@ -178,20 +174,17 @@ def populate_new_edge(split_shape, remaining_shape, cell_size: float, divisor: f
         # Define starting point with an x-value that will be common for all polygons
         starting_point = polygon_min_x
 
-        print("Hit main if")
         # Create list of corners
         corners = []
 
-        divisor_range = ((divisor - 0.1),(divisor + 0.1)) # Define tolerance
+        divisor_range = ((divisor - 0.1),(divisor + 0.1))  # Define tolerance NB: must find appropriate value
 
         # Find all corners of split shape. Corner = point that intersects divisor
         for point in split_shape:
             if divisor_range[0] < point.y < divisor_range[1]:
                 corners.append(point)
-                print("Hit point in split_shape")
 
         if not corners or len(corners) < 2:
-            print("no corners :(")
             return None
 
         # Sort corners in ascending order (left to right) based on x coordinate
@@ -207,7 +200,6 @@ def populate_new_edge(split_shape, remaining_shape, cell_size: float, divisor: f
         while starting_point < sorted_corners[-1].x:
             split_shape.insert(0, (starting_point, divisor))  # NB may have to add/subtract small offset of 0.00001
             remaining_shape.insert(-1, (starting_point, divisor))  # Prepend new point to shape
-            print("Hit point insertion")
 
             starting_point += cell_size