diff --git a/VERSION b/VERSION
index 437459cd..fad066f8 100644
--- a/VERSION
+++ b/VERSION
@@ -1 +1 @@
-2.5.0
+2.5.0
\ No newline at end of file
diff --git a/opendm/cutline.py b/opendm/cutline.py
index 413b0846..18be8587 100644
--- a/opendm/cutline.py
+++ b/opendm/cutline.py
@@ -1,27 +1,46 @@
 import os
 import shutil
+import rasterio
+import fiona
+import numpy as np
+import math
 from opendm import log
 from opendm import io
 from opendm import concurrency 
 from opendm import get_image_size
 from opendm import system
-import math
 
-def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrency=1, tmpdir=None, scale=1):
+from skimage.feature import canny
+from skimage.draw import line
+from skimage.graph import route_through_array
+from shapely.geometry import LineString, mapping, shape
+from shapely.ops import polygonize, unary_union
+
+def write_raster(data, file):
+    profile = {
+        'driver': 'GTiff',
+        'width': data.shape[1],
+        'height': data.shape[0],
+        'count': 1,
+        'dtype': 'float32',
+        'transform': None,
+        'nodata': None,
+        'crs': None
+    }
+
+    with rasterio.open(file, 'w', **profile) as wout:
+        wout.write(data, 1)
+
+def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrency=1, scale=1):
     if io.file_exists(orthophoto_file) and io.file_exists(crop_area_file):
-        from opendm.grass_engine import grass
         log.ODM_INFO("Computing cutline")
 
-        if tmpdir and not io.dir_exists(tmpdir):
-            system.mkdir_p(tmpdir)
-
         scale = max(0.0001, min(1, scale))
         scaled_orthophoto = None
-
         if scale < 1:
             log.ODM_INFO("Scaling orthophoto to %s%% to compute cutline" % (scale * 100))
 
-            scaled_orthophoto = os.path.join(tmpdir, os.path.basename(io.related_file_path(orthophoto_file, postfix=".scaled")))
+            scaled_orthophoto = io.related_file_path(orthophoto_file, postfix=".scaled")
             # Scale orthophoto before computing cutline
             system.run("gdal_translate -outsize {}% 0 "
                 "-co NUM_THREADS={} "
@@ -33,36 +52,126 @@ def compute_cutline(orthophoto_file, crop_area_file, destination, max_concurrenc
                 orthophoto_file,
                 scaled_orthophoto
             ))
+            
             orthophoto_file = scaled_orthophoto
-
-        try:
-            ortho_width,ortho_height = get_image_size.get_image_size(orthophoto_file, fallback_on_error=False)
-            log.ODM_INFO("Orthophoto dimensions are %sx%s" % (ortho_width, ortho_height))
-            number_lines = int(max(8, math.ceil(min(ortho_width, ortho_height) / 256.0)))
-        except:
-            log.ODM_INFO("Cannot compute orthophoto dimensions, setting arbitrary number of lines.")
-            number_lines = 32
         
-        log.ODM_INFO("Number of lines: %s" % number_lines)
+        # open raster
+        f =  rasterio.open(orthophoto_file)
+        rast = f.read(1) # First band only
+        height, width = rast.shape
+        number_lines = int(max(8, math.ceil(min(width, height) / 256.0)))
+        line_hor_offset = int(width / number_lines)
+        line_ver_offset = int(height / number_lines)
 
-        gctx = grass.create_context({'auto_cleanup' : False, 'tmpdir': tmpdir})
-        gctx.add_param('orthophoto_file', orthophoto_file)
-        gctx.add_param('crop_area_file', crop_area_file)
-        gctx.add_param('number_lines', number_lines)
-        gctx.add_param('max_concurrency', max_concurrency)
-        gctx.add_param('memory', int(concurrency.get_max_memory_mb(300)))
-        gctx.set_location(orthophoto_file)
+        if line_hor_offset <= 2 or line_ver_offset <= 2:
+            log.ODM_WARNING("Cannot compute cutline, orthophoto is too small (%sx%spx)" % (width, height))
+            return
 
-        cutline_file = gctx.execute(os.path.join("opendm", "grass", "compute_cutline.grass"))
-        if cutline_file != 'error':
-            if io.file_exists(cutline_file):
-                shutil.move(cutline_file, destination)
-                log.ODM_INFO("Generated cutline file: %s --> %s" % (cutline_file, destination))
-                gctx.cleanup()
-                return destination
+        crop_f = fiona.open(crop_area_file, 'r')
+        if len(crop_f) == 0:
+            log.ODM_WARNING("Crop area is empty, cannot compute cutline")
+            return
+        crop_poly = shape(crop_f[1]['geometry'])
+        crop_f.close()
+
+        linestrings = []
+
+        # Compute canny edges on first band
+        edges = canny(rast)
+
+        def compute_linestrings(direction):
+            log.ODM_INFO("Computing %s cutlines" % direction)
+            # Initialize cost map
+            cost_map = np.full((height, width), 1, dtype=np.float32)
+
+            # Write edges to cost map
+            cost_map[edges==True] = 0 # Low cost
+
+            # Write "barrier, floor is lava" costs
+            if direction == 'vertical':
+                lines = [((i, 0), (i, height - 1)) for i in range(line_hor_offset, width - line_hor_offset, line_hor_offset)]
+                points = []
+                pad_x = int(line_hor_offset / 2.0)
+                for i in range(0, len(lines)):
+                    a,b = lines[i]
+                    points.append(((a[0] - pad_x , a[1]), (b[0] - pad_x, b[1])))
+                a,b = lines[-1]
+                points.append(((a[0] + pad_x , a[1]), (b[0] + pad_x, b[1])))
             else:
