[GRASS-SVN] r57849 - in grass-addons/grass7/raster: . r.local.relief
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Sep 26 17:18:12 PDT 2013
Author: wenzeslaus
Date: 2013-09-26 17:18:12 -0700 (Thu, 26 Sep 2013)
New Revision: 57849
Added:
grass-addons/grass7/raster/r.local.relief/
grass-addons/grass7/raster/r.local.relief/Makefile
grass-addons/grass7/raster/r.local.relief/r.local.relief.html
grass-addons/grass7/raster/r.local.relief/r.local.relief.py
Log:
r.local.relief: initial commit with r.fillnulls implementation (co-author: Eric Goddard)
Added: grass-addons/grass7/raster/r.local.relief/Makefile
===================================================================
--- grass-addons/grass7/raster/r.local.relief/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.local.relief/Makefile 2013-09-27 00:18:12 UTC (rev 57849)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM=r.local.relief
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Property changes on: grass-addons/grass7/raster/r.local.relief/Makefile
___________________________________________________________________
Added: svn:mime-type
+ text/x-makefile
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.local.relief/r.local.relief.html
===================================================================
--- grass-addons/grass7/raster/r.local.relief/r.local.relief.html (rev 0)
+++ grass-addons/grass7/raster/r.local.relief/r.local.relief.html 2013-09-27 00:18:12 UTC (rev 57849)
@@ -0,0 +1,62 @@
+<h2>DESCRIPTION</h2>
+
+<p><em>r.local.relief</em> generates a local relief model (LRM) from lidar-derived high-resolution DEMs. Local relief models enhance the visibility of small-scale surface features by removing large-scale landforms from the DEM.</p>
+
+<p>Generating the LRM is accomplished in 7 steps (Hesse 2010:69):</p>
+<ol>
+ <li>Creation of the DEM from the LIDAR data. Buildings, trees and other objects on the earth's surface should be removed.</li>
+ <li>Apply a low pass filter to the DEM. The low pass filter approximates the large-scale landscape forms. The kernel size of the low pass filter determines the scale of features that will be visible in the LRM. A default kernel size of 25 is used.</li>
+ <li>Subtract the low-pass filter result from the DEM to get the local relief.</li>
+ <li>Extract the zero contour lines from the difference map.</li>
+ <li>Extract the point elevations from the DEM along the zero contour lines.</li>
+ <li>Create a purged DEM from the point elevations extracted in the last step. This layer represents the large-scale landscape forms that will be removed to expose the local relief in the final step. v.surf.bspline is used to interpolate the points to a raster.</li>
+ <li>Subtract the purged DEM from the original DEM to get the local relief model.
+</ol>
+Note that these steps are not completely valid for Vaclav's implementation which uses <a href="r.fillnulls.html">r.fillnulls</a> module.
+
+<h2>OUTPUT</h2>
+<p>
+The final local relief model gets the name specified by the <em>output</em> parameter.
+When the <em>-i</em> is specified <em>r.local.relief</em> creates six additional output files one for each step in the LRM creation process. The names consist of the user-specified <em>output</em> parameter and the following suffixes:</p>
+<ul>
+ <li><tt>_smooth_elevation</tt>: The result of running the low pass filter on the DEM.</li>
+ <li><tt>_subtracted_smooth_elevation:s</tt> The result of subtracting the low pass filter map from the DEM.</li>
+ <li><tt>_vector_contours:s</tt> The zero contours extracted from the subtracted smooth elevation map.</li>
+ <li><tt>raster_contours</tt>: The rasterized zero contours.</li>
+ <li><tt>raster_contours_with_values</tt>: The rasterized zero contours with the values from elevation map.</li>
+ <li><tt>_purged_elevation:s</tt> The raster interpolated from the raster map based on zero contours map.</li>
+</ul>
+Note that this is not implemented yet, you must change the variable save_intermediates in the code itself.
+
+<h2>EXAMPLE</h2>
+
+Basic example using the default kernel size of 5:
+<div class="code"><pre>
+r.local.relief input=elevation output=lrm5
+</pre></div>
+Example with a custom kernel size of 11:
+<div class="code"><pre>
+r.local.relief input=elevation output_prefix=lrm11 neighborhood_size=11
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<em>
+<a href="r.shaded.relief.html">r.shaded.relief</a>,
+<a href="r.shaded.pca.html">r.shaded.pca</a> (in GRASS Addons)
+</em>
+
+
+<h2>REFERENCES</h2>
+<ul>
+ <li>Hesse, Ralf (2010). LiDAR-derived Local Relief Models - a new tool for archaeological prospection. <em>Archaeological Prospection</em> 17:67-72.</li>
+ <li>Bennett, Rebecca (2011). <em>Archaeological Remote Sensing: Visualization and Analysis of grass-dominated environments using laser scanning and digital spectra data.</em> Unpublished PhD Thesis. Electronic Document, <a href="http://eprints.bournemouth.ac.uk/20459/1/Bennett%2C_Rebecca_PhD_Thesis_2011.pdf">http://eprints.bournemouth.ac.uk/20459/1/Bennett%2C_Rebecca_PhD_Thesis_2011.pdf</a>, Accessed 25 February 2013. (provided bash script with <a href="v.surf.bspline.html">v.surf.bspline</a>-based implementation)</li>
+</ul>
+
+<h2>AUTHOR</h2>
+<p>
+Vaclav Petras, <a href="http://gis.ncsu.edu/osgeorel/">NCSU OSGeoREL</a>,<br/>
+Eric Goddard
+</p>
+
+<p><i>Last changed: $Date$</i>
Property changes on: grass-addons/grass7/raster/r.local.relief/r.local.relief.html
___________________________________________________________________
Added: svn:mime-type
+ text/html
Added: svn:keywords
+ Author Date Id
Added: svn:eol-style
+ native
Added: grass-addons/grass7/raster/r.local.relief/r.local.relief.py
===================================================================
--- grass-addons/grass7/raster/r.local.relief/r.local.relief.py (rev 0)
+++ grass-addons/grass7/raster/r.local.relief/r.local.relief.py 2013-09-27 00:18:12 UTC (rev 57849)
@@ -0,0 +1,208 @@
+#!/usr/bin/env python
+
+############################################################################
+#
+# MODULE: r.local.relief
+# AUTHOR(S): Vaclav Petras <wenzeslaus gmail.com>,
+# Eric Goddard <egoddard memphis.edu>
+# PURPOSE: Create a local relief models from elevation map
+# 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 a local relief models from elevation map.
+#% keywords: raster
+#% keywords: elevation
+#% keywords: terrain
+#% keywords: relief
+#%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 local relief map
+#% required: yes
+#%end
+#%option
+#% key: neighborhood_size
+#% type: integer
+#% description: Neighborhood size used when smoothing the elevation model
+#% options: 0-
+#% answer: 3
+#%end
+
+# TODO: The step 5 (contours to smooth elevation) can be replaced by using
+# vector tools (v.surf.*), this needs further testing.
+# Note, that quality of interpolation highly determines the result.
+
+# TODO: implement save intermediate results flag
+# (code is ready, but consider also outputting only some useful ones)
+
+import os
+import atexit
+
+import grass.script as gscript
+import grass.script.core as gcore
+
+
+RREMOVE = []
+VREMOVE = []
+
+
+def cleanup():
+ if RREMOVE or VREMOVE:
+ gcore.info(_("Cleaning temporary maps..."))
+ for rast in RREMOVE:
+ gscript.run_command('g.remove', rast=rast, quiet=True)
+ for vect in VREMOVE:
+ gscript.run_command('g.remove', vect=vect, quiet=True)
+
+
+def create_tmp_map_name(name):
+ return '{mod}_{pid}_{map_}_tmp'.format(mod='r_shaded_pca',
+ pid=os.getpid(),
+ map_=name)
+
+
+def create_persistent_map_name(basename, name):
+ return '{basename}_{name}'.format(basename=basename,
+ name=name)
+
+
+def check_map_name(name, mapset, element_type):
+ if gscript.find_file(name, element=element_type, mapset=mapset)['file']:
+ gscript.fatal(_("Raster map <%s> already exists. "
+ "Change the base name or allow overwrite.") % name)
+
+
+def main():
+ atexit.register(cleanup)
+ options, unused = gscript.parser()
+
+ elevation_input = options['input']
+ local_relief_output = options['output']
+ neighborhood_size = int(options['neighborhood_size'])
+ fill_method = 'cubic'
+
+ color_table = 'differences'
+ save_intermediates = True
+
+ if save_intermediates:
+ def local_create(name):
+ """create_persistent_map_name with hard-coded first argument"""
+ basename = local_relief_output.split('@')[0]
+ return create_persistent_map_name(basename=basename, name=name)
+ create_map_name = local_create
+ else:
+ create_map_name = create_tmp_map_name
+
+ smooth_elevation = create_map_name('smooth_elevation')
+ subtracted_smooth_elevation = create_map_name('subtracted_smooth_elevation')
+ vector_contours = create_map_name('vector_contours')
+ raster_contours = create_map_name('raster_contours')
+ raster_contours_with_values = create_map_name('raster_contours_with_values')
+ purged_elevation = create_map_name('purged_elevation')
+
+ if not save_intermediates:
+ RREMOVE.append(smooth_elevation)
+ RREMOVE.append(subtracted_smooth_elevation)
+ VREMOVE.append(vector_contours)
+ RREMOVE.append(raster_contours)
+ RREMOVE.append(raster_contours_with_values)
+ RREMOVE.append(purged_elevation)
+
+ # check even for the temporary maps
+ # (although, in ideal world, we should always fail if some of them exists)
+ if not gcore.overwrite():
+ check_map_name(smooth_elevation, gscript.gisenv()['MAPSET'], 'cell')
+ check_map_name(subtracted_smooth_elevation, gscript.gisenv()['MAPSET'], 'cell')
+ check_map_name(vector_contours, gscript.gisenv()['MAPSET'], 'vect')
+ check_map_name(raster_contours, gscript.gisenv()['MAPSET'], 'cell')
+ check_map_name(raster_contours_with_values, gscript.gisenv()['MAPSET'], 'cell')
+ check_map_name(purged_elevation, gscript.gisenv()['MAPSET'], 'cell')
+
+ # algorithm according to Hesse 2010 (LiDAR-derived Local Relief Models)
+ # step 1 (point cloud to digital elevation model) omitted
+
+ # step 2
+ gscript.info(_("Smoothing using r.neighbors..."))
+ gscript.run_command('r.neighbors', input=elevation_input,
+ output=smooth_elevation,
+ size=neighborhood_size,
+ overwrite=gcore.overwrite())
+
+ # step 3
+ gscript.info(_("Subtracting smoothed from original elevation..."))
+ gscript.mapcalc('{c} = {a} - {b}'.format(c=subtracted_smooth_elevation,
+ a=elevation_input,
+ b=smooth_elevation,
+ overwrite=gcore.overwrite()))
+
+ # step 4
+ gscript.info(_("Finding zero contours in elevation difference map..."))
+ gscript.run_command('r.contour', input=subtracted_smooth_elevation,
+ output=vector_contours,
+ levels=[0],
+ overwrite=gcore.overwrite())
+
+ # step 5
+ gscript.info(_("Extracting z value from the elevation"
+ " for difference zero contours..."))
+ gscript.run_command('v.to.rast', input=vector_contours,
+ output=raster_contours,
+ type='line', use='val', value=1,
+ overwrite=gcore.overwrite())
+ gscript.mapcalc('{c} = {a} * {b}'.format(c=raster_contours_with_values,
+ a=raster_contours,
+ b=elevation_input,
+ overwrite=gcore.overwrite()))
+
+ gscript.info(_("Interpolating elevation between"
+ " difference zero contours..."))
+ gscript.run_command('r.fillnulls', input=raster_contours_with_values,
+ output=purged_elevation,
+ method=fill_method,
+ overwrite=gcore.overwrite())
+
+ # step 6
+ gscript.info(_("Subtracting purged from original elevation..."))
+ gscript.mapcalc('{c} = {a} - {b}'.format(c=local_relief_output,
+ a=elevation_input,
+ b=purged_elevation,
+ overwrite=gcore.overwrite()))
+
+ if save_intermediates:
+ # same color table as input
+ gscript.run_command('r.colors', map=smooth_elevation,
+ raster=elevation_input, quiet=True)
+ gscript.run_command('r.colors', map=raster_contours_with_values,
+ raster=elevation_input, quiet=True)
+ gscript.run_command('r.colors', map=purged_elevation,
+ raster=elevation_input, quiet=True)
+
+ # has only one color
+ gscript.run_command('r.colors', map=raster_contours,
+ color='grey', quiet=True)
+
+ # same color table as output
+ gscript.run_command('r.colors', map=subtracted_smooth_elevation,
+ color=color_table, quiet=True)
+
+ gscript.run_command('r.colors', map=local_relief_output,
+ color=color_table, quiet=True)
+
+
+if __name__ == "__main__":
+ main()
Property changes on: grass-addons/grass7/raster/r.local.relief/r.local.relief.py
___________________________________________________________________
Added: svn:executable
+ *
Added: svn:mime-type
+ text/x-python
Added: svn:eol-style
+ native
More information about the grass-commit
mailing list