#!/usr/bin/env python # ############################################################################## # # MODULE: r.skyview # # AUTHOR(S): Anna Petrasova (kratochanna gmail.com) # Vaclav Petras (colorize addition) # # PURPOSE: Implementation of Sky-View Factor visualization technique # # COPYRIGHT: (C) 2013-2014 by the GRASS Development Team # # This program is free software under the GNU General Public # License (version 2). Read the file COPYING that comes with GRASS # for details. # ############################################################################## #%module #% description: Computes Sky-View Factor visualization technique #% keyword: raster #% keyword: visualization #%end #%option G_OPT_R_INPUT #%end #%option G_OPT_R_OUTPUT #%end #%option #% key: ndir #% description: Number of directions (8 to 32 recommended) #% type: integer #% required: yes #% answer: 16 #% options: 2-360 #%end #%option #% key: maxdistance #% description: The maximum distance to consider when finding the horizon height #% type: double #% required: no #%end #%option #% key: color_source #% type: string #% label: Source raster for colorization #% description: Input and color_input are taken from input and color_input options respectively. The rest is computed using r.slope.aspect #% descriptions: input; use the raster from the input option;color_input;use the raster from the color_input option;slope;compute and use slope;aspect;compute and use aspect;dxy;compute and use second order partial derivative dxy #% multiple: no #% required: no #% options: input, color_input, slope, aspect, dxy #% answer: input #% guisection: Colorize #%end #%option G_OPT_R_INPUT #% key: color_input #% required: no #% description: Custom raster map to be used for colorization #% guisection: Colorize #%end #%option #% key: color_table #% type: string #% label: Color table for colorization raster (preset color table by default) #% description: If empty, the color table of the created raster is used (not used at all for input and color_input) #% multiple: no #% required: no #% options: reds, blues, greens, oranges, sepia, aspectcolr #% guisection: Colorize #%end #%option G_OPT_R_OUTPUT #% key: colorized_output #% required: no #% description: Colorized sky-view factor #% guisection: Colorize #%end #%option #% key: basename #% type: string #% multiple: no #% description: Set the basename for the intermediate maps #%end #%flag #% key: n #% label: Invert color table for colorization raster #% description: Ignored for input and color_input #% guisection: Colorize #%end import sys import os import atexit from grass.exceptions import CalledModuleError import grass.script.core as gcore import grass.script.raster as grast from grass.pygrass.messages import get_msgr # TODO: also used for r.slope.aspect result TMP_NAME = 'tmp_horizon_' + str(os.getpid()) CLEANUP = True def cleanup(): if CLEANUP: gcore.verbose(_("Cleaning temporary maps...")) gcore.run_command('g.remove', flags='f', type='raster', pattern=TMP_NAME + "*", quiet=True) def main(): elev = options['input'] output = options['output'] n_dir = int(options['ndir']) global TMP_NAME, CLEANUP if options['basename']: TMP_NAME = options['basename'] CLEANUP = False colorized_output = options['colorized_output'] colorize_color = options['color_table'] if colorized_output: color_raster_tmp = TMP_NAME + "_color_raster" else: color_raster_tmp = None color_raster_type = options['color_source'] color_input = options['color_input'] if color_raster_type == 'color_input' and not color_input: gcore.fatal(_("Provide raster name in color_input option")) if color_raster_type != 'color_input' and color_input: gcore.fatal(_("The option color_input is not needed" " when not using it as source for color")) # this would be needed only when no value would allowed if not color_raster_type and color_input: color_raster_type = 'color_input' # enable for convenience if color_raster_type == 'aspect' \ and colorize_color \ and colorize_color not in ['default', 'aspectcolr']: gcore.warning(_("Using possibly inappropriate color table <{}>" " for aspect".format(colorize_color))) horizon_step = 360. / n_dir msgr = get_msgr() # checks if there are already some maps old_maps = _get_horizon_maps() if old_maps: if not gcore.overwrite(): CLEANUP = False msgr.fatal(_("You have to first check overwrite flag or remove" " the following maps:\n" "{names}").format(names=','.join(old_maps))) else: msgr.warning(_("The following maps will be overwritten: {names}" ).format(names=','.join(old_maps))) if not gcore.overwrite() and color_raster_tmp: check_map_name(color_raster_tmp) try: params = {} if options['maxdistance']: params['maxdistance'] = options['maxdistance'] gcore.run_command('r.horizon', elevation=elev, step=horizon_step, output=TMP_NAME, flags='d', **params) msgr.message(_("Computing sky view factor ...")) new_maps = _get_horizon_maps() expr = '"{out}" = 1 - (sin("{first}") '.format(first=new_maps[0], out=output) for horizon in new_maps[1:]: expr += '+ sin("{name}") '.format(name=horizon) expr += ") / {n}.".format(n=len(new_maps)) grast.mapcalc(exp=expr) gcore.run_command('r.colors', map=output, color='grey') except CalledModuleError: msgr.fatal(_("r.horizon failed to compute horizon elevation " "angle maps. Please report this problem to developers.")) return 1 if colorized_output: if color_raster_type == 'slope': gcore.run_command('r.slope.aspect', elevation=elev, slope=color_raster_tmp) elif color_raster_type == 'aspect': gcore.run_command('r.slope.aspect', elevation=elev, aspect=color_raster_tmp) elif color_raster_type == 'dxy': gcore.run_command('r.slope.aspect', elevation=elev, dxy=color_raster_tmp) elif color_raster_type == 'color_input': color_raster_tmp = color_input else: color_raster_tmp = elev # don't modify user's color table for inputs if colorize_color \ and color_raster_type not in ['input', 'color_input']: rcolors_flags = '' if flags['n']: rcolors_flags += 'n' gcore.run_command('r.colors', map=color_raster_tmp, color=colorize_color, flags=rcolors_flags) gcore.run_command('r.shade', shade=output, color=color_raster_tmp, output=colorized_output) grast.raster_history(colorized_output) grast.raster_history(output) return 0 def _get_horizon_maps(): return gcore.list_grouped('rast', pattern=TMP_NAME + "*")[gcore.gisenv()['MAPSET']] def check_map_name(name): # cell means any raster in this context # mapset needs to retrieved in very call, ok for here if gcore.find_file(name, element='cell', mapset=gcore.gisenv()['MAPSET'])['file']: gcore.fatal(_("Raster map <%s> already exists. " "Remove the existing map or allow overwrite.") % name) if __name__ == "__main__": options, flags = gcore.parser() atexit.register(cleanup) sys.exit(main())