diff --git a/app/lib/pages/widgets/cloropleth_map.dart b/app/lib/pages/widgets/cloropleth_map.dart
index 6cd8d968bcef384623060fc22276515ce62c43a7..9c59c4fd33b7759eaf2883783c657c982d3aa1e4 100644
--- a/app/lib/pages/widgets/cloropleth_map.dart
+++ b/app/lib/pages/widgets/cloropleth_map.dart
@@ -59,82 +59,14 @@ class _ChoroplethMapState extends State<ChoroplethMap> {
                 primaryValueMapper: (int index) => iceThicknessList[index].subDivID,
                 shapeColorValueMapper: (int index) => iceThicknessList[index].thickness,
               ),
-              color: Colors.lightBlueAccent, // Map color
+              color: Colors.blue.shade400, // Map color
               zoomPanBehavior: _zoomPanBehavior,
-              strokeColor: Colors.orange,
+              strokeColor: Colors.black,
             ),
           ],
         ),
-        CustomPaint( // Render polygons clipped on map
-          size: Size.infinite,
-          painter: OverlayPainter(mapData: widget.relation),
-        ),
       ],
     );
   }
 }
 
-
-/// OverlayPainter is a class extending CustomPainter to
-class OverlayPainter extends CustomPainter {
-  final Uint8List mapData;
-
-  const OverlayPainter({
-    required this.mapData, // Pass map shape
-  });
-
-  @override
-  void paint(Canvas canvas, Size size) {
-    // NB: test polygon
-    Path testPol = Path();
-    testPol.moveTo(60.8015, 10.5908);
-    testPol.lineTo(60.5146, 10.8270);
-    testPol.lineTo(60.7115, 11.3242);
-    testPol.close();
-
-    // Decode map data
-    var mapShape = json.decode(utf8.decode(mapData));
-
-    // Extract polygons and render them
-    for (var feature in mapShape['features']) {
-      var geometry = feature['geometry'];
-      var coordinates = geometry['coordinates'];
-
-      for (var ring in coordinates) {
-        var points = <Offset>[];
-        for (var point in ring) {
-          points.add(Offset(point[0], point[1]));
-        }
-        var path = Path()..addPolygon(points, true);
-
-        // Render each map shape
-        canvas.drawPath(path, Paint()..color = Colors.lightBlueAccent);
-      }
-    }
-
-    // Now clip the rendering to the map shape
-    var mapBoundsPath = Path();
-    for (var feature in mapShape['features']) {
-      var geometry = feature['geometry'];
-      var coordinates = geometry['coordinates'];
-
-      for (var ring in coordinates) {
-        var points = <Offset>[];
-        for (var point in ring) {
-          points.add(Offset(point[0], point[1]));
-        }
-        mapBoundsPath.addPolygon(points, true);
-      }
-    }
-    canvas.clipPath(mapBoundsPath);
-
-    // Render test shape over map, not clipped to map
-    canvas.drawPath(testPol, Paint()..color = Colors.pink.shade900);
-  }
-
-
-  @override
-  bool shouldRepaint(covariant CustomPainter oldDelegate) {
-    return false;
-  }
-}
diff --git a/server/map/__pycache__/get_relation.cpython-311.pyc b/server/map/__pycache__/get_relation.cpython-311.pyc
index 1a49b7006d65c76c71c4ebaa80c673617677d1bd..f8284994a9998d150a327d0e5a06b482735e9e22 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 6c7db6417b38e9e7429aa5358192b6e6edd95eb5..360d24515e07077387aa623ec34d298b0de22f6a 100644
--- a/server/map/get_relation.py
+++ b/server/map/get_relation.py
@@ -1,7 +1,10 @@
 import geopandas as gpd
-from shapely.geometry import Polygon, Point, LineString
+from shapely.geometry import Polygon, Point, LineString, MultiLineString, MultiPolygon
 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):
@@ -19,14 +22,32 @@ def get_relation(self, body_of_water: str):
         print("Failed to convert to polygons")
         return
 
+    ########################################################################################
     print("Performing div call")
-    tiles = divide_relation(polygons, 0.01)  # NB: any size under 0.02 will require patience
-    print("Executed 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)
 
-    # NB: temp test plot of map
-    tiles.plot(color='blue', edgecolor='orange', linewidth=0.5)
+    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)
+
+    tiles = gpd.GeoDataFrame(geometry=[new_polygon])
+
+    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()
 
@@ -39,6 +60,7 @@ 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
@@ -75,12 +97,10 @@ def divide_relation(polygons, cell_size):
     return tiles_gdf
 
 
-def divide_relation_2(polygons, cell_size):
-    polygon = Point(0, 0).buffer(2).difference(Point(0, 0).buffer(1))
-    line1 = LineString([(0, 0), (3, 3)])
-    line2 = LineString([(0, 0), (3, -3)])
-
-    line1_pol = line1.buffer(1e-3)
-    line2_pol = line2.buffer(1e-3)
+######################### Approach 2 #################################
 
-    new_polygon = polygon.difference(line1_pol).difference(line2_pol)
+def cut_polygon_by_line(polygon, line):
+    merged = linemerge([polygon.boundary, line])
+    borders = unary_union(merged)
+    polygons = polygonize(borders)
+    return list(polygons)