[GRASS-SVN] r36834 - grass-addons/raster/r.surf.volcano

[svn_grass at osgeo.org]

Author: hamish
Date: 2009-04-21 05:40:47 -0400 (Tue, 21 Apr 2009)
New Revision: 36834

Modified:
   grass-addons/raster/r.surf.volcano/r.surf.volcano
Log:
add a flag to use a gaussian surface instead of 1/d^n

Modified: grass-addons/raster/r.surf.volcano/r.surf.volcano
===================================================================
--- grass-addons/raster/r.surf.volcano/r.surf.volcano	2009-04-21 08:27:43 UTC (rev 36833)
+++ grass-addons/raster/r.surf.volcano/r.surf.volcano	2009-04-21 09:40:47 UTC (rev 36834)
@@ -50,20 +50,48 @@
 #% type: integer
 #% required: no
 #% options: 1-25
-#% description: Friction of distance, (the 'n' in 1/d^n)
+#% label: Friction of distance, (the 'n' in 1/d^n)
+#% description: Higher values generate steeper slopes. (not used in Gaussian curve mode)
 #% answer: 6
 #%End
+
+#%Flag
+#% key: g
+#% description: Use a Gaussian curve instead of 1/(d^n) for radial basis function
+#%End
+# normalized Gaussian curve:  f(x) = a * e^( (x-b)^2 / 2*(c^2) )
+#  parameters: a = 1/(sigma*sqrt(2pi)), b = mu, c = sigma
+#  mu is mean value and sigma^2 is variance.
+#  so we only need b,c. and b can be locked at 0 here. so user only needs
+#  to give sigma (width)
+#  thus r.mapcalc expr could look like
+#   SIGMA=1.0   # beware variable namespace collision with roughness factor
+#   r.mapcalc doesn't have a pi(), so define it:  PI=3.14159265358979323846
+#   f(distance) = ( 1 / ($SIGMA*sqrt(2*$PI)) ) * exp( -1* $DIST^2 / (2 * $SIGMA^2) )
+#
+## FIXME: ok, it isn't really kurtosis but it's similar and I couldn't
+##        think of a better name.
 #%Option
-#% key: sigma
+#% key: kurtosis
 #% type: double
 #% required: no
-#% description: Surface roughness factor
+#% label: Peak steepness (used with Gaussian curve mode)
+#% description: Nearer to 0 is flatter, higher values generate steeper slopes.
 #% answer: 1.0
 #%End
+
 #%Flag
 #% key: r
 #% description: Roughen the surface
 #%End
+#%Option
+#% key: sigma
+#% type: double
+#% required: no
+#% label: Surface roughness factor
+#% description: Nearer to 0 is smoother, higher values make a rougher surface.
+#% answer: 1.0
+#%End
 
 ## TODO
 ##%Flag
@@ -71,20 +99,6 @@
 ##% description: Use (d^n)*log(d) instead of 1/(d^n) for radial basis function
 ##%End
 
-## TODO
-##%Flag
-##% key: g
-##% description: Use a Gaussian curve instead of 1/(d^n) for radial basis function
-##%End
-# normalized Gaussian curve:  f(x) = a * e^( (x-b)^2 / 2*(c^2) )
-#  parameters: a = 1/(sigma*sqrt(2pi)), b = mu, c = sigma
-#  mu is mean value and sigma^2 is variance.
-#  so we only need b,c. and b can be locked at 0 here. so user only needs
-#  to give sigma (width)
-#  thus r.mapcalc expr could look like
-#   SIGMA=1.0   # beware variable namespace collision with roughness factor
-#   r.mapcalc doesn't have a pi(), so define it:  PI=3.14159265358979323846
-#   f(distance) = ( 1 / ($SIGMA*sqrt(2*$PI)) ) * exp( $DIST^2 / (2 * $SIGMA^2) )
 
 
 if [ -z "$GISBASE" ] ; then
