[GRASS-SVN] r57552 - in grass-addons/grass7/raster: . r.shaded.pca
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Aug 29 18:29:30 PDT 2013
Author: wenzeslaus
Date: 2013-08-29 18:29:30 -0700 (Thu, 29 Aug 2013)
New Revision: 57552
Added:
grass-addons/grass7/raster/r.shaded.pca/
grass-addons/grass7/raster/r.shaded.pca/Makefile
grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.html
grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.py
Log:
r.shaded.pca: initial version of module for shade-pca-composition
Added: grass-addons/grass7/raster/r.shaded.pca/Makefile
===================================================================
--- grass-addons/grass7/raster/r.shaded.pca/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.shaded.pca/Makefile 2013-08-30 01:29:30 UTC (rev 57552)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.shaded.pca
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.shaded.pca/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.html
===================================================================
--- grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.html (rev 0)
+++ grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.html 2013-08-30 01:29:30 UTC (rev 57552)
@@ -0,0 +1,71 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.shaded.pca</em> is a tool for the generation
+of RGB composite of the three main components
+of PCA created from different hill shades
+(created by <a href="r.shaded.relief.html">r.shaded.relief</a>).
+
+<h3>Input parameters explanation</h3>
+
+Input parameters are the same as for
+<a href="r.shaded.relief.html">r.shaded.relief</a> module except for
+an <i>azimuth</i> parameter which is replaced by <i>nazimuths</i> parameter
+(we need to specify number of different azimuths rather than one)
+and for an <i>nprocs</i> parameter which adds the possibility to run the
+shades creation (<a href="r.shaded.relief.html">r.shaded.relief</a>)
+in parallel. However, the speed of <a href="i.pca.html">r.i.pca</a>
+limits the overall speed of this module.
+
+In order to provide simple interface, it is not possible to customize
+principal component analyses which uses the default settings of the
+<a href="i.pca.html">r.i.pca</a> module.
+
+<h3>Output parameters explanation</h3>
+
+The the standard output map is an RGB composition of first three principal
+components where components are assigned to red, green and blue colors
+in this order.
+
+If you want to create your own RGB composition, HIS composition or do another
+analyses you can specify the <i>pca_shades_basename</i> parameter. If this
+parameter is specified, the module outputs the PCA maps as created during
+the process by <a href="i.pca.html">r.i.pca</a>.
+
+Moreover, if you would like to add one of the shades to your composition,
+you can specify the <i>shades_basename</i> parameter then the module will
+output also the hill shade maps as created during
+the process by <a href="r.shaded.relief.html">r.shaded.relief</a>.
+One of the shades can be used to subtract the intensity channel in HIS
+composition or just as an overlay in your visualization tool.
+
+<h2>EXAMPLE</h2>
+
+<div class="code"><pre>
+# basic example with changed vertical exaggeration
+r.shaded.pca input=elevation output=elevation_pca_shaded zmult=100
+
+# example of more complicated settings
+# including output shades and principal component maps
+r.shaded.pca input=elevation output=elevation_pca_shaded \
+ zmult=100 altitude=15 nazimuths=16 nprocs=4 \
+ shades_basename=elevation_pca_shaded_shades pca_shades_basename=elevation_pca_shaded_pcs
+</pre></div>
+
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.shaded.relief.html">r.shaded.relief</a>,
+<a href="i.pca.html">i.pca</a>
+</em>
+
+<h2>REFERENCES</h2>
+
+Devereux, B. J., Amable, G. S., & Crow, P. P. (2008). Visualisation of LiDAR terrain models for archaeological feature detection. Antiquity, 82(316), 470-479.
+
+<h2>AUTHOR</h2>
+
+Vaclav Petras<br>
+
+<p>
+<i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.py
===================================================================
--- grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.py (rev 0)
+++ grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.py 2013-08-30 01:29:30 UTC (rev 57552)
@@ -0,0 +1,262 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: r.shaded.pca
+# AUTHOR(S): Vaclav Petras
+# PURPOSE: Creates RGB composition from PCA of hill shades
+# COPYRIGHT: (C) 2013 by 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
+#% description: Creates shades from various directions and combines then into RGB composition.
+#% keywords: raster
+#% keywords: elevation
+#% keywords: terrain
+#% keywords: 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: zmult
+#% 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', rast=rast, quiet=True)
+ for pattern in MREMOVE:
+ grass.run_command('g.mremove', rast='%s*' % pattern,
+ flags='f', 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.shaded.relief', input=elevation_input,
+ output=shades_basename + suffix,
+ azimuth=azimuth, zmult=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['zmult'])
+ 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.shade.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.shade.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_prefix=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)
+
+ if pca_basename_user:
+ set_color_table(pca_maps, map_=shade_maps[0])
+
+
+if __name__ == "__main__":
+ atexit.register(cleanup)
+ main()
Property changes on: grass-addons/grass7/raster/r.shaded.pca/r.shaded.pca.py
___________________________________________________________________
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list