소스 검색

add support to v.surf.spline

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@50128 15284696-431f-4ddb-bdfa-cd5b030d7da7
Luca Delucchi 13 년 전
부모
커밋
33a6be8dc4
2개의 변경된 파일52개의 추가작업 그리고 22개의 파일을 삭제
  1. 13 7
      scripts/r.fillnulls/r.fillnulls.html
  2. 39 15
      scripts/r.fillnulls/r.fillnulls.py

+ 13 - 7
scripts/r.fillnulls/r.fillnulls.html

@@ -2,8 +2,8 @@
 
 
 <em>r.fillnulls</em> fills NULL pixels (no data areas) in input map and
 <em>r.fillnulls</em> fills NULL pixels (no data areas) in input map and
 stores filled map to a new output map. The fill data are interpolated 
 stores filled map to a new output map. The fill data are interpolated 
-from the no data area boundaries buffer using <em>v.surf.rst</em>
-spline interpolation.
+from the no data area boundaries buffer using <em>v.surf.rst</em>  or 
+<em>v.surf.bspline</em> spline interpolation.
 
 
 <h2>NOTES</h2>
 <h2>NOTES</h2>
 
 
@@ -11,7 +11,7 @@ Each area boundary buffer is set to three times the map resolution to get nomina
 three points around the edge. This way the algorithm interpolates into the hole with
 three points around the edge. This way the algorithm interpolates into the hole with
 a trained slope and curvature at the edges, in order to avoid that such a flat plane
 a trained slope and curvature at the edges, in order to avoid that such a flat plane
 is generated in a hole.
 is generated in a hole.
-<p>During the interpolation following warning may occur:<p>
+<p>During the interpolation following warning may occur when using the RST method:<p>
 <tt>
 <tt>
 Warning: strip exists with insufficient data<br>
 Warning: strip exists with insufficient data<br>
 Warning: taking too long to find points for interpolation--please change
 Warning: taking too long to find points for interpolation--please change
@@ -24,13 +24,16 @@ may pay attention to below notes.
 
 
 <h2>NOTES</h2>
 <h2>NOTES</h2>
 
 
-The algorithm is based on <em>v.surf.rst</em>
+When using the default RST method, the algorithm is based on <em>v.surf.rst</em>
 regularized splines with tension interpolation module which interpolates the
 regularized splines with tension interpolation module which interpolates the
 raster cell values for NULL data areas from the boundary values of the NULL
 raster cell values for NULL data areas from the boundary values of the NULL
 data area. An eventual raster MASK is respected during the NULL data area(s)
 data area. An eventual raster MASK is respected during the NULL data area(s)
 filling. The interpolated values are patched into the NULL data area(s) of
 filling. The interpolated values are patched into the NULL data area(s) of
 the input map and saved into a new raster map.
 the input map and saved into a new raster map.
 
 
+Otherwise, either the bilinear or bicubic method can be selected (based on
+<em>v.surf.bspline</em>). 
+
 <h2>WARNING</h2>
 <h2>WARNING</h2>
 
 
 Depending on the shape of the NULL data area(s) problems may occur due to an
 Depending on the shape of the NULL data area(s) problems may occur due to an
@@ -38,7 +41,8 @@ insufficient number of input cell values for the interpolation process. Most
 problems will occur if a NULL data area reaches a large amount of the map
 problems will occur if a NULL data area reaches a large amount of the map
 boundary. The user will have to carefully check the result using
 boundary. The user will have to carefully check the result using
 <em>r.mapcalc</em> (generating a difference map to the
 <em>r.mapcalc</em> (generating a difference map to the
-input map) and/or <em>d.what.rast</em> to query individual cell values.
+input map and applying the "differences" color table with <em>r.colors</em>)
+and/or <em>d.what.rast</em> to query individual cell values.
 
 
 <h2>EXAMPLE</h2>
 <h2>EXAMPLE</h2>
 
 
@@ -65,7 +69,8 @@ d.rast elev_srtm_30m_complete
 
 
 <em>
 <em>
 <a href="r.fill.dir.html">r.fill.dir</a>, 
 <a href="r.fill.dir.html">r.fill.dir</a>, 
-<a href="r.mapcalc.html">r.mapcalc</a>, 
+<a href="r.mapcalc.html">r.mapcalc</a>,
+<a href="v.surf.bspline.html">v.surf.bspline</a>, 
 <a href="v.surf.rst.html">v.surf.rst</a>
 <a href="v.surf.rst.html">v.surf.rst</a>
 </em>
 </em>
 
 
@@ -91,7 +96,8 @@ II. Application to Terrain Modeling and Surface Geometry Analysis,
 <i>Mathematical Geology</i> 25, 657-667.
 <i>Mathematical Geology</i> 25, 657-667.
 
 
 <h2>AUTHORS</h2>
 <h2>AUTHORS</h2>
-r.fillnulls: Markus Neteler, University of Hannover<p>and authors of v.surf.rst<br>
+r.fillnulls: Markus Neteler, University of Hannover  and Fondazione Edmund Mach
+<p>and authors of v.surf.rst<br>
 Improvement by Hamish Bowman, NZ
 Improvement by Hamish Bowman, NZ
 
 
 <p><i>Last changed: $Date$</i>
 <p><i>Last changed: $Date$</i>

+ 39 - 15
scripts/r.fillnulls/r.fillnulls.py

@@ -3,15 +3,16 @@
 ############################################################################
 ############################################################################
 #
 #
 # MODULE:	r.fillnulls
 # MODULE:	r.fillnulls