@@ -100,6 +114,7 @@
 PEAK="$GIS_OPT_PEAK"
 CRATER="$GIS_OPT_CRATER"
 SIGMA="$GIS_OPT_SIGMA"
+KURTOSIS="$GIS_OPT_KURTOSIS"
 
 
 # test for overwrite as r.mapcalc doesn't
@@ -119,16 +134,40 @@
    sqrt( (x() - $center_easting)^2 + (y() - $center_northing)^2 )"
 
 
-g.message -v "Normalizing cost map with 1 in the center and 0 at outer edge ..."
+g.message -v "Normalizing cost map ..."
 eval `r.info -r "volc.dist_units.$$"`
-r.mapcalc "volc.dist_norm.$$ = ($max - volc.dist_units.$$) / $max"
+if [ "$GIS_FLAG_G" -eq 0 ] ; then
+   # Normalize with 1 in the center and 0 at outer edge
+   r.mapcalc "volc.dist_norm.$$ = ($max - volc.dist_units.$$) / $max"
 
-
-g.message -v "Creating IDW surface ..."
-r.mapcalc "volc.peak.$$ = ($PEAK + abs($CRATER) ) \
+   g.message -v "Creating IDW surface ..."
+   r.mapcalc "volc.peak.$$ = ($PEAK + abs($CRATER) ) \
 	* pow( volc.dist_norm.$$, $FRICTION )"
+else
+   # Normalize with 0 in the center and 1 at outer edge
+   r.mapcalc "volc.dist_norm.$$ = volc.dist_units.$$ / $max"
 
+   g.message -v "Creating Gaussian surface ..."
+   SIGMA_C=1.0   # beware variable namespace collision with roughness factor
+   # r.mapcalc doesn't have a pi(), so define it:
+   PI=3.14159265358979323846
 
+   # f(distance) = ( 1 / ($SIGMA*sqrt(2*$PI)) ) * exp( -1* $DIST^2 / (2 * $SIGMA^2) )
+
+   ## FIXME: the 10*kurtosis stuff is a completely bogus hack!
+   r.mapcalc "volc.gauss.$$ = \
+     ( 1 / ( $SIGMA_C * sqrt(2 * $PI) ) ) \
+      * exp( -1* (10 * $KURTOSIS * volc.dist_norm.$$)^2 / (2 * $SIGMA_C^2) )"
+
+   eval `r.info -r "volc.gauss.$$"`
+   g.message -v "Normalizing Gaussian surface ..."
+   r.mapcalc "volc.peak.$$ = \
+     ( ($PEAK + abs($CRATER) ) / $max ) * volc.gauss.$$"
+
+   g.remove rast="volc.gauss.$$" --quiet
+fi
+
+
 if [ "$GIS_FLAG_R" -eq 0 ] ; then
    g.rename rast="volc.peak.$$,volc.surf.$$" --quiet
 else
@@ -170,14 +209,21 @@
 r.support map="$GIS_OPT_OUTPUT" \
     description="generated by r.surf.volcano" \
     source1="Peak height = $PEAK" \
-    source2="Friction of distance = $FRICTION" \
     title="Artificial surface resembling a seamount or cone volcano"
 
 if [ "$GIS_FLAG_R" -eq 1 ] ; then
    r.support map="$GIS_OPT_OUTPUT" \
-     history="Surface roughness used a Gaussian deviate with sigma of $SIGMA."
+      history="Surface roughness used a Gaussian deviate with sigma of $SIGMA."
 fi
 
+if [ "$GIS_FLAG_G" -eq 0 ] ; then
+   r.support map="$GIS_OPT_OUTPUT" \
+      source2="IDW surface with friction of distance = $FRICTION"
+else
+   r.support map="$GIS_OPT_OUTPUT" \
+      source2="Gaussian surface with pseudo-kurtosis factor = $KURTOSIS"
+fi
+
 g.message -v "Done."
 
 exit