#!/usr/bin/env python ############################################################################ # # MODULE: r.shaded.pca # AUTHOR(S): Vaclav Petras # PURPOSE: Creates RGB composition from PCA of hill shades # COPYRIGHT: (C) 2013-2014 by Vaclav Petras and 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 #% label: Creates relief shades from various directions and combines them into RGB composition. #% description: The combined shades highlight terrain features which wouldn't be visible using standard shading technique. #% keyword: raster #% keyword: elevation #% keyword: terrain #% keyword: visualization #%end #%option #% type: string #% gisprompt: old,cell,raster #% key: input #% description: Name of the input elevation raster map #% required: yes #%end #%option #% type: string #% gisprompt: new,cell,raster #% key: output #% description: Name for output PCA shaded relief map #% required: yes #%end #%option #% key: altitude #% type: double #% description: Altitude of the sun in degrees above the horizon #% options: 0-90 #% answer: 30 #%end #%option #% key: nazimuths #% type: integer #% description: The number of azimuths (suggested values are 4, 8, 16, 32) #% answer: 8 #%end #%option #% key: zscale #% type: double #% description: Factor for exaggerating relief #% answer: 1 #%end #%option #% key: scale #% type: double #% description: Azimuth of the sun in degrees to the east of north #% answer: 1 #%end #%option #% key: units #% type: string #% description: Elevation units (overrides scale factor) #% options: intl,survey #% descriptions: intl;international feet;survey;survey feet #%end #%option #% key: shades_basename #% type: string #% label: Base name for output shades map #% description: A base of the name of shades maps for all azimuths. An underscore ('_') and a azimuth will be added to the base name. When empty, no maps will be outputted (although they need to be generated). #%end #%option #% key: pca_shades_basename #% type: string #% label: Base name for output PCA shades map #% description: A base of the name of PCA shades maps. An underscore ('_') and a azimuth will be added to the base name. When empty, no maps will be outputted (although they need to be generated). #%end #%option #% key: nprocs #% type: integer #% description: Number of r.shade.relief processes to run in parallel #% options: 1- #% answer: 1 #%end import os import atexit from multiprocessing import Process import grass.script as grass import grass.script.core as core REMOVE = [] MREMOVE = [] def cleanup(): if REMOVE or MREMOVE: core.info(_("Cleaning temporary maps...")) for rast in REMOVE: grass.run_command("g.remove", flags="f", type="raster", name=rast, quiet=True) for pattern in MREMOVE: grass.run_command( "g.remove", flags="f", type="raster", pattern="%s*" % pattern, quiet=True ) def is_grass_7(): if core.version()["version"].split(".")[0] == "7": return True return False def create_tmp_map_name(name): return "{mod}_{pid}_{map_}_tmp".format( mod="r_shaded_pca", pid=os.getpid(), map_=name ) # add latitude map def run_r_shaded_relief( elevation_input, shades_basename, altitude, azimuth, z_exaggeration, scale, units, suffix, ): params = {} if units: params.update({"units": units}) grass.run_command( "r.relief", input=elevation_input, output=shades_basename + suffix, azimuth=azimuth, zscale=z_exaggeration, scale=scale, altitude=altitude, overwrite=core.overwrite(), quiet=True, **params ) def set_color_table(rasters, map_): if is_grass_7(): grass.run_command("r.colors", map=rasters, raster=map_, quiet=True) else: for rast in rasters: grass.run_command("r.colors", map=rast, raster=map_, quiet=True) def check_map_names(basename, mapset, suffixes): for suffix in suffixes: map_ = "%s%s%s" % (basename, "_", suffix) if grass.find_file(map_, element="cell", mapset=mapset)["file"]: grass.fatal( _( "Raster map <%s> already exists. " "Change the base name or allow overwrite." ) % map_ ) def frange(x, y, step): # we want to be close to range behaviour, so no <= while x < y: yield x x += step def main(): options, flags = grass.parser() elevation_input = options["input"] pca_shade_output = options["output"] altitude = float(options["altitude"]) number_of_azimuths = int(options["nazimuths"]) z_exaggeration = float(options["zscale"]) scale = float(options["scale"]) units = options["units"] shades_basename = options["shades_basename"] pca_basename = pca_basename_user = options["pca_shades_basename"] nprocs = int(options["nprocs"]) full_circle = 360 # let's use floats here and leave the consequences to the user smallest_angle = float(full_circle) / number_of_azimuths azimuths = list(frange(0, full_circle, smallest_angle)) if not shades_basename: shades_basename = create_tmp_map_name("shade") MREMOVE.append(shades_basename) if not pca_basename: pca_basename = pca_shade_output + "_pca" pca_maps = [pca_basename + "." + str(i) for i in range(1, number_of_azimuths + 1)] if not pca_basename_user: REMOVE.extend(pca_maps) # here we check all the posible if not grass.overwrite(): check_map_names(shades_basename, grass.gisenv()["MAPSET"], suffixes=azimuths) check_map_names( pca_basename, grass.gisenv()["MAPSET"], suffixes=range(1, number_of_azimuths), ) grass.info(_("Running r.relief in a loop...")) count = 0 # Parallel processing proc_list = [] proc_count = 0 suffixes = [] all_suffixes = [] core.percent(0, number_of_azimuths, 1) for azimuth in azimuths: count += 1 core.percent(count, number_of_azimuths, 10) suffix = "_" + str(azimuth) proc_list.append( Process( target=run_r_shaded_relief, args=( elevation_input, shades_basename, altitude, azimuth, z_exaggeration, scale, units, suffix, ), ) ) proc_list[proc_count].start() proc_count += 1 suffixes.append(suffix) all_suffixes.append(suffix) if ( proc_count == nprocs or proc_count == number_of_azimuths or count == number_of_azimuths ): proc_count = 0 exitcodes = 0 for proc in proc_list: proc.join() exitcodes += proc.exitcode if exitcodes != 0: core.fatal(_("Error during r.relief computation")) # Empty process list proc_list = [] suffixes = [] # FIXME: how percent really works? # core.percent(1, 1, 1) shade_maps = [shades_basename + suf for suf in all_suffixes] grass.info(_("Running r.pca...")) # not quiet=True to get percents grass.run_command( "i.pca", input=shade_maps, output=pca_basename, overwrite=core.overwrite() ) grass.info( _("Creating RGB composite from " "PC1 (red), PC2 (green), PC3 (blue) ...") ) grass.run_command( "r.composite", red=pca_maps[0], green=pca_maps[1], blue=pca_maps[2], output=pca_shade_output, overwrite=core.overwrite(), quiet=True, ) grass.raster_history(pca_shade_output) if pca_basename_user: set_color_table(pca_maps, map_=shade_maps[0]) if __name__ == "__main__": atexit.register(cleanup) main()