fix region alignment, table key column, -nan, enable translations, and speed up the module

git-svn-id: https://svn.osgeo.org/grass/grass/trunk@44444 15284696-431f-4ddb-bdfa-cd5b030d7da7

scripts/v.rast.stats/v.rast.stats.py

@@ -3,7 +3,9 @@
 ############################################################################
 #
 # MODULE:	v.rast.stats
-# AUTHOR(S):	Markus Neteler, converted to Python by Glynn Clements
+# AUTHOR(S):	Markus Neteler
+#		converted to Python by Glynn Clements
+#		speed up by Markus Metz
 # PURPOSE:	Calculates univariate statistics from a GRASS raster map
 #		only for areas covered by vector objects on a per-category base
 # COPYRIGHT:	(C) 2005-2010 by the GRASS Development Team
@@ -78,7 +80,7 @@ def has_column(vector, col):
     return 
 
 def cleanup():
-    grass.run_command('g.remove', rast = '%s_%s' % (vector, tmpname), quiet = True)
+    grass.run_command('g.remove', rast = rastertmp, quiet = True)
     grass.run_command('g.remove', rast = 'MASK', quiet = True, stderr = nuldev)
     if mask_found:
 	grass.message(_("Restoring previous MASK..."))
@@ -87,7 +89,7 @@ def cleanup():
 #	grass.try_remove(f)
 
 def main():
-    global tmp, sqltmp, tmpname, nuldev, vector, mask_found
+    global tmp, sqltmp, tmpname, nuldev, vector, mask_found, rastertmp
     mask_found = False
     #### setup temporary files
     tmp = grass.tempfile()
@@ -119,47 +121,42 @@ def main():
 
     vector = vs[0]
 
-    #check the input raster map
+    rastertmp = "%s_%s" % (vector, tmpname)
+
+    # check the input raster map
     if not grass.find_file(raster, 'cell')['file']:
 	grass.fatal(_("Raster map <%s> not found") % raster)
 
-    #check presence of raster MASK, put it aside
+    # check presence of raster MASK, put it aside
     mask_found = bool(grass.find_file('MASK', 'cell')['file'])
     if mask_found:
 	grass.message(_("Raster MASK found, temporarily disabled"))
 	grass.run_command('g.rename', rast = ('MASK', tmpname + "_origmask"), quiet = True)
 
-    #get RASTER resolution of map which we want to query:
-    #fetch separated to permit for non-square cells (latlong etc)
-    kv = grass.raster_info(map = raster)
-    nsres = kv['nsres']
-    ewres = kv['ewres']
-
-    #save current settings:
+    # save current settings:
     grass.use_temp_region()
 
-    #Temporarily setting raster resolution to $RASTER resolution
-    #keep boundary settings
-    grass.run_command('g.region', flags = 'a', nsres = nsres, ewres = ewres)
+    # Temporarily aligning region resolution to $RASTER resolution
+    # keep boundary settings
+    grass.run_command('g.region', align = raster)
 
