i.cutlines.py

#!/usr/bin/env python3

############################################################################
#
# 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 gs
from grass.pygrass.modules.grid.grid import GridModule
from grass.pygrass.vector import VectorTopo
from grass.pygrass.vector import geometry as geom


def cleanup():
    gs.message(_("Erasing temporary files..."))
    for temp_map, maptype in temp_maps:
        if gs.find_file(temp_map, element=maptype)["name"]:
            gs.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
        gs.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:
            gs.run_command(
                "v.patch",
                input_=existing_cutlines,
                output=existingcutlinesmap,
                quiet=True,
                overwrite=True,
            )
            existing_cutlines = existingcutlinesmap

        gs.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])

    gs.message(
        _("Creating edge map using <%s> edgedetection algorithm")
        % edge_detection_algorithm
    )
    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:
            gs.run_command("i.zc", **kwargs)

    elif edge_detection_algorithm == "canny":
        if not gs.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'")
            gs.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,
                flags="n",
                **kwargs,
            )
            grd.run()
        else:
            gs.run_command("i.edge", flags="n", **kwargs)

    else:
        gs.fatal("Only zero-crossing and Canny available as edge detection algorithms.")

    region = gs.region()
    gs.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 = max(int(number_lines * (ewrange / nsrange)), 1)
    else:
        vnumber_lines = number_lines
        hnumber_lines = max(int(number_lines * (nsrange / ewrange)), 1)

    # 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"]

    gs.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)
    gs.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"]

    gs.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)
    gs.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 = gs.start_command(
            "r.cost",
            input_=vbasemap,
            startcoordinates=vstartpoints,
            stopcoordinates=vstoppoints,
            output=vcumcost,
            outdir=vdir,
            memory=pmemory,
            quiet=True,
            overwrite=True,
        )

        rch = gs.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:
        gs.run_command(
            "r.cost",
            input_=vbasemap,
            startcoordinates=vstartpoints,
            stopcoordinates=vstoppoints,
            output=vcumcost,
            outdir=vdir,
            memory=memory,
            quiet=True,
            overwrite=True,
        )

        gs.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 = gs.start_command(
            "r.drain",
            input_=hcumcost,
            direction=hdir,
            startcoordinates=hstoppoints,
            output=hlines,
            flags="d",
            quiet=True,
            overwrite=True,
        )

        rdv = gs.start_command(
            "r.drain",
            input_=vcumcost,
            direction=vdir,
            startcoordinates=vstoppoints,
            output=vlines,
            flags="d",
            quiet=True,
            overwrite=True,
        )

        rdh.wait()
        rdv.wait()

    else:
        gs.run_command(
            "r.drain",
            input_=hcumcost,
            direction=hdir,
            startcoordinates=hstoppoints,
            output=hlines,
            flags="d",
            quiet=True,
            overwrite=True,
        )

        gs.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])
    gs.run_command(
        "r.patch",
        input_=[hlines, vlines],
        output=temp_raster_tile_borders,
        quiet=True,
        overwrite=True,
    )

    gs.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
    gs.use_temp_region()
    gs.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])
    gs.run_command(
        "v.in.region", output=region_map, type_="line", quiet=True, overwrite=True
    )

    gs.del_temp_region()

    gs.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])
    gs.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])
    gs.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])

    gs.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
    gs.run_command(
        "v.category",
        input_=temp_vector_polygons1,
        op="del",
        cat="-1",
        output=temp_vector_polygons2,
        quiet=True,
        overwrite=True,
    )

    # Transform lines to boundaries
    gs.run_command(
        "v.type",
        input_=temp_vector_polygons2,
        from_type="line",
        to_type="boundary",
        output=temp_vector_polygons3,
        quiet=True,
        overwrite=True,
    )

    # Add centroids
    gs.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:
        gs.run_command(
            "v.clean",
            input_=temp_vector_polygons4,
            tool=["rmdangle", "rmarea"],
            threshold=[-1, min_tile_size],
            output=tiles,
            quiet=True,
            overwrite=True,
        )
    else:
        gs.run_command(
            "g.copy", vect=[temp_vector_polygons4, tiles], quiet=True, overwrite=True
        )

    gs.vector_history(tiles)


if __name__ == "__main__":
    options, flags = gs.parser()
    atexit.register(cleanup)
    main()