#!/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()