Geografic-Information-System/scripts/r.grow/r.grow.py

#!/usr/bin/env python3
#
############################################################################
#
# MODULE:    r.grow
# AUTHOR(S): Glynn Clements
# PURPOSE:   Fast replacement for r.grow using r.grow.distance
#
# COPYRIGHT: (C) 2008 by Glynn Clements
#
#   This program is free software under the GNU General Public
#   License (>=v2). Read the file COPYING that comes with GRASS
#   for details.
#
#############################################################################

# %Module
# % description: Generates a raster map layer with contiguous areas grown by one cell.
# % keyword: raster
# % keyword: distance
# % keyword: proximity
# %end
# %flag
# % key: m
# % description: Radius is in map units rather than cells
# %end
# %option G_OPT_R_INPUT
# %end
# %option G_OPT_R_OUTPUT
# %end
# %option
# % key: radius
# % type: double
# % required: no
# % multiple: no
# % description: Radius of buffer in raster cells
# % answer: 1.01
# %end
# %option
# % key: metric
# % type: string
# % required: no
# % multiple: no
# % options: euclidean,maximum,manhattan
# % description: Metric
# % answer: euclidean
# %end
# %option
# % key: old
# % type: integer
# % required: no
# % multiple: no
# % description: Value to write for input cells which are non-NULL (-1 => NULL)
# %end
# %option
# % key: new
# % type: integer
# % required: no
# % multiple: no
# % description: Value to write for "grown" cells
# %end

import os
import atexit
import math

import grass.script as grass
from grass.exceptions import CalledModuleError


# what to do in case of user break:
def cleanup():
    for map in [temp_dist, temp_val]:
        if map:
            grass.run_command('g.remove', flags='fb', quiet=True,
                              type='rast', name=map)


def main():
    global temp_dist, temp_val

    input = options['input']
    radius = float(options['radius'])
    metric = options['metric']
    old = options['old']
    new = options['new']
    mapunits = flags['m']

    tmp = str(os.getpid())

    temp_dist = "r.grow.tmp.%s.dist" % tmp

    shrink = False
    if radius < 0.0:
        shrink = True
        radius = -radius

    if new == '' and not shrink:
        temp_val = "r.grow.tmp.%s.val" % tmp
        new = '"%s"' % temp_val
    else:
        temp_val = None

    if old == '':
        old = '"%s"' % input

    if not mapunits:
        kv = grass.region()
        scale = math.sqrt(float(kv['nsres']) * float(kv['ewres']))
        radius *= scale

    if metric == 'euclidean':
        metric = 'squared'
        radius = radius * radius

    # check if input file exists
    if not grass.find_file(input)['file']:
        grass.fatal(_("Raster map <%s> not found") % input)

    # Workaround for r.mapcalc bug #3475
    # Mapcalc will fail if output is a fully qualified map name
    out_name = options['output'].split('@')
    if len(out_name) == 2:
        if out_name[1] != grass.gisenv()['MAPSET']:
            grass.fatal(_("Output can be written only to the current mapset"))
        output = out_name[0]
    else:
        output = out_name[0]

    if not shrink:
        try:
            grass.run_command('r.grow.distance', input=input, metric=metric,
                              distance=temp_dist, value=temp_val)
        except CalledModuleError:
            grass.fatal(_("Growing failed. Removing temporary maps."))

        grass.mapcalc(
            '$output = if(!isnull("$input"),$old,if($dist < $radius,$new,null()))',
            output=output, input=input, radius=radius,
            old=old, new=new, dist=temp_dist)
    else:
        # shrink
        try:
            grass.run_command('r.grow.distance', input=input, metric=metric,
                              distance=temp_dist, value=temp_val, flags='n')
        except CalledModuleError:
            grass.fatal(_("Shrinking failed. Removing temporary maps."))

        grass.mapcalc(
            "$output = if(isnull($dist), $old, if($dist < $radius,null(),$old))",
            output=output, radius=radius, old=old, dist=temp_dist)

    grass.run_command('r.colors', map=output, raster=input)

    # write cmd history:
    grass.raster_history(output)


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