[GRASS-SVN] r71185 - grass-addons/grass7/raster/r.mblend

[svn_grass at osgeo.org]

Author: Lads
Date: 2017-06-14 00:29:46 -0700 (Wed, 14 Jun 2017)
New Revision: 71185

Modified:
   grass-addons/grass7/raster/r.mblend/r.mblend.py
Log:


Modified: grass-addons/grass7/raster/r.mblend/r.mblend.py
===================================================================
--- grass-addons/grass7/raster/r.mblend/r.mblend.py	2017-06-13 21:49:26 UTC (rev 71184)
+++ grass-addons/grass7/raster/r.mblend/r.mblend.py	2017-06-14 07:29:46 UTC (rev 71185)
@@ -48,6 +48,10 @@
 #% multiple: no
 #% required: no
 #%end
+#%flag
+#% key: a
+#% description: Assign the average difference between the two rasters to the far edge (instead of zero). 
+#%end
 
 import os
 import atexit
@@ -55,6 +59,7 @@
 import grass.script as gscript
 
 index = 0
+far_edge_value = '0'
 TMP_MAPS = []
 WEIGHT_MAX = 10000
 COL_VALUE = 'value'
@@ -77,6 +82,7 @@
 
 
 def main():
+    global far_edge_value
 
     options, flags = gscript.parser()
     high = options['high']
@@ -84,6 +90,7 @@
     output = options['output']
     far_edge = float(options['far_edge'])
     inter_points = int(options['inter_points'])
+    use_average_differences = flags['a']
 
     if(high is None or high == ""):
         gscript.error(_("[r.mblend] ERROR: high is a mandatory parameter."))
@@ -151,6 +158,17 @@
                       " rasters"))
     gscript.run_command('r.to.vect', input=diff, output=diff_points,
                         type='point')
+    
+    # Compute average of the differences if flag -a was passed
+    if use_average_differences:
+        p = gscript.pipe_command('r.univar', map=diff)
+        result = {}
+        for line in p.stdout:
+            vector = line.split(": ")
+            if vector[0] == "mean":
+                print("Found it: " + vector[1])
+                far_edge_value = vector[1]
+        p.wait()
 
     # Obtain edge points of the high resolution raster
     interpol_area_buff = getTemporaryIdentifier()
@@ -208,34 +226,28 @@
     points_edges = getTemporaryIdentifier()
     gscript.message(_("[r.mblend] Setting far edge weights to zero"))
     gscript.run_command('v.db.update', map=weight_points_edge,
-                        column=COL_VALUE, value='0')
+                        column=COL_VALUE, value=far_edge_value)
     gscript.message(_("[r.mblend] Patching the two edges"))
     gscript.run_command('v.patch',
                         input=weight_points_edge + ',' + diff_points_edge,
                         output=points_edges, flags='e')
 
     # Interpolate stitching raster
-    stitching_full = getTemporaryIdentifier()
-    interpol_area_mask = getTemporaryIdentifier()
-    stitching = getTemporaryIdentifier()
+    smoothing = getTemporaryIdentifier()
+    # Consign region to interpolation area
+    gscript.run_command('g.region', vector=interpol_area_buff)
     gscript.message(_("[r.mblend] Interpolating smoothing surface. This" +
                       " might take a while..."))
     gscript.run_command('v.surf.idw', input=points_edges, column=COL_VALUE,
-                        output=stitching_full, power=2, npoints=inter_points)
-    # Create mask
-    gscript.message(_("[r.mblend] Creating mask for the interpolation area"))
-    gscript.run_command('v.to.rast', input=interpol_area,
-                        output=interpol_area_mask, use='val', value=1)
-    # Crop to area of interest
-    gscript.message(_("[r.mblend] Cropping the mask"))
-    gscript.mapcalc(stitching + ' = if(' + interpol_area_mask + ',' +
-                    stitching_full + ')')
+                        output=smoothing, power=2, npoints=inter_points)
+    # Reset region to full extent
+    gscript.run_command('g.region', raster=high + "," + low)
 
     # Apply stitching
     smooth_low_res = getTemporaryIdentifier()
     # Sum to low res
     gscript.message(_("[r.mblend] Applying smoothing surface"))
-    gscript.mapcalc(smooth_low_res + ' = ' + low_res_inter + ' + ' + stitching)
+    gscript.mapcalc(smooth_low_res + ' = ' + low_res_inter + ' + ' + smoothing)
     # Add both rasters
     try:
         gscript.message(_("[r.mblend] Joining result into a single raster"))