GRASS execution POC

opendm/concurrency.py

@@ -9,3 +9,10 @@ def get_max_memory(minimum = 5, use_at_most = 0.5):
     """
     return max(minimum, (100 - virtual_memory().percent) * use_at_most)
 
+def get_max_memory_mb(minimum = 100, use_at_most = 0.5):
+    """
+    :param minimum minimum value to return (return value will never be lower than this)
+    :param use_at_most use at most this fraction of the available memory. 0.5 = use at most 50% of available memory
+    :return value of memory to use in megabytes.
+    """
+    return max(minimum, (virtual_memory().available / 1024 / 1024) * use_at_most)

opendm/grass/generate_cutlines.grass

@@ -0,0 +1,37 @@
+# orthophoto_files: comma-separated GeoTIFF file paths
+# 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 contours file. Otherwise it prints "error"
+
+# Split string using ',' separator
+IFS=',' read -ra DST <<< "${orthophoto_files}"
+ORTHOPHOTO_FILES=("$${DST[@]}")
+
+i=0
+existing_cutlines=""
+for orthophoto_file in "$${ORTHOPHOTO_FILES[@]}"; do
+
+    # Import orthophoto
+    r.external input=$$orthophoto_file output=ortho$$i --overwrite
+
+    # Generate cutlines
+    i.cutlines --overwrite input=ortho$$i output=cutline$$i number_lines=4 edge_detection=zc processes=${max_concurrency} memory=${memory}
+    
+    # Prepend cutline to list of cutlines
+    existing_cutlines="cutline$$i,$$existing_cutlines"
+
+    # Next
+    i=$$[i+1]
+done
+
+last_cutline="cutline$$[i-1]"
+
+# Export
+v.out.ogr input="$$last_cutline" output="cutline.gpkg" format=GPKG
+
+if [ -e "cutline.gpkg" ]; then
+    echo "$$(pwd)/cutline.gpkg"
+else
+    echo "error"
+fi

opendm/grass_engine.py

@@ -89,8 +89,10 @@ class GrassContext:
 
         # Execute it
         log.ODM_INFO("Executing grass script from {}: {} -c {} location --exec sh 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("scripts/grass_addons"))
         p = subprocess.Popen([self.grass_binary, '-c', self.location, 'location', '--exec', 'sh', 'script.sh'],
-                             cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
+                             cwd=self.get_cwd(), stdout=subprocess.PIPE, stderr=subprocess.PIPE, env=env)
         out, err = p.communicate()
 
         out = out.decode('utf-8').strip()

scripts/grass_addons/i.cutlines.py

@@ -0,0 +1,640 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE:	i.cutlines
+# AUTHOR(S):	Moritz Lennert, with help of Stefanos Georganos
+#
+# 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 * (ewrange / nsrange))
+    else:
+        vnumber_lines = number_lines
+        hnumber_lines = int(number_lines * (nsrange / ewrange))
+
+    # 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()

scripts/splitmerge.py

@@ -7,6 +7,7 @@ from opendm import io
 from opendm import system
 from opendm.dem import pdal
 from opensfm.large import metadataset
+from opendm import concurrency
 from pipes import quote
 
 class ODMSplitStage(types.ODM_Stage):
@@ -84,7 +85,7 @@ class ODMSplitStage(types.ODM_Stage):
                     submodel_name = os.path.basename(os.path.abspath(sp_octx.path("..")))
 
                     # Aligned reconstruction is in reconstruction.aligned.json
-                    # We need to replace reconstruction.json with it
+                    # We need to rename it to reconstruction.json
 
                     aligned_recon = sp_octx.path('reconstruction.aligned.json')
                     main_recon = sp_octx.path('reconstruction.json')
@@ -120,17 +121,34 @@ class ODMSplitStage(types.ODM_Stage):
 
 class ODMMergeStage(types.ODM_Stage):
     def process(self, args, outputs):
-        from opendm import grass_engine
+        from opendm.grass_engine import grass 
 
         tree = outputs['tree']
         reconstruction = outputs['reconstruction']
 
         if outputs['large']:
             # Merge point clouds
-            all_point_clouds = get_submodel_paths(tree.submodels_path, "odm_georeferencing", "odm_georeferenced_model.laz")
-            pdal.merge_point_clouds(all_point_clouds, tree.odm_georeferencing_model_laz, args.verbose)
+            # all_point_clouds = get_submodel_paths(tree.submodels_path, "odm_georeferencing", "odm_georeferenced_model.laz")
+            # pdal.merge_point_clouds(all_point_clouds, tree.odm_georeferencing_model_laz, args.verbose)
 
-            # Merge orthophoto
+            # Merge orthophotos
+            all_orthophotos = get_submodel_paths(tree.submodels_path, "odm_orthophoto", "odm_orthophoto.tif")
+            if len(all_orthophotos) > 1:
+                gctx = grass.create_context({'auto_cleanup' : False})
+
+                gctx.add_param('orthophoto_files', ",".join(map(quote, all_orthophotos)))
+                gctx.add_param('max_concurrency', args.max_concurrency)
+                gctx.add_param('memory', concurrency.get_max_memory_mb(300))
+                gctx.set_location(all_orthophotos[0])
+
+                cutline_file = gctx.execute(os.path.join("opendm", "grass", "generate_cutlines.grass"))
+            
+            elif len(all_orthophotos) == 1:
+                # Simply copy
+                log.ODM_WARNING("A single orthophoto was found between all submodels.")
+                shutil.copyfile(all_orthophotos[0], tree.odm_orthophoto_tif)
+            else:
+                log.ODM_WARNING("No orthophotos were found in any of the submodels. No orthophoto will be generated.")
 
             # TODO: crop ortho if necessary