diff --git a/server/map/__pycache__/get_relation.cpython-311.pyc b/server/map/__pycache__/get_relation.cpython-311.pyc
index 812af27f4b3e7b4696370419efaf670e47b9b2ff..76c8ff54db3931e76b092b850b48ecad05d4aa01 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 b17ab008349f74516ffb79ca63635334ce33f5d9..15d0590ec75cb8f1268d4316be59c42b1d1a8a88 100644
--- a/server/map/get_relation.py
+++ b/server/map/get_relation.py
@@ -13,7 +13,7 @@ def get_relation(self, body_of_water: str):
     # Filter only polygons, exclude points and other feature types to reduce response size
     polygon_data = geo_data[geo_data['geometry'].geom_type == 'Polygon']
 
-    if not polygon_data:
+    if polygon_data.empty:
         raise ValueError("Failed to extract polygon data from file")
 
     # Extract coordinates from polygons and create polygon objects
@@ -26,22 +26,30 @@ def get_relation(self, body_of_water: str):
 
     # Divide all polygons into sections
     for polygon in polygons:
-        # Initialise grid
-        grid_lines = create_grid(polygon, cell_size=0.1)
 
-        hrz_lines = grid_lines[0]  # Horizontal lines
-        vrt_lines = grid_lines[1]  # Vertical lines
+        # Divide the length and with of polygon into a grid
+        grid_lines = create_grid_coords(polygon, cell_size=0.1)
+
+        hrz_lines = grid_lines[0]  # Horizontal coordinates
+        vrt_lines = grid_lines[1]  # Vertical coordinates
 
         # Cut polygon into horizontal sections
         for hrz_line in hrz_lines:
-            # Splits shape into upper and lower section by hrz_line
-            cut_poly_1 = cut_polygon_in_two(polygon, hrz_line)
 
+            # Split shape into upper and lower section as hrz_line as divider
+            cut_poly_1 = cut_polygon_in_two(polygon, hrz_line, True)
             polygon_part = cut_poly_1[0]  # Save upper horizontal section
 
+            if not polygon_part or not cut_poly_1[0]:
+                continue
+
             # Cut each horizontal section into vertical sections
             for vrt_line in vrt_lines:
-                cut_poly_2 = cut_polygon_in_two(polygon_part, vrt_line)
+                cut_poly_2 = cut_polygon_in_two(polygon_part, vrt_line, False)
+
+                if not cut_poly_2[0]:
+                    continue
+
                 divided_map.extend(cut_poly_2[0])  # Append part to final list of shapes
 
                 # Set polygon_part to the remaining, un-split, horizontal section for next iteration
@@ -70,66 +78,50 @@ def get_relation(self, body_of_water: str):
     self.wfile.write(tiles_json.encode('utf-8'))
 
 
-def cut_polygon_in_two(polygon, line):
-    # Check for invalid parameters
-    if not isinstance(polygon, Polygon) or not isinstance(line, LineString):
-        print("Inputs must be Shapely Polygon and LineString obects")
-
+# Takes a polygon and divides its coordinates into two shapes, where divisor is a
+# coordinate that defines the point of division
+def cut_polygon_in_two(polygon: Polygon, divisor: float, horizontal: bool):
     # Extract polygon exterior coordinates
     exterior_coords = list(polygon.exterior.coords)
 
-    # Initialize lists to store points for the two shapes
-    shape_1 = []
-    shape_2 = []
+    # Initialize lists to store coordinates of new shapes after split
+    split_shape = []
+    remaining_shape = []
 
-    # Split points into separate shapes based on their location relative to the line
+    # Loop through points and check which side of the division line they are
     for point in exterior_coords:
         point = Point(point)  # Convert coordinates to Shapely Point object
-        if line.distance(point) < 1e-10:  # Floating point error tolerance
-            if line.contains(point):  # Determine if point belongs to shape 1 or shape 2
-                shape_1.append(point)
+        if horizontal:  # Horizontal split
+            if point.y < divisor:
+                split_shape.append(point)
             else:
-                shape_2.append(point)
-
-    # Check if shapes are empty
-    if not shape_1 or not shape_2:
-        raise ValueError("One or both shapes are empty after cutting")
-
-    # Create Polygon objects from the lists of points
-    poly_1 = Polygon(shape_1)
-    poly_2 = Polygon(shape_2)
+                remaining_shape.append(point)
+        else:  # Vertical split
+            if point.x < divisor:
+                split_shape.append(point)
+            else:
+                remaining_shape.append(point)
 
-    return poly_1, poly_2
+    # Check if polygons have enough coordinates
+    if len(split_shape) < 3 or len(remaining_shape) < 3:
+        print("Not enough coordinates to create valid polygons")
+        return None, None
 
+    # Append the first coordinate of the shapes to create closed loop
+    split_shape.append(split_shape[0]) # NB: may not be necessary?
+    remaining_shape.append(remaining_shape[0])
 
-# From https://stackoverflow.com/questions/39338550/cut-a-polygon-with-two-lines-in-shapely
-# Splits a polygon into two parts based on a dividing line
-def cut_polygon_by_line(polygon, 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)
+    return Polygon(split_shape), Polygon(remaining_shape)  # Return split polygons as Shapely Polygon objects
 
 
-# Create a grid of lines spaced by a cell size
-def create_grid(polygon, cell_size):
+# 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):
     min_x, min_y, max_x, max_y = polygon.bounds
     x_coords = np.arange(min_x, max_x, cell_size)
     y_coords = np.arange(min_y, max_y, cell_size)
-    horizontal_lines = []
-    vertical_lines = []
-
-    # Vertical lines from top to bottom
-    for x in x_coords:
-        vertical_lines.append(LineString([(x, min_y), (x, max_y)]))
 
-    # Horizontal lines from top to bottom
-    for y in y_coords:
-        horizontal_lines.append(LineString([(min_x, y), (max_x, y)]))
+    if len(x_coords) == 0 or len(y_coords) == 0:
+        raise ValueError("No grid points generated")
 
-    if not horizontal_lines or not vertical_lines:
-        raise ValueError("List of horizontal or vertical lines is empty")
+    return x_coords, y_coords
 
-    return horizontal_lines, vertical_lines