-                log.ODM_WARNING("Unexpected script result: %s. No cutline file has been generated." % cutline_file)
-        else:
-            log.ODM_WARNING("Could not generate orthophoto cutline. An error occured when running GRASS. No orthophoto will be generated.")
+                lines = [((0, j), (width - 1, j)) for j in range(line_ver_offset, height - line_ver_offset, line_ver_offset)]
+                points = []
+                pad_y = int(line_ver_offset / 2.0)
+                for i in range(0, len(lines)):
+                    a,b = lines[i]
+                    points.append(((a[0] , a[1] - pad_y), (b[0], b[1] - pad_y)))
+                a,b = lines[-1]
+                points.append(((a[0] , a[1] + pad_y), (b[0], b[1] + pad_y)))
+                    
+            for a, b in lines:
+                rr,cc = line(*a, *b)
+                cost_map[cc, rr] = 9999 # Lava
+            
+            # Calculate route
+            for a, b in points:
+                line_coords, cost = route_through_array(cost_map, (a[1], a[0]), (b[1], b[0]), fully_connected=True, geometric=True)
+
+                # Convert to geographic
+                geo_line_coords = [f.xy(*c) for c in line_coords]
+
+                # Simplify
+                ls = LineString(geo_line_coords)
+                linestrings.append(ls.simplify(0.05, preserve_topology=False))
+                
+        compute_linestrings('vertical')
+        compute_linestrings('horizontal')
+
+                    
+        # Generate polygons and keep only those inside the crop area
+        log.ODM_INFO("Generating polygons... this could take a bit.")
+        polygons = []
+        for p in polygonize(unary_union(linestrings)):
+            if crop_poly.contains(p):
+                polygons.append(p)
+
+        # This should never happen
+        if len(polygons) == 0:
+            log.ODM_WARNING("No polygons, cannot compute cutline")
+            return
+        
+        log.ODM_INFO("Merging polygons")
+        cutline_polygons = unary_union(polygons)
+        largest_cutline = cutline_polygons[0]
+        max_area = largest_cutline.area
+        for p in cutline_polygons:
+            if p.area > max_area:
+                max_area = p.area
+                largest_cutline = p
+        
+        log.ODM_INFO("Largest cutline found: %s m^2" % max_area)
+
+        meta = {
+            'crs': {'init': str(f.crs).lower() },
+            'driver': 'GPKG',
+            'schema': {
+                'properties': {},
+                'geometry': 'Polygon'
+            }
+        }
+
+        # Remove previous
+        if os.path.exists(destination):
+            os.remove(destination)
+        
+        with fiona.open(destination, 'w', **meta) as sink:
+            sink.write({
+                'geometry': mapping(largest_cutline),
+                'properties': {}
+            })
+        f.close()
+        log.ODM_INFO("Wrote %s" % destination)
+
+        # Cleanup
+        if scaled_orthophoto is not None and os.path.exists(scaled_orthophoto):
+            os.remove(scaled_orthophoto)
     else:
         log.ODM_WARNING("We've been asked to compute cutline, but either %s or %s is missing. Skipping..." % (orthophoto_file, crop_area_file))
