Changeset 44444

Timestamp:

Nov 27, 2010, 12:16:04 AM (14 years ago)

Author:

mmetz

Message:

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

File:

• grass/trunk/scripts/v.rast.stats/v.rast.stats.py (modified) (12 diffs)

grass/trunk/scripts/v.rast.stats/v.rast.stats.py

@@ -4 +4 @@
 #
 # 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
@@ -79 +81 @@
 
 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:
@@ -88 +90 @@
 
 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
@@ -120 +122 @@
     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:
@@ -130 +134 @@
         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)
-
-    #prepare raster MASK
-    if grass.run_command('v.to.rast', input = vector, output = "%s_%s" % (vector, tmpname),
+    # 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 = 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:
@@ -155 +153 @@
     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)]
@@ -168 +165 @@
     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.'))
@@ -192 +189 @@
     addcols = []
     for i in basecols + extracols:
-        #check if column already present
+        # check if column already present
         currcolumn = ("%s_%s" % (colprefix, i))
         if dbfdriver:
@@ -199 +196 @@
         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":
@@ -213 +210 @@
             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?
@@ -224 +230 @@
     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']
-
-        if flags['e'] and dbfdriver:
-            percvar = 'percentile_' + percentile
-            vars[perccol] = vars[percvar]
-
-        for var in basecols + extracols:
-            value = vars[var]
-            if value.lower() == 'nan':
+
+    # do the stats
+    p = grass.pipe_command('r.univar', flags = 't' + 'g' + extstat, map = raster, 
+                      zones = rastertmp, percentile = percentile, fs = ';')
+
+    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))
-
-        currnum += 1
-
+            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
+
+        f.write(" WHERE %s=%s;\n" % (fi['key'], vars[0]))
+
+    p.wait()
     f.close()
 
@@ -261 +271 @@
     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))
+    else:
+        grass.warning(_("Failed to upload statistics to attribute table of vector map <%s>.") % vector)
     
-    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))
     
     sys.exit(exitcode)