[GRASS-SVN] r50984 -
grass/branches/releasebranch_6_4/scripts/r.fillnulls
svn_grass at osgeo.org
svn_grass at osgeo.org
Sun Mar 4 06:55:32 EST 2012
Author: martinl
Date: 2012-03-04 03:55:32 -0800 (Sun, 04 Mar 2012)
New Revision: 50984
Modified:
grass/branches/releasebranch_6_4/scripts/r.fillnulls/description.html
grass/branches/releasebranch_6_4/scripts/r.fillnulls/r.fillnulls
Log:
r.fillnulls: major backport from devbr6 (see http://trac.osgeo.org/grass/wiki/Grass6Planning#GRASS6.4.3)
Modified: grass/branches/releasebranch_6_4/scripts/r.fillnulls/description.html
===================================================================
--- grass/branches/releasebranch_6_4/scripts/r.fillnulls/description.html 2012-03-04 11:52:56 UTC (rev 50983)
+++ grass/branches/releasebranch_6_4/scripts/r.fillnulls/description.html 2012-03-04 11:55:32 UTC (rev 50984)
@@ -2,8 +2,8 @@
<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
-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> interpolation.
<h2>NOTES</h2>
@@ -12,7 +12,7 @@
a trained slope and curvature at the edges, in order to avoid that such a flat plane
is generated in a hole.
<p>
-During the interpolation following warning may occur:<p>
+During the interpolation following warning may occur when using the RST method:<p>
<tt>
Warning: strip exists with insufficient data<br>
@@ -27,13 +27,16 @@
<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
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)
filling. The interpolated values are patched into the NULL data area(s) of
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>
Depending on the shape of the NULL data area(s) problems may occur due to an
@@ -41,7 +44,8 @@
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
<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>
@@ -69,6 +73,7 @@
<em>
<a href="r.fill.dir.html">r.fill.dir</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>
</em>
@@ -94,7 +99,7 @@
<i>Mathematical Geology</i> 25, 657-667.
<h2>AUTHORS</h2>
-r.fillnulls: Markus Neteler, University of Hannover<p>
+r.fillnulls: Markus Neteler, University of Hannover and Fondazione Edmund Mach<p>
and authors of v.surf.rst<br>
Improvement by Hamish Bowman, NZ
Modified: grass/branches/releasebranch_6_4/scripts/r.fillnulls/r.fillnulls
===================================================================
--- grass/branches/releasebranch_6_4/scripts/r.fillnulls/r.fillnulls 2012-03-04 11:52:56 UTC (rev 50983)
+++ grass/branches/releasebranch_6_4/scripts/r.fillnulls/r.fillnulls 2012-03-04 11:55:32 UTC (rev 50984)
@@ -3,14 +3,15 @@
############################################################################
#
# 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 6.0 by Markus Neteler
# Ring improvements by Hamish Bowman
+# v.surf.bspline support by Markus Neteler
# PURPOSE: fills NULL (no data areas) in raster maps
# The script respects a user mask (MASK) if present.
#
-# COPYRIGHT: (C) 2001,2004-2005 by the GRASS Development Team
+# COPYRIGHT: (C) 2001-2012 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
@@ -20,21 +21,21 @@
#%Module
-#% description: Fills no-data areas in raster maps using v.surf.rst splines interpolation
+#% description: Fills no-data areas in raster maps using spline interpolation.
#% keywords: raster, elevation, interpolation
#%End
#%option
#% key: input
#% gisprompt: old,cell,raster
#% type: string
-#% description: Raster map in which to fill nulls
+#% description: Name of input raster map in which to fill nulls
#% required : yes
#%end
#%option
#% key: output
#% gisprompt: new,cell,raster
#% type: string
-#% description: Output raster map with nulls filled by interpolation from surrounding values
+#% description: Name for output raster map with nulls filled by interpolation
#% required : yes
#%end
#%option
@@ -51,6 +52,14 @@
#% required : no
#% answer : 0.1
#%end
+#%option
+#% key: method
+#% type: string
+#% description: Interpolation method
+#% options:bilinear,bicubic,rst
+#% answer:rst
+#% required:yes
+#%end
if test "$GISBASE" = ""; then
echo "You must be in GRASS GIS to run this program." >&2
@@ -94,6 +103,11 @@
g.remove rast=MASK --quiet > /dev/null
g.rename "$USERMASK",MASK --quiet > /dev/null
fi
+
+ # remove temporary region
+ unset WIND_OVERRIDE
+ g.remove region="tmp_rfillnulls.$$" --quiet
+
exit 1
}
# shell check for user break (signal list: trap -l)
@@ -115,14 +129,18 @@
g.message -e "<$GIS_OPT_INPUT> does not exist! Aborting."
exit 1
fi
+case "$GIS_OPT_METHOD" in
+ rst | bilinear | bicubic) ;;
+ *) GIS_OPT_METHOD=rst ;;
+esac
#check if a MASK is already present:
MASKTMP=mask.$$
USERMASK="usermask_$MASKTMP"
eval `g.findfile element=cell file=MASK`
if [ "$file" ] ; then
- g.message "A user raster mask (MASK) is present. Saving it..."
- g.rename MASK,"$USERMASK" --quiet > /dev/null
+ g.message "A user raster mask (MASK) is present. Saving it..."
+ g.rename MASK,"$USERMASK" --quiet > /dev/null
fi
#make a mask of NULL cells
@@ -134,7 +152,7 @@
g.message "Locating and isolating NULL areas..."
#creating 0/1 map:
-r.mapcalc "$TMP1 = if(isnull($GIS_OPT_INPUT),1,null())"
+r.mapcalc "$TMP1 = if(isnull(\"$GIS_OPT_INPUT\"),1,null())"
#generate a ring:
# the buffer is set to three times the map resolution so you get nominally
@@ -147,6 +165,7 @@
# check-me: quoting ok?
RES="`echo $nsres $ewres | awk '{printf "%f", ($1+$2) * 3. / 2.}'`" # avg*3
+# maybe use r.grow here..
r.buffer input="$TMP1" dist="$RES" out="$TMP1.buf"
if [ $? -ne 0 ] ; then
g.message -e "$0 abandoned. Removing temporary map, restoring user mask if needed:"
@@ -162,7 +181,11 @@
# Use this outline (raster border) for interpolating the fill data:
VECTTMP=vecttmp_fillnulls_$$
g.message "Creating interpolation points..."
-r.to.vect input="$GIS_OPT_INPUT" output="$VECTTMP" feature=point
+
+## use the -b flag to avoid topology building on big jobs?
+## no, can't, 'g.region vect=' currently wants to see level 2
+r.to.vect -z input="$GIS_OPT_INPUT" output="$VECTTMP" feature=point
+
if [ $? -ne 0 ] ; then
g.message -e "$0 abandoned. Removing temporary maps, restoring user mask if needed:"
g.remove rast=MASK,"$TMP1","$TMP1.buf","${TMP1}_filled" --quiet > /dev/null
@@ -170,6 +193,8 @@
g.rename "$USERMASK",MASK --quiet > /dev/null
exit 1
fi
+# this mask not needed any more
+g.remove rast=MASK --quiet
#count number of points to control segmax parameter for interpolation:
# check-me: quoting?
@@ -185,8 +210,8 @@
#restoring user's mask, if present:
eval `g.findfile element=cell file=$USERMASK`
if [ "$file" ] ; then
- g.message "Restoring user mask (MASK)..."
- g.rename "$USERMASK",MASK --quiet > /dev/null
+ g.message "Restoring user mask (MASK)..."
+ g.rename "$USERMASK",MASK --quiet > /dev/null
fi
#cleanup
@@ -195,38 +220,71 @@
exit 1
fi
-g.message "Note: The following warnings may be ignored."
-#remove internal MASK first -- WHY???? MN 10/2005
-g.remove MASK > /dev/null
+if [ "$GIS_OPT_METHOD" = "rst" ] ; then
-eval `g.findfile element=cell file=$USERMASK`
-if [ "$file" ] ; then
- g.message "Using user mask while interpolating"
- RST_CMD="v.surf.rst zcol=value tension=$GIS_OPT_TENSION smooth=$GIS_OPT_SMOOTH maskmap=$USERMASK"
-else
- RST_CMD="v.surf.rst zcol=value tension=$GIS_OPT_TENSION smooth=$GIS_OPT_SMOOTH"
-fi
+ g.message "Using RST interpolation"
-SEGMAX=600
-if [ "$POINTSNUMBER" -ge "$SEGMAX" ] ; then
- g.message "Using segmentation for interpolation..."
- $RST_CMD input="$VECTTMP" elev="${TMP1}_filled" #2> /dev/null
+ ### zooming made it take slightly longer for bspline, so only use for rst
+ # setup internal region
+ g.region save="tmp_rfillnulls.$$"
+ WIND_OVERRIDE="tmp_rfillnulls.$$"
+ export WIND_OVERRIDE
+ g.region vect="$VECTTMP" align="$GIS_OPT_INPUT"
+
+ eval `g.findfile element=cell file=$USERMASK`
+ if [ "$file" ] ; then
+ g.message "Using user mask while interpolating"
+ RST_CMD="v.surf.rst -z layer=0 tension=$GIS_OPT_TENSION smooth=$GIS_OPT_SMOOTH maskmap=$USERMASK"
+ else
+ RST_CMD="v.surf.rst -z layer=0 tension=$GIS_OPT_TENSION smooth=$GIS_OPT_SMOOTH"
+ fi
+
+ g.message "Note: The following warnings may be ignored."
+
+ SEGMAX=600
+ if [ "$POINTSNUMBER" -ge "$SEGMAX" ] ; then
+ g.message "Using segmentation for interpolation..."
+ $RST_CMD input="$VECTTMP" elev="${TMP1}_filled" #2> /dev/null
+ else
+ #if less than $SEGMAX points, use no segmentation for speedup:
+ g.message "Not using segmentation for interpolation as it's not needed..."
+ $RST_CMD input="$VECTTMP" elev="${TMP1}_filled" segmax="$SEGMAX" #2> /dev/null
+ fi
+ g.message "Note: Above warnings may be ignored."
+
+ #restoring user's mask, if present:
+ eval `g.findfile element=cell file=$USERMASK`
+ if [ "$file" ] ; then
+ g.message "Restoring user mask (MASK)..."
+ g.rename "$USERMASK",MASK --quiet > /dev/null
+ fi
+
+ # restore real region
+ unset WIND_OVERRIDE
+ g.remove region="tmp_rfillnulls.$$" --quiet
+
else
- #if less than $SEGMAX points, use no segmentation for speedup:
- g.message "Using no segmentation for interpolation as not needed..."
- $RST_CMD input="$VECTTMP" elev="${TMP1}_filled" segmax="$SEGMAX" #2> /dev/null
-fi
+ g.message "Using $GIS_OPT_METHOD (v.surf.bspline) interpolation"
-g.message "Note: Above warnings may be ignored."
+ #restoring user's mask, if present
+ # pointless for now, since v.surf.bspline ignores MASK, but someday...
+ eval `g.findfile element=cell file=$USERMASK`
+ if [ "$file" ] ; then
+ g.message "Restoring user mask (MASK)..."
+ g.rename "$USERMASK",MASK --quiet
+ fi
-#restoring user's mask, if present:
-eval `g.findfile element=cell file=$USERMASK`
-if [ "$file" ] ; then
- g.message "Restoring user mask (MASK)..."
- g.rename "$USERMASK",MASK --quiet > /dev/null
+ eval `g.region -g`
+ SIE=`echo "$ewres" | awk '{printf "%f", $1 * 3.}'`
+ SIN=`echo "$nsres" | awk '{printf "%f", $1 * 3.}'`
+
+ v.surf.bspline input="$VECTTMP" raster="${TMP1}_filled" layer=0 \
+ method="$GIS_OPT_METHOD" sie="$SIE" sin="$SIN" lambda_i=0.005
+
fi
+
#patch orig and fill map
g.message "Patching fill data into NULL areas..."
# we can use --o here as g.parser already checks on startup
More information about the grass-commit
mailing list