diff --git a/server/map_handler/add_new_lake.py b/server/map_handler/add_new_lake.py
index 46a2a4809c573466b97acb25efef456fb783e928..90c8d2484a5f2c3a2eb06784ea5b7188d0862d26 100644
--- a/server/map_handler/add_new_lake.py
+++ b/server/map_handler/add_new_lake.py
@@ -29,20 +29,21 @@ def cut_map(self, cursor, lake_name: str, cell_size_in_km: float = 0.5):
         if len(polygons) <= 1:
             raise Exception("Failed to convert JSON object to Shapely Polygons")
 
-        # List to store all the map tiles
-        divided_map = []
-
         # Select an arbitrary x and y value from within the polygon
         bounds = polygons[0].bounds
         start_x, start_y, _, _ = bounds
 
         cell_width = cell_size_in_km / 111.3200
         cell_height = cell_width / cos(start_x * 0.01745)
+        # List to store new GeoJSON feature objects
+        features = []   # List to store new GeoJSON feature objects
+        sub_div_id = 0  # Tracker to create unique subdivision ids
+        divided_map = []    # Object for plotting the tiles
 
         # Process all polygons
         for polygon in polygons:
             # Generate a grid based on the calculated cell size
-            lines = create_grid(polygon, cell_width*2, cell_height)
+            lines = create_grid(polygon, cell_width * 2, cell_height)
             lines.append(polygon.boundary)
             # Merge the grid lines into a single grid object
             lines = unary_union(lines)
@@ -51,48 +52,40 @@ def cut_map(self, cursor, lake_name: str, cell_size_in_km: float = 0.5):
 
             # Combine the polygon and the grid to form the subdivisions
             for line in lines:
-                if line.intersects(polygon):
-                    divided_map.append(line.intersection(polygon))
-
+                if line.intersects(polygon):  # Add the geometry which intersects the gird
+                    tile = (line.intersection(polygon))
+                    divided_map.append(tile)
+
+                    # Calculate tile center based on bounds, and round down to two decimals
+                    min_x, min_y, max_x, max_y = tile.bounds
+                    center = round(max_y - (max_y - min_y), 6), round(max_x - (max_x - min_x), 6)
+
+                    rounded_coordinates = []
+                    if isinstance(tile, Polygon):
+                        for coords in tile.exterior.coords:
+                            rounded_coords = (round(coords[0], 4), round(coords[1], 4))
+                            rounded_coordinates.append(rounded_coords)
+
+                    rounded_tile = Polygon(rounded_coordinates)
+                    geometry = rounded_tile.__geo_interface__
+
+                    if not geometry['coordinates']:
+                        continue  # Skip empty tiles
+
+                    # Create new feature object
+                    tile_feature = {
+                        'type': 'Feature',
+                        'properties': {
+                            'sub_div_id': str(sub_div_id),
+                            'sub_div_center': center
+                        },
+                        'geometry': geometry
+                    }
+                    # Append new feature object to list, and increment sub_div_id for next iteration
+                    features.append(tile_feature)
+                    sub_div_id += 1
             break  # NB test break
 
-        # List to store new GeoJSON feature objects
-        features = []
-
-        # Tracker variables
-        sub_div_id = 0
-
-        # Create subdivisions for each map tile
-        for tile in divided_map:
-            # Calculate tile center based on bounds, and round down to two decimals
-            min_x, min_y, max_x, max_y = tile.bounds
-            center = round(max_y - (max_y - min_y), 6), round(max_x - (max_x - min_x), 6)
-
-            rounded_coordinates = []
-            if isinstance(tile, Polygon):
-                for coords in tile.exterior.coords:
-                    rounded_coords = (round(coords[0], 4), round(coords[1], 4))
-                    rounded_coordinates.append(rounded_coords)
-
-            rounded_tile = Polygon(rounded_coordinates)
-            geometry = rounded_tile.__geo_interface__
-
-            if not geometry['coordinates']:
-                continue  # Skip empty tiles
-
-            # Create new feature object
-            tile_feature = {
-                'type': 'Feature',
-                'properties': {
-                    'sub_div_id': str(sub_div_id),
-                    'sub_div_center': center
-                },
-                'geometry': geometry
-            }
-            # Append new feature object to list, and increment sub_div_id for next iteration
-            features.append(tile_feature)
-            sub_div_id += 1
-
         # Create new GeoJSON object containing all the new feature objects
         feature_collection = {
             'type': 'FeatureCollection',
@@ -170,7 +163,6 @@ def create_grid(poly: Polygon, cell_width: float, cell_height: float):
     return grid_lines
 
 
-
 def write_json_to_file(lake_name: str, map_data: dict):
     """
     Writes a divided map to a JSON file and updates all_lake_names.json
@@ -200,7 +192,6 @@ def write_json_to_file(lake_name: str, map_data: dict):
             json.dump(data, f, ensure_ascii=False, indent=2)
 
 
-
 # Plotting the map can take a considerable amount of time, especially when creating maps with many
 # subdivisions. Removing calls to plot_map will speed up the process, but it is recommended to plot the map
 # after each division to ensure that the map was divided as intended.