Geografic-Information-System/scripts/r.buffer.lowmem/r.buffer.lowmem.py

#!/usr/bin/env python3
#
############################################################################
#
# MODULE:    r.buffer.lowmem
# AUTHOR(S): Glynn Clements
# PURPOSE:   Low-memory replacement for r.buffer using r.grow.distance
#
# COPYRIGHT: (C) 2008, 2010 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
# % label: Creates a raster map showing buffer zones surrounding cells that contain non-NULL category values.
# % description: This is the low-memory alternative to the classic r.buffer module.
# % keyword: raster
# % keyword: buffer
# %end
# %flag
# % key: z
# % description: Ignore zero (0) data cells instead of NULL cells
# %end
# %option G_OPT_R_INPUT
# %end
# %option G_OPT_R_OUTPUT
# %end
# %option
# % key: distances
# % type: double
# % required: yes
# % multiple: yes
# % description: Distance zone(s)
# %end
# %option G_OPT_M_UNITS
# % options: meters,kilometers,feet,miles,nautmiles
# % description: Units of distance
# % answer: meters
# %end

import os
import atexit
import grass.script as grass
from grass.script.utils import encode


scales = {
    "meters": 1.0,
    "kilometers": 1000.0,
    "feet": 0.3048,
    "miles": 1609.344,
    "nautmiles": 1852.0,
}

# what to do in case of user break:


def cleanup():
    if grass.find_file(temp_src)["file"]:
        grass.run_command(
            "g.remove", quiet=True, flags="fb", type="raster", name=temp_src
        )
    if grass.find_file(temp_dist)["file"]:
        grass.run_command(
            "g.remove", quiet=True, flags="fb", type="raster", name=temp_dist
        )


def main():
    global temp_dist, temp_src

    input = options["input"]
    output = options["output"]
    distances = options["distances"]
    units = options["units"]
    zero = flags["z"]

    tmp = str(os.getpid())
    temp_dist = "r.buffer.tmp.%s.dist" % tmp
    temp_src = "r.buffer.tmp.%s.src" % tmp

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

    scale = scales[units]

    distances = distances.split(",")
    distances1 = [scale * float(d) for d in distances]
    distances2 = [d * d for d in distances1]

    s = grass.read_command("g.proj", flags="j")
    kv = grass.parse_key_val(s)
    if kv["+proj"] == "longlat":
        metric = "geodesic"
    else:
        metric = "squared"

    grass.run_command(
        "r.grow.distance", input=input, metric=metric, distance=temp_dist, flags="m"
    )

    if zero:
        exp = "$temp_src = if($input == 0,null(),1)"
    else:
        exp = "$temp_src = if(isnull($input),null(),1)"

    grass.message(_("Extracting buffers (1/2)..."))
    grass.mapcalc(exp, temp_src=temp_src, input=input)

    exp = "$output = if(!isnull($input),$input,%s)"
    if metric == "squared":
        for n, dist2 in enumerate(distances2):
            exp %= "if($dist <= %f,%d,%%s)" % (dist2, n + 2)
    else:
        for n, dist2 in enumerate(distances1):
            exp %= "if($dist <= %f,%d,%%s)" % (dist2, n + 2)
    exp %= "null()"

    grass.message(_("Extracting buffers (2/2)..."))
    grass.mapcalc(exp, output=output, input=temp_src, dist=temp_dist)

    p = grass.feed_command("r.category", map=output, separator=":", rules="-")
    msg = "1:distances calculated from these locations\n"
    p.stdin.write(encode(msg))
    d0 = "0"
    for n, d in enumerate(distances):
        msg = "%d:%s-%s %s\n" % (n + 2, d0, d, units)
        p.stdin.write(encode(msg))
        d0 = d
    p.stdin.close()
    p.wait()

    grass.run_command("r.colors", map=output, color="rainbow")

    # write cmd history:
    grass.raster_history(output)


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