-# AUTHOR(S):	Markus Neteler <neteler itc it>
+# AUTHOR(S):	Markus Neteler
 #               Updated to GRASS 5.7 by Michael Barton
 #               Updated to GRASS 5.7 by Michael Barton
 #               Updated to GRASS 6.0 by Markus Neteler
 #               Updated to GRASS 6.0 by Markus Neteler
 #               Ring improvements by Hamish Bowman
 #               Ring improvements by Hamish Bowman
 #               Converted to Python by Glynn Clements
 #               Converted to Python by Glynn Clements
+#               Add support to v.surf.bspline by Luca Delucchi
 # PURPOSE:	fills NULL (no data areas) in raster maps
 # PURPOSE:	fills NULL (no data areas) in raster maps
 #               The script respects a user mask (MASK) if present.
 #               The script respects a user mask (MASK) if present.
 #
 #
-# COPYRIGHT:	(C) 2001,2004-2005,2008 by the GRASS Development Team
+# COPYRIGHT:	(C) 2001-2012 by the GRASS Development Team
 #
 #
 #		This program is free software under the GNU General Public
 #		This program is free software under the GNU General Public
 #		License (>=v2). Read the file COPYING that comes with GRASS
 #		License (>=v2). Read the file COPYING that comes with GRASS
@@ -21,7 +22,7 @@
 
 
 
 
 #%module
 #%module
-#% description: Fills no-data areas in raster maps using splines interpolation.
+#% description: Fills no-data areas in raster maps using spline interpolation.
 #% keywords: raster
 #% keywords: raster
 #% keywords: elevation
 #% keywords: elevation
 #% keywords: interpolation
 #% keywords: interpolation
@@ -44,6 +45,14 @@
 #% required : no
 #% required : no
 #% answer : 0.1
 #% answer : 0.1
 #%end
 #%end
+#%option
+#% key: method
+#% type: string
+#% description: Interpolation method
+#% required : yes
+#% options : bilinear,bicubic,rst
+#% answer : rst
+#%end
 
 
 import sys
 import sys
 import os
 import os
@@ -75,6 +84,7 @@ def main():
     output = options['output']
     output = options['output']
     tension = options['tension']
     tension = options['tension']
     smooth = options['smooth']
     smooth = options['smooth']
+    method = options['method']
 
 
     mapset = grass.gisenv()['MAPSET']
     mapset = grass.gisenv()['MAPSET']
     unique = str(os.getpid())
     unique = str(os.getpid())
@@ -140,30 +150,44 @@ def main():
     # remove internal MASK first -- WHY???? MN 10/2005
     # remove internal MASK first -- WHY???? MN 10/2005
     grass.run_command('g.remove', quiet = True, rast = 'MASK')
     grass.run_command('g.remove', quiet = True, rast = 'MASK')
 
 
+    # print message is a usermask it was present
     if grass.find_file(usermask, mapset = mapset)['file']:
     if grass.find_file(usermask, mapset = mapset)['file']:
 	grass.message(_("Using user mask while interpolating"))
 	grass.message(_("Using user mask while interpolating"))
 	maskmap = usermask
 	maskmap = usermask
     else:
     else:
 	maskmap = None
 	maskmap = None
 
 
-    segmax = 600
-    if pointsnumber > segmax:
-	grass.message(_("Using segmentation for interpolation..."))
-	segmax = None
-    else:
-	grass.message(_("Using no segmentation for interpolation as not needed..."))
-
-    grass.run_command('v.surf.rst', input = vecttmp, elev = tmp1 + '_filled',
-		      zcol = 'value', tension = tension, smooth = smooth,
-		      maskmap = maskmap, segmax = segmax)
-
-    grass.message(_("Note: Above warnings may be ignored."))
+    #check if method is rst to use v.surf.rst
+    if method == 'rst':
+        # set the max number before segmantation
+        segmax = 600
+        if pointsnumber > segmax:
+            grass.message(_("Using segmentation for interpolation..."))
+            segmax = None
+        else:
+            grass.message(_("Using no segmentation for interpolation as not needed..."))
+        # launch v.surf.rst    
+	grass.message(_("Using RST interpolation..."))
+	grass.run_command('v.surf.rst', input = vecttmp, elev = tmp1 + '_filled',
+			zcol = 'value', tension = tension, smooth = smooth,
+			maskmap = maskmap, segmax = segmax)
+
+	grass.message(_("Note: Above warnings may be ignored."))
 
 
     # restoring user's mask, if present:
     # restoring user's mask, if present:
     if grass.find_file(usermask, mapset = mapset)['file']:
     if grass.find_file(usermask, mapset = mapset)['file']:
 	grass.message(_("Restoring user mask (MASK)..."))
 	grass.message(_("Restoring user mask (MASK)..."))
 	grass.run_command('g.rename', quiet = True, rast = (usermask, 'MASK'))
 	grass.run_command('g.rename', quiet = True, rast = (usermask, 'MASK'))
 
 
+    #check if method is different from rst to use v.surf.bspline
+    if method != 'rst':
+	grass.message(_("Using %s (v.surf.bspline) interpolation") % method)
+	reg = grass.region()
+	# launch v.surf.bspline
+	grass.run_command('v.surf.bspline', input = vecttmp, layer = 1,
+			raster = tmp1 + '_filled', method = method, column = 'value',
+			sie = reg['ewres'], sin = reg['nsres'])
+
     # patch orig and fill map
     # patch orig and fill map
     grass.message(_("Patching fill data into NULL areas..."))
     grass.message(_("Patching fill data into NULL areas..."))
     # we can use --o here as g.parser already checks on startup
     # we can use --o here as g.parser already checks on startup