[GRASS-SVN] r72613 - grass-addons/grass7/raster/r.subdayprecip.design

[svn_grass at osgeo.org]

Author: martinl
Date: 2018-04-10 02:48:14 -0700 (Tue, 10 Apr 2018)
New Revision: 72613

Modified:
   grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.html
   grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.py
Log:
r.subdayprecip.design: implement reduction for large areas

Modified: grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.html
===================================================================
--- grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.html	2018-04-09 17:24:56 UTC (rev 72612)
+++ grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.html	2018-04-10 09:48:14 UTC (rev 72613)
@@ -14,6 +14,18 @@
 
 <h2>NOTES</h2>
 
+By default areas above <b>area_size</b> limit are not processed
+(calculated value is set to -1 in such case). Processing of large
+areas can be enabled by specifying <b>-k</b> flag. In this case calculated
+avarage value for large areas is reduced by coefficient below.
+
+<pre>
+  k = exp(-0.08515989 * (x^2) - 0.001344925 * (x^4))
+
+  where x = log10(area_size_km2) - 0.9
+</pre>
+
+<p>
 Subday design precipitation series are important for hydrological
 modelling and soil erosion problems in a small catchment scale when
 designing common measures for promoting water retention, landscape

Modified: grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.py
===================================================================
--- grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.py	2018-04-09 17:24:56 UTC (rev 72612)
+++ grass-addons/grass7/raster/r.subdayprecip.design/r.subdayprecip.design.py	2018-04-10 09:48:14 UTC (rev 72613)
@@ -8,7 +8,7 @@
 #
 # PURPOSE:      Computes subday design precipitation totals.
 #
-# COPYRIGHT:    (C) 2015 Martin Landa and GRASS development team
+# COPYRIGHT:    (C) 2015-2018 Martin Landa and GRASS development team
 #
 #               This program is free software under the GNU General
 #               Public License (>=v2). Read the file COPYING that
@@ -24,12 +24,12 @@
 #%end
 
 #%option G_OPT_V_MAP
-#% label: Vector map of location under analysis
+#% label: Name of input vector map of location under analysis
 #%end
 
 #%option G_OPT_R_INPUTS
 #% key: return_period
-#% description: Rainfall raster maps of required return period
+#% description: Name of input rainfall raster maps of required return period
 #% options: N2,N5,N10,N20,N50,N100
 #%end
 
@@ -48,11 +48,18 @@
 #% answer: 20
 #%end
 
+#%flag
+#% key: k
+#% description: Process also large areas above area_size limit by applying reduction
+#%end
+
 import os
 import sys
+import math
 
 import grass.script as grass 
 from grass.pygrass.modules import Module
+from grass.pygrass.vector import VectorTopo
 from grass.exceptions import CalledModuleError
 
 def coeff(name, rl):
@@ -120,6 +127,22 @@
 
     return a, c
 
+
+def area_size_reduction(map_name, field_name, area_col_name):
+    vmap = VectorTopo(map_name)
+    vmap.open('rw')
+
+    cats = [] # TODO: do it better
+    for feat in vmap.viter('areas'):
+        if feat.attrs['cat'] not in cats:
+            x = math.log10(float(feat.attrs[area_col_name]) )- 0.9
+            k = math.exp(-0.08515989 * pow(x, 2) - 0.001344925 * pow(x, 4))
+            feat.attrs[field_name] *= k
+            cats.append(feat.attrs['cat'])
+
+    vmap.table.conn.commit()  
+    vmap.close()
+    
 def main():
     # check if the map is in the current mapset
     mapset = grass.find_file(opt['map'], element='vector')['mapset']
@@ -151,7 +174,7 @@
                        where='{} > {}'.format(area_col_name, opt['area_size']),
                        stdout_=grass.PIPE)
         large_areas = len(areas.outputs.stdout.splitlines())
-        if large_areas > 0:
+        if large_areas > 0 and not flg['k']:
             grass.warning('{} areas larger than size limit will be skipped from computation'.format(large_areas))
 
     # extract multi values to points
@@ -208,9 +231,15 @@
                column=field_name, query_column=expression)
 
         if check_area_size:
-            Module('v.db.update', map=opt['map'],
-                   column=field_name, value='-1',
-                   where='{} > {}'.format(area_col_name, opt['area_size']))
+            args = {}
+            if flg['k']:
+                area_size_reduction(opt['map'], field_name, area_col_name)
+            else:
+                Module('v.db.update', map=opt['map'],
+                       column=field_name,
+                       value='-1',
+                       where='{} > {}'.format(area_col_name, opt['area_size']),
+                       **args)
 
         # remove unused column
         Module('v.db.dropcolumn', map=opt['map'],