[GRASS-SVN] r57972 - grass-addons/grass7/raster/r.local.relief
svn_grass at osgeo.org
svn_grass at osgeo.org
Thu Oct 10 12:45:07 PDT 2013
Author: wenzeslaus
Date: 2013-10-10 12:45:07 -0700 (Thu, 10 Oct 2013)
New Revision: 57972
Modified:
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: adding bspline interpolation (author: Eric Goddard)
Modified: grass-addons/grass7/raster/r.local.relief/r.local.relief.html
===================================================================
--- grass-addons/grass7/raster/r.local.relief/r.local.relief.html 2013-10-10 11:56:20 UTC (rev 57971)
+++ grass-addons/grass7/raster/r.local.relief/r.local.relief.html 2013-10-10 19:45:07 UTC (rev 57972)
@@ -9,21 +9,36 @@
<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 input DEM's elevation values along the zero contour lines.</li>
- <li>Create a purged DEM by interpolating the null values between the rasterized contours generated in the previous step. This layer represents the large-scale landforms that will be removed to expose the local relief in the final step. This step is performed by the <a href="r.fillnulls.html">r.fillnulls</a> module using cubic spline interpolation.</li>
- <li>Subtract the purged DEM from the original DEM to get the local relief model.
+ <li>Create a purged DEM by interpolating the null values between the rasterized contours generated in the previous step. This layer represents the large-scale landforms that will be removed to expose the local relief in the final step.</li>
+ <li>Subtract the purged DEM from the original DEM to get the local relief model.</li>
</ol>
+<p>
+The interpolation step is performed by the <a href="r.fillnulls.html">r.fillnulls</a> module by default (using cubic interpolation).
+If this is not working on your data, you can use <em>-v</em> flag to use <a href="v.surf.bspline.html">v.surf.bspline</a> cubic interpolation instead (this might be slower on some types of data).
+</p>
<h2>OUTPUT</h2>
<p>
The final local relief model is named according to the <em>output</em> parameter.
-When the <em>-i</em> flag is specified, <em>r.local.relief</em> creates four additional output files representing the intermediate steps in the LRM generation process. The names of the intermediate fiels are composed of the user-specified <em>output</em> parameter and the following suffixes:</p>
+When the <em>-i</em> flag is specified, <em>r.local.relief</em> creates additional output files representing the intermediate steps in the LRM generation process. The names and number of the intermediate files vary depending on whether <a href="r.fillnulls.html">r.fillnulls</a> (default) or <a href="v.surf.bspline.html">v.surf.bspline</a> (specified by using the <em>-v</em> flag) is used for interpolation. The intermediate maps are composed of the user-specified <em>output</em> parameter and suffixes describing the intermediate map.</p>
+
+<p>Without using the <em>-v</em> flag (<a href="r.fillnulls.html">r.fillnulls</a> interpolation), intermediate maps have 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</tt>: The result of subtracting the low pass filter map from the DEM.</li>
- <li><tt>raster_contours_with_values</tt>: The rasterized zero contours with the values from elevation map.</li>
- <li><tt>_purged_elevation</tt>: The raster interpolated from the raster map based on zero contours map.</li>
+ <li><tt>_raster_contours_with_values</tt>: The rasterized zero contours with the values from elevation map.</li>
+ <li><tt>_purged_elevation</tt>: The raster interpolated from the _raster_contours_with_values map based that represents the large-scale landforms.</li>
</ul>
+<p>With using the <em>-v</em> flag (<a href="v.surf.bspline.html">v.surf.bspline</a> interpolation), intermediate maps have 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</tt>: The result of subtracting the low pass filter map from the DEM.</li>
+ <li><tt>_vector_contours</tt>: The zero contours extracted from the DEM.</li>
+ <li><tt>_contour_points</tt>: The points extractacted along the zero contour lines with the input DEM elevation values.</li>
+ <li><tt>_purged_elevation</tt>: The raster interpolated from the _contour_points map that represents the large-scale landforms.</li>
+</ul>
+
<h2>EXAMPLE</h2>
Basic example using the default neighborhood size of 11:
@@ -38,6 +53,10 @@
<div class="code"><pre>
r.local.relief -i input=elevation output=lrm11
</pre></div>
+Example using the default neighborhood size of 11 with bspline interpolation and saving the intermediate maps:
+<div class="code"><pre>
+r.local.relief -i -v input=elevation output=lrm11
+</pre></div>
<h2>SEE ALSO</h2>
Modified: grass-addons/grass7/raster/r.local.relief/r.local.relief.py
===================================================================
--- grass-addons/grass7/raster/r.local.relief/r.local.relief.py 2013-10-10 11:56:20 UTC (rev 57971)
+++ grass-addons/grass7/raster/r.local.relief/r.local.relief.py 2013-10-10 19:45:07 UTC (rev 57972)
@@ -40,6 +40,7 @@
#%option
#% key: neighborhood_size
#% type: integer
+#% label: Smoothing neighborhood size
#% description: Neighborhood size used when smoothing the elevation model
#% options: 0-
#% answer: 11
@@ -48,12 +49,12 @@
#% key: i
#% description: Save intermediate maps
#%end
+#%flag
+#% key: v
+#% label: Use bspline interpolation
+#% description: Uses v.surf.bspline cubic interpolation instead of r.fillnulls cubic interpolation.
+#%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.
-
-
import os
import atexit
@@ -99,9 +100,15 @@
local_relief_output = options['output']
neighborhood_size = int(options['neighborhood_size'])
save_intermediates = flags['i']
+ bspline = flags['v'] # when bspline == False, r.fillnulls is used
+ # constants
fill_method = 'cubic'
- color_table = 'differences'
+ # color table changed from difference to grey to histogram equalized-grey
+ # It does make more sense to use that since many archaeologists use the same
+ # color scheme for magnetometry and gpr data.
+ color_table = 'grey'
+ rcolors_flags = 'e'
if save_intermediates:
def local_create(name):
@@ -115,21 +122,29 @@
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 bspline:
+ contour_points = create_map_name('contour_points')
+ else:
+ raster_contours = create_map_name('raster_contours')
+
# if saving intermediates, keep only 1 contour layer
if save_intermediates:
- RREMOVE.append(raster_contours)
- VREMOVE.append(vector_contours)
+ if not bspline:
+ RREMOVE.append(raster_contours)
+ VREMOVE.append(vector_contours)
else:
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)
+ if bspline:
+ VREMOVE.append(contour_points)
+ else:
+ RREMOVE.append(raster_contours)
+ RREMOVE.append(raster_contours_with_values)
# check even for the temporary maps
# (although, in ideal world, we should always fail if some of them exists)
@@ -137,9 +152,12 @@
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')
+ if bspline:
+ check_map_name(contour_points, gscript.gisenv()['MAPSET'], 'vect')
+ else:
+ check_map_name(raster_contours, gscript.gisenv()['MAPSET'], 'cell')
# algorithm according to Hesse 2010 (LiDAR-derived Local Relief Models)
# step 1 (point cloud to digital elevation model) omitted
@@ -165,25 +183,51 @@
levels=[0],
overwrite=gcore.overwrite())
+ # Diverge here if using bspline interpolation
# 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()))
+ if bspline:
+ # Extract points from vector contours
+ gscript.run_command('v.to.points', _input=vector_contours, llayer='1',
+ _type='line', output=contour_points, dmax='10',
+ 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())
+ # Extract original dem elevations at point locations
+ gscript.run_command('v.what.rast', _map=contour_points,
+ raster=elevation_input, layer='2', column='along')
+ # Get mean distance between points to optimize spline interpolation
+ mean_dist = gscript.parse_command('v.surf.bspline', flags='e',
+ _input=contour_points,
+ raster_output=purged_elevation,
+ layer='2', column='along',
+ method='cubic')
+ spline_step = round(float(mean_dist.keys()[0].split(" ")[-1])) * 2
+
+ gscript.info(_("Interpolating purged surface using a spline step value"
+ " of {s}".format(s=spline_step)))
+ gscript.run_command('v.surf.bspline', _input=contour_points,
+ raster_output=purged_elevation, layer='2',
+ column='along', method='cubic',
+ overwrite=gcore.overwrite())
+ else:
+ 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,
@@ -191,24 +235,24 @@
b=purged_elevation,
overwrite=gcore.overwrite()))
+ # set color tables
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)
+ if not bspline:
+ gscript.run_command('r.colors', map=raster_contours_with_values,
+ raster=elevation_input, 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,
+ gscript.run_command('r.colors', flags=rcolors_flags, map=local_relief_output,
color=color_table, quiet=True)
More information about the grass-commit
mailing list