diff --git a/opendm/grass/addons/i.cutlinesmod.py b/opendm/grass/addons/i.cutlinesmod.py
deleted file mode 100755
index bd9f9d46..00000000
--- a/opendm/grass/addons/i.cutlinesmod.py
+++ /dev/null
@@ -1,641 +0,0 @@
-#!/usr/bin/env python3
-
-############################################################################
-#
-# MODULE:	i.cutlinesmod
-# AUTHOR(S):	Moritz Lennert, with help of Stefanos Georganos, modified by
-#               Piero Toffanin
-#
-# PURPOSE:	Create tiles the borders of which do not cut across semantically
-#               meaningful objects
-# COPYRIGHT:	(C) 1997-2018 by the GRASS Development Team
-#
-#		This program is free software under the GNU General Public
-#		License (>=v2). Read the file COPYING that comes with GRASS
-#		for details.
-#############################################################################
-
-#%Module
-#% description: Creates semantically meaningful tile borders
-#% keyword: imagery
-#% keyword: tiling
-#%end
-#
-#%option G_OPT_R_INPUT
-#% description: Raster map to use as input for tiling
-#% required: yes
-#%end
-#
-#%option G_OPT_V_OUTPUT
-#% description: Name of output vector map with cutline polygons
-#%end
-#
-#%option
-#% key: number_lines
-#% type: integer
-#% description: Number of tile border lines in each direction
-#% required: yes
-#%end
-#
-#%option
-#% key: edge_detection
-#% type: string
-#% description: Edge detection algorithm to use
-#% options: zc,canny
-#% answer: zc
-#% required: yes
-#%end
-#
-#%option G_OPT_V_INPUTS
-#% key: existing_cutlines
-#% label: Input vector maps with existing cutlines
-#% required: no
-#%end
-#
-#%option
-#% key: no_edge_friction
-#% type: integer
-#% description: Additional friction for non-edge pixels
-#% required: yes
-#% answer: 5
-#%end
-#
-#%option
-#% key: lane_border_multiplier
-#% type: integer
-#% description: Multiplier for borders of lanes compared to non-edge pixels
-#% required: yes
-#% answer: 10
-#%end
-#
-#%option
-#% key: min_tile_size
-#% type: integer
-#% description: Minimum size of tiles in map units
-#% required: no
-#%end
-#
-#%option
-#% key: zc_threshold
-#% type: double
-#% label: Sensitivity of Gaussian filter (i.zc)
-#% answer: 1
-#% required: no
-#% guisection: Zero-crossing
-#%end
-#
-#%option
-#% key: zc_width
-#% type: integer
-#% label: x-y extent of the Gaussian filter (i.zc)
-#% answer: 9
-#% required: no
-#% guisection: Zero-crossing
-#%end
-#
-#%option
-#% key: canny_low_threshold
-#% type: double
-#% label: Low treshold for edges (i.edge)
-#% answer: 3
-#% required: no
-#% guisection: Canny
-#%end
-#
-#%option
-#% key: canny_high_threshold
-#% type: double
-#% label: High treshold for edges (i.edge)
-#% answer: 10
-#% required: no
-#% guisection: Canny
-#%end
-#
-#%option
-#% key: canny_sigma
-#% type: double
-#% label: Kernel radius (i.edge)
-#% answer: 2
-#% required: no
-#% guisection: Canny
-#%end
-#
-#%option
-#% key: tile_width
-#% type: integer
-#% description: Width of tiles for tiled edge detection (pixels)
-#% required: no
-#% guisection: Parallel processing
-#%end
-#
-#%option
-#% key: tile_height
-#% type: integer
-#% description: Height of tiles for tiled edge detection (pixels)
-#% required: no
-#% guisection: Parallel processing
-#%end
-#
-#%option
-#% key: overlap
-#% type: integer
-#% description: Overlap between tiles for tiled edge detection (pixels)
-#% required: no
-#% answer: 1
-#% guisection: Parallel processing
-#%end
-#
-#%option
-#% key: processes
-#% type: integer
-#% description: Number of parallel processes
-#% answer: 1
-#% required: yes
-#% guisection: Parallel processing
-#%end
-#
-#%option
-#% key: memory
-#% type: integer
-#% description: RAM memory available (in MB)
-#% answer: 300
-#% required: yes
-#%end
-#
-#%rules
-#% collective: tile_width, tile_height, overlap
-#%end
-
-import os
-import atexit
-import grass.script as gscript
-from grass.pygrass.modules.grid.grid import GridModule
-from grass.pygrass.vector import VectorTopo
-from grass.pygrass.vector import geometry as geom
-
-def cleanup():
-    gscript.message(_("Erasing temporary files..."))
-    for temp_map, maptype in temp_maps:
-        if gscript.find_file(temp_map, element=maptype)['name']:
-            gscript.run_command('g.remove', flags='f', type=maptype,
-                              name=temp_map, quiet=True)
-
-
-def listzip(input1, input2):
-    # python3 compatible
-    out = zip(input1, input2)
-    if not isinstance(out, list):
-        out = list(zip(input1, input2))
-    return out
-
-
-def main():
-    inputraster = options['input']
-    number_lines = int(options['number_lines'])
-    edge_detection_algorithm = options['edge_detection']
-    no_edge_friction = int(options['no_edge_friction'])
-    lane_border_multiplier = int(options['lane_border_multiplier'])
-    min_tile_size = None
-    if options['min_tile_size']:
-        min_tile_size = float(options['min_tile_size'])
-    existing_cutlines = None
-    if options['existing_cutlines']:
-        existing_cutlines = options['existing_cutlines'].split(',')
-    tiles = options['output']
-    memory = int(options['memory'])
-    tiled = False
-
-    if options['tile_width']:
-        tiled = True
-        gscript.message(_("Using tiles processing for edge detection"))
-        width = int(options['tile_width'])
-        height = int(options['tile_height'])
-        overlap = int(options['overlap'])
-
-    processes = int(options['processes'])
-
-    global temp_maps
-    temp_maps = []
-    r = 'raster'
-    v = 'vector'
-
-    if existing_cutlines:
-        existingcutlinesmap = 'temp_icutlines_existingcutlinesmap_%i' % os.getpid()
-        if len(existing_cutlines) > 1:
-            gscript.run_command('v.patch',
-                                input_=existing_cutlines,
-                                output=existingcutlinesmap,
-                                quiet=True,
-                                overwrite=True)
-            existing_cutlines=existingcutlinesmap
-
-        gscript.run_command('v.to.rast',
-                            input_=existing_cutlines,
-                            output=existingcutlinesmap,
-                            use='val',
-                            type_='line,boundary',
-                            overwrite=True,
-                            quiet=True)
-
-        temp_maps.append([existingcutlinesmap, r])
-
-    temp_edge_map = "temp_icutlines_edgemap_%d" % os.getpid()
-    temp_maps.append([temp_edge_map, r])
-
-    gscript.message(_("Creating edge map"))
-    if edge_detection_algorithm == 'zc':
-        kwargs = {'input' : inputraster,
-                  'output' : temp_edge_map,
-                  'width_' : int(options['zc_width']),
-                  'threshold' : float(options['zc_threshold']),
-                  'quiet' : True}
-
-        if tiled:
-            grd = GridModule('i.zc',
-                             width=width,
-                             height=height,
-                             overlap=overlap,
-                             processes=processes,
-                             split=False,
-                             **kwargs)
-            grd.run()
-        else:
-            gscript.run_command('i.zc',
-                                **kwargs)
-
-    elif edge_detection_algorithm == 'canny':
-        if not gscript.find_program('i.edge', '--help'):
-                message = _("You need to install the addon i.edge to use ")
-                message += _("the Canny edge detector.\n")
-                message += _(" You can install the addon with 'g.extension i.edge'")
-                gscript.fatal(message)
-
-        kwargs = {'input' : inputraster,
-                  'output' : temp_edge_map,
-                  'low_threshold' : float(options['canny_low_threshold']),
-                  'high_threshold' : float(options['canny_high_threshold']),
-                  'sigma' : float(options['canny_sigma']),
-                  'quiet' : True}
-
-
-        if tiled:
-            grd = GridModule('i.edge',
-                             width=width,
-                             height=height,
-                             overlap=overlap,
-                             processes=processes,
-                             split=False,
-                             **kwargs)
-            grd.run()
-        else:
-            gscript.run_command('i.edge',
-                                **kwargs)
-
-    else:
-        gscript.fatal("Only zero-crossing and Canny available as edge detection algorithms.")
-
-    region = gscript.region()
-    gscript.message(_("Finding cutlines in both directions"))
-
-    nsrange = float(region.n - region.s - region.nsres)
-    ewrange = float(region.e - region.w - region.ewres)
-
-    if nsrange > ewrange:
-        hnumber_lines = number_lines
-        vnumber_lines = int(number_lines * (nsrange / ewrange))
-    else:
-        vnumber_lines = number_lines
-        hnumber_lines = int(number_lines * (ewrange / nsrange))
-
-    # Create the lines in horizonal direction
-    nsstep = float(region.n - region.s - region.nsres) / hnumber_lines
-    hpointsy = [((region.n - i * nsstep) - region.nsres / 2.0) for i in range(0, hnumber_lines+1)]
-    hlanepointsy = [y - nsstep / 2.0 for y in hpointsy]
-    hstartpoints = listzip([region.w + 0.2 * region.ewres] * len(hpointsy), hpointsy)
-    hstoppoints = listzip([region.e - 0.2 * region.ewres] * len(hpointsy), hpointsy)
-    hlanestartpoints = listzip([region.w + 0.2 * region.ewres] * len(hlanepointsy), hlanepointsy)
-    hlanestoppoints = listzip([region.e - 0.2 * region.ewres] * len(hlanepointsy), hlanepointsy)
-
-    hlanemap = 'temp_icutlines_hlanemap_%i' % os.getpid()
-    temp_maps.append([hlanemap, v])
-    temp_maps.append([hlanemap, r])
-
-    os.environ['GRASS_VERBOSE'] = '0'
-    new = VectorTopo(hlanemap)
-    new.open('w')
-    for line in listzip(hlanestartpoints,hlanestoppoints):
-        new.write(geom.Line(line), cat=1)
-    new.close()
-    del os.environ['GRASS_VERBOSE']
-
-    gscript.run_command('v.to.rast',
-                        input_=hlanemap,
-                        output=hlanemap,
-                        use='val',
-                        type_='line',
-                        overwrite=True,
-                        quiet=True)
-
-    hbasemap = 'temp_icutlines_hbasemap_%i' % os.getpid()
-    temp_maps.append([hbasemap, r])
-
-    # Building the cost maps using the following logic
-    # - Any pixel not on an edge, nor on an existing cutline gets a
-    # no_edge_friction cost, or no_edge_friction_cost x 10  if there are
-    # existing cutlines
-    # - Any pixel on an edge gets a cost of 1 if there are no existing cutlines,
-    # and a cost of no_edge_friction if there are
-    # - A lane line gets a very high cost (lane_border_multiplier x cost of no
-    # edge pixel - the latter depending on the existence of cutlines).
-
-    mapcalc_expression = "%s = " % hbasemap
-    mapcalc_expression += "if(isnull(%s), " % hlanemap
-    if existing_cutlines:
-        mapcalc_expression += "if(%s == 0 && isnull(%s), " % (temp_edge_map, existingcutlinesmap)
-        mapcalc_expression += "%i, " % (no_edge_friction * 10)
-        mapcalc_expression += "if(isnull(%s), %s, 1))," % (existingcutlinesmap, no_edge_friction)
-        mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction * 10)
-    else:
-        mapcalc_expression += "if(%s == 0, " % temp_edge_map
-        mapcalc_expression += "%i, " % no_edge_friction
-        mapcalc_expression += "1), "
-        mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction)
-    gscript.run_command('r.mapcalc',
-                        expression=mapcalc_expression,
-                        quiet=True,
-                        overwrite=True)
-
-    hcumcost = 'temp_icutlines_hcumcost_%i' % os.getpid()
-    temp_maps.append([hcumcost, r])
-    hdir = 'temp_icutlines_hdir_%i' % os.getpid()
-    temp_maps.append([hdir, r])
-
-
-    # Create the lines in vertical direction
-    ewstep = float(region.e - region.w - region.ewres) / vnumber_lines
-    vpointsx = [((region.e - i * ewstep) - region.ewres / 2.0) for i in range(0, vnumber_lines+1)]
-    vlanepointsx = [x + ewstep / 2.0 for x in vpointsx]
-    vstartpoints = listzip(vpointsx, [region.n - 0.2 * region.nsres] * len(vpointsx))
-    vstoppoints = listzip(vpointsx, [region.s + 0.2 * region.nsres] * len(vpointsx))
-    vlanestartpoints = listzip(vlanepointsx, [region.n - 0.2 * region.nsres] * len(vlanepointsx))
-    vlanestoppoints = listzip(vlanepointsx, [region.s + 0.2 * region.nsres] * len(vlanepointsx))
-
-
-    vlanemap = 'temp_icutlines_vlanemap_%i' % os.getpid()
-    temp_maps.append([vlanemap, v])
-    temp_maps.append([vlanemap, r])
-
-    os.environ['GRASS_VERBOSE'] = '0'
-    new = VectorTopo(vlanemap)
-    new.open('w')
-    for line in listzip(vlanestartpoints,vlanestoppoints):
-        new.write(geom.Line(line), cat=1)
-    new.close()
-    del os.environ['GRASS_VERBOSE']
-
-    gscript.run_command('v.to.rast',
-                        input_=vlanemap,
-                        output=vlanemap,
-                        use='val',
-                        type_='line',
-                        overwrite=True,
-                        quiet=True)
-
-    vbasemap = 'temp_icutlines_vbasemap_%i' % os.getpid()
-    temp_maps.append([vbasemap, r])
-    mapcalc_expression = "%s = " % vbasemap
-    mapcalc_expression += "if(isnull(%s), " % vlanemap
-    if existing_cutlines:
-        mapcalc_expression += "if(%s == 0 && isnull(%s), " % (temp_edge_map, existingcutlinesmap)
-        mapcalc_expression += "%i, " % (no_edge_friction  * 10)
-        mapcalc_expression += "if(isnull(%s), %s, 1))," % (existingcutlinesmap, no_edge_friction)
-        mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction * 10)
-    else:
-        mapcalc_expression += "if(%s == 0, " % temp_edge_map
-        mapcalc_expression += "%i, " % no_edge_friction
-        mapcalc_expression += "1), "
-        mapcalc_expression += "%i)" % (lane_border_multiplier * no_edge_friction)
-    gscript.run_command('r.mapcalc',
-                        expression=mapcalc_expression,
-                        quiet=True,
-                        overwrite=True)
-
-    vcumcost = 'temp_icutlines_vcumcost_%i' % os.getpid()
-    temp_maps.append([vcumcost, r])
-    vdir = 'temp_icutlines_vdir_%i' % os.getpid()
-    temp_maps.append([vdir, r])
-
-    if processes > 1:
-        pmemory = memory / 2.0
-        rcv = gscript.start_command('r.cost',
-                                    input_=vbasemap,
-                                    startcoordinates=vstartpoints,
-                                    stopcoordinates=vstoppoints,
-                                    output=vcumcost,
-                                    outdir=vdir,
-                                    memory=pmemory,
-                                    quiet=True,
-                                    overwrite=True)
-
-        rch = gscript.start_command('r.cost',
-                                    input_=hbasemap,
-                                    startcoordinates=hstartpoints,
-                                    stopcoordinates=hstoppoints,
-                                    output=hcumcost,
-                                    outdir=hdir,
-                                    memory=pmemory,
-                                    quiet=True,
-                                    overwrite=True)
-        rcv.wait()
-        rch.wait()
-
-    else:
-        gscript.run_command('r.cost',
-                            input_=vbasemap,
-                            startcoordinates=vstartpoints,
-                            stopcoordinates=vstoppoints,
-                            output=vcumcost,
-                            outdir=vdir,
-                            memory=memory,
-                            quiet=True,
-                            overwrite=True)
-
-        gscript.run_command('r.cost',
-                            input_=hbasemap,
-                            startcoordinates=hstartpoints,
-                            stopcoordinates=hstoppoints,
-                            output=hcumcost,
-                            outdir=hdir,
-                            memory=memory,
-                            quiet=True,
-                            overwrite=True)
-
-    hlines = 'temp_icutlines_hlines_%i' % os.getpid()
-    temp_maps.append([hlines, r])
-    vlines = 'temp_icutlines_vlines_%i' % os.getpid()
-    temp_maps.append([vlines, r])
-
-    if processes > 1:
-        rdh = gscript.start_command('r.drain',
-                                    input_=hcumcost,
-                                    direction=hdir,
-                                    startcoordinates=hstoppoints,
-                                    output=hlines,
-                                    flags='d',
-                                    quiet=True,
-                                    overwrite=True)
-
-
-        rdv = gscript.start_command('r.drain',
-                                    input_=vcumcost,
-                                    direction=vdir,
-                                    startcoordinates=vstoppoints,
-                                    output=vlines,
-                                    flags='d',
-                                    quiet=True,
-                                    overwrite=True)
-
-        rdh.wait()
-        rdv.wait()
-
-    else:
-        gscript.run_command('r.drain',
-                            input_=hcumcost,
-                            direction=hdir,
-                            startcoordinates=hstoppoints,
-                            output=hlines,
-                            flags='d',
-                            quiet=True,
-                            overwrite=True)
-
-
-        gscript.run_command('r.drain',
-                            input_=vcumcost,
-                            direction=vdir,
-                            startcoordinates=vstoppoints,
-                            output=vlines,
-                            flags='d',
-                            quiet=True,
-                            overwrite=True)
-
-    # Combine horizonal and vertical lines
-    temp_raster_tile_borders = 'temp_icutlines_raster_tile_borders_%i' % os.getpid()
-    temp_maps.append([temp_raster_tile_borders, r])
-    gscript.run_command('r.patch',
-                        input_=[hlines,vlines],
-                        output=temp_raster_tile_borders,
-                        quiet=True,
-                        overwrite=True)
-
-    gscript.message(_("Creating vector polygons"))
-
-    # Create vector polygons
-
-    # First we need to shrink the region a bit to make sure that all vector
-    # points / lines fall within the raster
-    gscript.use_temp_region()
-    gscript.run_command('g.region',
-                        s=region.s+region.nsres,
-                        e=region.e-region.ewres,
-                        quiet=True)
-
-    region_map = 'temp_icutlines_region_map_%i' % os.getpid()
-    temp_maps.append([region_map, v])
-    temp_maps.append([region_map, r])
-    gscript.run_command('v.in.region',
-                        output=region_map,
-                        type_='line',
-                        quiet=True,
-                        overwrite=True)
-
-    gscript.del_temp_region()
-
-    gscript.run_command('v.to.rast',
-                        input_=region_map,
-                        output=region_map,
-                        use='val',
-                        type_='line',
-                        quiet=True,
-                        overwrite=True)
-
-    temp_raster_polygons = 'temp_icutlines_raster_polygons_%i' % os.getpid()
-    temp_maps.append([temp_raster_polygons, r])
-    gscript.run_command('r.patch',
-                        input_=[temp_raster_tile_borders,region_map],
-                        output=temp_raster_polygons,
-                        quiet=True,
-                        overwrite=True)
-
-    temp_raster_polygons_thin = 'temp_icutlines_raster_polygons_thin_%i' % os.getpid()
-    temp_maps.append([temp_raster_polygons_thin, r])
-    gscript.run_command('r.thin',
-                        input_=temp_raster_polygons,
-                        output=temp_raster_polygons_thin,
-                        quiet=True,
-                        overwrite=True)
-
-    # Create a series of temporary map names as we have to go
-    # through several steps until we reach the final map.
-    temp_vector_polygons1 = 'temp_icutlines_vector_polygons1_%i' % os.getpid()
-    temp_maps.append([temp_vector_polygons1, v])
-    temp_vector_polygons2 = 'temp_icutlines_vector_polygons2_%i' % os.getpid()
-    temp_maps.append([temp_vector_polygons2, v])
-    temp_vector_polygons3 = 'temp_icutlines_vector_polygons3_%i' % os.getpid()
-    temp_maps.append([temp_vector_polygons3, v])
-    temp_vector_polygons4 = 'temp_icutlines_vector_polygons4_%i' % os.getpid()
-    temp_maps.append([temp_vector_polygons4, v])
-
-    gscript.run_command('r.to.vect',
-                        input_=temp_raster_polygons_thin,
-                        output=temp_vector_polygons1,
-                        type_='line',
-                        flags='t',
-                        quiet=True,
-                        overwrite=True)
-
-    # Erase all category values from the lines
-    gscript.run_command('v.category',
-                        input_=temp_vector_polygons1,
-                        op='del',
-                        cat='-1',
-                        output=temp_vector_polygons2,
-                        quiet=True,
-                        overwrite=True)
-
-    # Transform lines to boundaries
-    gscript.run_command('v.type',
-                        input_=temp_vector_polygons2,
-                        from_type='line',
-                        to_type='boundary',
-                        output=temp_vector_polygons3,
-                        quiet=True,
-                        overwrite=True)
-
-    # Add centroids
-    gscript.run_command('v.centroids',
-                        input_=temp_vector_polygons3,
-                        output=temp_vector_polygons4,
-                        quiet=True,
-                        overwrite=True)
-
-    # If a threshold is given erase polygons that are too small
-    if min_tile_size:
-        gscript.run_command('v.clean',
-                            input_=temp_vector_polygons4,
-                            tool='rmarea',
-                            threshold=min_tile_size,
-                            output=tiles,
-                            quiet=True,
-                            overwrite=True)
-    else:
-        gscript.run_command('g.copy',
-                            vect=[temp_vector_polygons4,tiles],
-                            quiet=True,
-                            overwrite=True)
-
-    gscript.vector_history(tiles)
-
-if __name__ == "__main__":
-    options, flags = gscript.parser()
-    atexit.register(cleanup)
-    main()
diff --git a/opendm/grass/compute_cutline.grass b/opendm/grass/compute_cutline.grass
deleted file mode 100644
index 2a062ef2..00000000
--- a/opendm/grass/compute_cutline.grass
+++ /dev/null
@@ -1,38 +0,0 @@
-# orthophoto_file: input GeoTIFF raster file
-# crop_area_file: input vector polygon file delimiting the safe area for processing
-# number_lines: number of cutlines on the smallest side of the orthophoto for computing the final cutline
-# max_concurrency: maximum number of parallel processes to use
-# memory: maximum MB of memory to use
-# ------
-# output: If successful, prints the full path to the cutlines file. Otherwise it prints "error"
-
-# Import orthophoto (green band only)
-r.external band=2 input="${orthophoto_file}" output=ortho --overwrite
-
-# Import crop area
-v.in.ogr input="${crop_area_file}" output=crop_area --overwrite
-
-g.region vector=crop_area
-
-# Generate cutlines
-i.cutlinesmod.py --overwrite input=ortho output=cutline number_lines=${number_lines} edge_detection=zc no_edge_friction=20 lane_border_multiplier=1000000 tile_width=1024 tile_height=1024 overlap=20 processes=${max_concurrency} memory=${memory}
-
-#v.out.ogr input=cutline output="cutline_raw.gpkg" format=GPKG
-
-# Select cutlines that are within crop area
-v.select ainput=cutline binput=crop_area output=result operator=within
-
-# Export
-v.out.ogr input=result output="result.gpkg" format=GPKG --overwrite
-
-# Merge all geometries, select only the largest one (remove islands)
-ogr2ogr -f GPKG -overwrite -explodecollections -dialect SQLite -sql "SELECT ST_Union(geom) FROM result ORDER BY ST_AREA(geom) DESC LIMIT 1" cutline.gpkg result.gpkg
-
-# Add new line output in case the last command didn't.
-echo ""
-
-if [ -e "cutline.gpkg" ]; then
-    echo "$$(pwd)/cutline.gpkg"
-else
-    echo "error"
-fi
diff --git a/opendm/grass_engine.py b/opendm/grass_engine.py
deleted file mode 100644
index 1f9b48e3..00000000
--- a/opendm/grass_engine.py
+++ /dev/null
@@ -1,143 +0,0 @@
-import shutil
-import tempfile
-import subprocess
-import os
-import sys
-import time
-from opendm import log
-from opendm import system
-import locale
-
-from string import Template
-
-class GrassEngine:
-    def __init__(self):
-        self.grass_binary = system.which('grass7') or \
-                            system.which('grass72') or \
-                            system.which('grass74') or \
-                            system.which('grass76') or \
-                            shutil.which('grass78') or \
-                            shutil.which('grass80')
-
-        if self.grass_binary is None:
-            log.ODM_WARNING("Could not find a GRASS 7 executable. GRASS scripts will not work.")
-        else:
-            log.ODM_INFO("Initializing GRASS engine using {}".format(self.grass_binary))
-
-    def create_context(self, serialized_context = {}):
-        if self.grass_binary is None: raise GrassEngineException("GRASS engine is unavailable")
-        return GrassContext(self.grass_binary, **serialized_context)
-
-
-class GrassContext:
-    def __init__(self, grass_binary, tmpdir = None, template_args = {}, location = None, auto_cleanup=True):
-        self.grass_binary = grass_binary
-        if tmpdir is None:
-            tmpdir = tempfile.mkdtemp('_grass_engine')
-        self.tmpdir = tmpdir
-        self.template_args = template_args
-        self.location = location
-        self.auto_cleanup = auto_cleanup
-
-    def get_cwd(self):
-        return self.tmpdir
-
-    def add_file(self, filename, source, use_as_location=False):
-        param = os.path.splitext(filename)[0] # filename without extension
-
-        dst_path = os.path.abspath(os.path.join(self.get_cwd(), filename))
-        with open(dst_path, 'w') as f:
-            f.write(source)
-        self.template_args[param] = dst_path
-
-        if use_as_location:
-            self.set_location(self.template_args[param])
-
-        return dst_path
-
-    def add_param(self, param, value):
-        self.template_args[param] = value
-
-    def set_location(self, location):
-        """
-        :param location: either a "epsg:XXXXX" string or a path to a geospatial file defining the location
-        """
-        if not location.lower().startswith('epsg:'):
-            location = os.path.abspath(location)
-        self.location = location
-
-    def execute(self, script):
-        """
-        :param script: path to .grass script
-        :return: script output
-        """
-        if self.location is None: raise GrassEngineException("Location is not set")
-
-        script = os.path.abspath(script)
-
-        # Create grass script via template substitution
-        try:
-            with open(script) as f:
-                script_content = f.read()
-        except FileNotFoundError:
-            raise GrassEngineException("Script does not exist: {}".format(script))
-
-        tmpl = Template(script_content)
-
-        # Write script to disk
-        if not os.path.exists(self.get_cwd()):
-            os.mkdir(self.get_cwd())
-
-        with open(os.path.join(self.get_cwd(), 'script.sh'), 'w') as f:
-            f.write(tmpl.substitute(self.template_args))
-
-        # Execute it
-        log.ODM_INFO("Executing grass script from {}: {} --tmp-location {} --exec bash script.sh".format(self.get_cwd(), self.grass_binary, self.location))
-        env = os.environ.copy()
-        env["GRASS_ADDON_PATH"] = env.get("GRASS_ADDON_PATH", "") + os.path.abspath(os.path.join("opendm/grass/addons"))
-        env["LC_ALL"] = "C.UTF-8"
-
-        filename = os.path.join(self.get_cwd(), 'output.log')
-        with open(filename, 'wb') as writer, open(filename, 'rb', 1) as reader:
-            p = subprocess.Popen([self.grass_binary, '--tmp-location', self.location, '--exec', 'bash', 'script.sh'],
-                                cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=writer, env=env)
-            
-            while p.poll() is None:
-                sys.stdout.write(reader.read().decode('utf8'))
-                time.sleep(0.5)
-            
-            # Read the remaining
-            sys.stdout.write(reader.read().decode('utf8'))
-
-            out, err = p.communicate()
-            out = out.decode('utf-8').strip()
-
-            if p.returncode == 0:
-                return out
-            else:
-                raise GrassEngineException("Could not execute GRASS script {} from {}: {}".format(script, self.get_cwd(), err))
-
-    def serialize(self):
-        return {
-            'tmpdir': self.tmpdir,
-            'template_args': self.template_args,
-            'location': self.location,
-            'auto_cleanup': self.auto_cleanup
-        }
-
-    def cleanup(self):
-        if os.path.exists(self.get_cwd()):
-            shutil.rmtree(self.get_cwd())
-
-    def __del__(self):
-        if self.auto_cleanup:
-            self.cleanup()
-
-class GrassEngineException(Exception):
-    pass
-
-def cleanup_grass_context(serialized_context):
-    ctx = grass.create_context(serialized_context)
-    ctx.cleanup()
-
-grass = GrassEngine()
\ No newline at end of file
diff --git a/requirements.txt b/requirements.txt
index ed1fb849..25f718bb 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -28,3 +28,4 @@ scikit-image==0.17.2
 scipy==1.5.4
 xmltodict==0.12.0
 fpdf2==2.2.0rc2