-    #prepare raster MASK
-    if grass.run_command('v.to.rast', input = vector, output = "%s_%s" % (vector, tmpname),
+    # prepare raster MASK
+    if grass.run_command('v.to.rast', input = vector, output = rastertmp,
 			 use = 'cat', quiet = True) != 0:
 	grass.fatal(_("An error occurred while converting vector to raster"))
 
-    #dump cats to file to avoid "too many argument" problem:
-    p = grass.pipe_command('r.category', map = '%s_%s' % (vector, tmpname), fs = ';', quiet = True)
+    # dump cats to file to avoid "too many argument" problem:
+    p = grass.pipe_command('r.category', map = rastertmp, fs = ';', quiet = True)
     cats = []
     for line in p.stdout:
 	cats.append(line.rstrip('\r\n').split(';')[0])
     p.wait()
 
-    #echo "List of categories found: $CATSLIST"
     number = len(cats)
     if number < 1:
 	grass.fatal(_("No categories found in raster map"))
 
-    #check if DBF driver used, in this case cut to 10 chars col names:
+    # check if DBF driver used, in this case cut to 10 chars col names:
     try:
         fi = grass.vector_db(map = vector)[int(layer)]
     except KeyError:
@@ -167,7 +164,7 @@ def main():
     # we need this for non-DBF driver:
     dbfdriver = fi['driver'] == 'dbf'
 
-    #Find out which table is linked to the vector map on the given layer
+    # Find out which table is linked to the vector map on the given layer
     if not fi['table']:
 	grass.fatal(_('There is no table connected to this map. Run v.db.connect or v.db.addtable first.'))
 
@@ -191,15 +188,15 @@ def main():
 
     addcols = []
     for i in basecols + extracols:
-	#check if column already present
+	# check if column already present
 	currcolumn = ("%s_%s" % (colprefix, i))
 	if dbfdriver:
 	    currcolumn = currcolumn[:10]
 
 	if currcolumn in grass.vector_columns(vector, layer).keys():
 	    if not flags['c']:
-		grass.fatal(("Cannot create column <%s> (already present). " % currcolumn) +
-			    "Use -c flag to update values in this column.")
+		grass.fatal((_("Cannot create column <%s> (already present). ") % currcolumn) +
+			    _("Use -c flag to update values in this column."))
 	else:
 	    if i == "n":
 		coltype = "INTEGER"
@@ -212,60 +209,76 @@ def main():
 	if grass.run_command('v.db.addcolumn', map = vector, columns = addcols) != 0:
 	    grass.fatal(_("Adding columns failed. Exiting."))
 
-    #loop over cats and calculate statistics:
+    # calculate statistics:
     grass.message(_("Processing data (%d categories)...") % number)
 
     # get rid of any earlier attempts
     grass.try_remove(sqltmp)
 
+    colnames = []
+    for var in basecols + extracols:
+	colname = '%s_%s' % (colprefix, var)
+	if dbfdriver:
+	    colname = colname[:10]
+	colnames.append(colname)
+
+    ntabcols = len(colnames)
+
     # do extended stats?
     if flags['e']:
 	extstat = 'e'
     else:
 	extstat = ""
-
+	
     f = file(sqltmp, 'w')
-    currnum = 1
-
-    for i in cats:
-	grass.message(_("Processing category %s (%d/%d)...") % (i, currnum, number))
-	grass.run_command('g.remove', rast = 'MASK', quiet = True, stderr = nuldev)
-	grass.mapcalc("MASK = if($name == $i, 1, null())",
-		      name = "%s_%s" % (vector, tmpname), i = i)
-
-        #n, min, max, range, mean, stddev, variance, coeff_var, sum
-	# How to test r.univar $? exit status? using -e creates the real possibility of out-of-memory errors
-	s = grass.read_command('r.univar', flags = 'g' + extstat, map = raster, percentile = percentile)
-	vars = grass.parse_key_val(s)
 
-	vars['cf_var'] = vars['coeff_var']
+    # do the stats
+    p = grass.pipe_command('r.univar', flags = 't' + 'g' + extstat, map = raster, 
+                      zones = rastertmp, percentile = percentile, fs = ';')
 
-	if flags['e'] and dbfdriver:
-	    percvar = 'percentile_' + percentile
-	    vars[perccol] = vars[percvar]
-
-	for var in basecols + extracols:
-	    value = vars[var]
-	    if value.lower() == 'nan':
+    first_line = 1
+    for line in p.stdout:
+	if first_line:
+	    first_line = 0
+	    continue
+
+	vars = line.rstrip('\r\n').split(';')
+
+	f.write("UPDATE %s SET" % fi['table'])
+	i = 2
+	first_var = 1
+	for colname in colnames:
+	    value = vars[i]
+	    # convert nan, +nan, -nan to NULL
+	    if value.lower().endswith('nan'):
 		value = 'NULL'
-	    colname = '%s_%s' % (colprefix, var)
-	    if dbfdriver:
-		colname = colname[:10]
-	    f.write("UPDATE %s SET %s=%s WHERE cat=%s;\n" % (fi['table'], colname, value, i))
+	    if not first_var:
+		f.write(" , ")
+	    else:
+		first_var = 0
+	    f.write(" %s=%s" % (colname, value))
+	    i += 1
+	    # skip n_null_cells, mean_of_abs, sum_of_abs
+	    if i == 3 or i == 8 or i == 13:
+		i += 1
 
-	currnum += 1
+	f.write(" WHERE %s=%s;\n" % (fi['key'], vars[0]))
 
+    p.wait()
     f.close()
 
     grass.message(_("Updating the database ..."))
     exitcode = grass.run_command('db.execute', input = sqltmp,
 				 database = fi['database'], driver = fi['driver'])
-    
+
     grass.run_command('g.remove', rast = 'MASK', quiet = True, stderr = nuldev)
 
     if exitcode == 0:
-	grass.message(("Statistics calculated from raster map <%s>" % raster) +
-		      (" and uploaded to attribute table of vector map <%s>." % vector))
+	grass.message((_("Statistics calculated from raster map <%s>") % raster) +
+		      (_(" and uploaded to attribute table of vector map <%s>.") % vector))
+    else:
+	grass.warning(_("Failed to upload statistics to attribute table of vector map <%s>.") % vector)
+    
     
     sys.exit(exitcode)