[GRASS-SVN] r44444 - grass/trunk/scripts/v.rast.stats
svn_grass at osgeo.org
svn_grass at osgeo.org
Sat Nov 27 03:16:04 EST 2010
Author: mmetz
Date: 2010-11-27 00:16:04 -0800 (Sat, 27 Nov 2010)
New Revision: 44444
Modified:
grass/trunk/scripts/v.rast.stats/v.rast.stats.py
Log:
fix region alignment, table key column, -nan, enable translations, and speed up the module
Modified: grass/trunk/scripts/v.rast.stats/v.rast.stats.py
===================================================================
--- grass/trunk/scripts/v.rast.stats/v.rast.stats.py 2010-11-27 07:38:16 UTC (rev 44443)
+++ grass/trunk/scripts/v.rast.stats/v.rast.stats.py 2010-11-27 08:16:04 UTC (rev 44444)
@@ -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 @@
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 @@
# 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 @@
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 @@
# 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 @@
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,61 +209,77 @@
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)
+ # do the stats
+ p = grass.pipe_command('r.univar', flags = 't' + 'g' + extstat, map = raster,
+ zones = rastertmp, percentile = percentile, fs = ';')
- #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)
+ first_line = 1
+ for line in p.stdout:
+ if first_line:
+ first_line = 0
+ continue
- vars['cf_var'] = vars['coeff_var']
+ vars = line.rstrip('\r\n').split(';')
- if flags['e'] and dbfdriver:
- percvar = 'percentile_' + percentile
- vars[perccol] = vars[percvar]
-
- for var in basecols + extracols:
- value = vars[var]
- if value.lower() == 'nan':
+ 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)
if __name__ == "__main__":
More information about the grass-commit
mailing list