+Shapely==1.7.1
\ No newline at end of file
diff --git a/snap/snapcraft.yaml b/snap/snapcraft.yaml
index 65430ee2..07a023c7 100644
--- a/snap/snapcraft.yaml
+++ b/snap/snapcraft.yaml
@@ -39,7 +39,6 @@ parts:
       - gdal-bin
       - gfortran # to build scipy
       - git
-      - grass-core
       - libboost-log-dev
       - libgdal-dev
       - libgeotiff-dev
@@ -57,7 +56,6 @@ parts:
       - swig3.0
     stage-packages:
       - gdal-bin
-      - grass-core
       - libboost-log1.71.0
       - libgdal26
       - libgeotiff5
diff --git a/stages/odm_orthophoto.py b/stages/odm_orthophoto.py
index f9fac5de..a96ac647 100644
--- a/stages/odm_orthophoto.py
+++ b/stages/odm_orthophoto.py
@@ -127,7 +127,6 @@ class ODMOrthoPhotoStage(types.ODM_Stage):
                                     bounds_file_path,
                                     cutline_file,
                                     args.max_concurrency,
-                                    tmpdir=os.path.join(tree.odm_orthophoto, "grass_cutline_tmpdir"),
                                     scale=0.25)
 
                     orthophoto.compute_mask_raster(tree.odm_orthophoto_tif, cutline_file,