[GRASS-SVN] r48058 - grass/branches/develbranch_6/scripts/v.rast.stats

svn_grass at osgeo.org svn_grass at osgeo.org
Fri Sep 2 05:02:58 EDT 2011


Author: mmetz
Date: 2011-09-02 02:02:58 -0700 (Fri, 02 Sep 2011)
New Revision: 48058

Modified:
   grass/branches/develbranch_6/scripts/v.rast.stats/v.rast.stats
Log:
v.rast.stats: update to grass-addons version

Modified: grass/branches/develbranch_6/scripts/v.rast.stats/v.rast.stats
===================================================================
--- grass/branches/develbranch_6/scripts/v.rast.stats/v.rast.stats	2011-09-02 08:52:46 UTC (rev 48057)
+++ grass/branches/develbranch_6/scripts/v.rast.stats/v.rast.stats	2011-09-02 09:02:58 UTC (rev 48058)
@@ -3,16 +3,15 @@
 ############################################################################
 #
 # MODULE:	v.rast.stats
-# AUTHOR(S):	Markus Neteler
+# AUTHOR(S):	Markus Neteler, Otto Dassau
 # PURPOSE:	Calculates univariate statistics from a GRASS raster map
 #		only for areas covered by vector objects on a per-category base
-# COPYRIGHT:	(C) 2005 by the GRASS Development Team
+# COPYRIGHT:	(C) 2005, 2011 by the GRASS Development Team
 #
 #		This program is free software under the GNU General Public
 #		License (>=v2). Read the file COPYING that comes with GRASS
 #		for details.
 #
-# TODO: do we need layer= management?
 #############################################################################
 
 #%Module
@@ -38,7 +37,7 @@
 #%option
 #% key: layer
 #% type: integer
-#% label: Layer number
+#% gisprompt: old_layer,layer,layer
 #% description: A single vector map can be connected to multiple database tables. This number determines which table to use
 #% answer: 1
 #% required : no
@@ -94,16 +93,21 @@
     g.message -e "Unable to create temporary files" 
     exit 1
 fi
+
+### variables for temp files
 SQLTMP="$TEMPFILE.sql"
+STATSTMP="$TEMPFILE.csv"
+COLNAMETMP="$TEMPFILE.col"
+
 # we need a random name
 TMPNAME=`basename "$TEMPFILE"`
 
 cleanup()
 {
-   #restore settings:
+   # restore settings:
    g.region region="$TMPNAME"
    g.remove region="$TMPNAME" --quiet
-   g.remove rast="${VECTOR}_${TMPNAME}" --quiet
+   g.remove rast="${VECTOR}_$TMPNAME" --quiet
    g.remove MASK --quiet 2>/dev/null
    if [ $MASKFOUND -eq 1 ] ; then
       g.message "Restoring previous MASK..."
@@ -122,7 +126,6 @@
 # shell check for user break (signal list: trap -l)
 trap "exitprocedure" 2 3 15
 
-
 RASTER="$GIS_OPT_RASTER"
 COLPREFIX="$GIS_OPT_COLPREFIX"
 
@@ -141,16 +144,17 @@
    exit 1
 else
    VECTOR=`echo "$GIS_OPT_VECTOR" | cut -f1 -d'@'`
+   VECTORFULL="${VECTOR}@${VECT_MAPSET}"
 fi
 
-#check the input raster map
+# check the input raster map
 eval `g.findfile element=cell file="$RASTER"`
 if [ ! "$file" ] ; then
    g.message -e "Raster map <$RASTER> not found"
    exit 1
 fi
 
-#check presence of raster MASK, put it aside
+# check presence of raster MASK, put it aside
 MASKFOUND=0
 eval `g.findfile element=cell file=MASK`
 if [ "$file" ] ; then
@@ -159,8 +163,8 @@
    MASKFOUND=1
 fi
 
-#get RASTER resolution of map which we want to query:
-#fetch separated to permit for non-square cells (latlong etc)
+# get RASTER resolution of map which we want to query:
+# fetch separated to permit for non-square cells (latlong etc)
 NSRES=`r.info -s "$RASTER" | grep nsres | cut -d'=' -f2`
 if [ $? -ne 0 ] ; then
    g.message -e "An error occurred reading the input raster map resolution."
@@ -169,25 +173,23 @@
 fi
 EWRES=`r.info -s "$RASTER" | grep ewres | cut -d'=' -f2`
 
-#save current settings:
+# save current settings:
 g.region save="$TMPNAME" --quiet
 
-#Temporarily setting raster resolution to $RASTER resolution
-#keep boundary settings
-g.region nsres="$NSRES" ewres="$EWRES" -a
+# temporarily aligning current region to raster <$RASTER>
+g.region align="$RASTER"
 
-#prepare raster MASK
-v.to.rast in="$VECTOR" out="${VECTOR}_$TMPNAME" use=cat --quiet
+# create raster from vector map for r.univar 
+v.to.rast in="$VECTORFULL" out="${VECTOR}_$TMPNAME" layer="$GIS_OPT_LAYER" use=cat type=area --quiet
 if [ $? -ne 0 ] ; then
    g.message -e "An error occurred while converting vector to raster"
    cleanup
    exit 1
 fi
 
-
-#dump cats to file to avoid "too many argument" problem:
+# dump cats to file to avoid "too many argument" problem:
 r.category "${VECTOR}_$TMPNAME" fs=';' --quiet | cut -d';' -f1 > "$TEMPFILE.cats"
-#echo "List of categories found: $CATSLIST"
+# echo "List of categories found: $CATSLIST"
 NUMBER=`cat "$TEMPFILE.cats" | wc -l | awk '{print $1}'`
 if [ $NUMBER -lt 1 ] ; then
    g.message -e "No categories found in raster map"
@@ -195,11 +197,10 @@
    exit 1
 fi
 
-
-#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:
 DBFDRIVER=0
 v.db.connect -gl map="$VECTOR" fs="|" layer="$GIS_OPT_LAYER" | \
-   cut -d'|' -f5 | grep -i 'dbf' -q
+   cut -d'|' -f5 | grep -i 'dbf' --quiet
 if [ $? -eq 0 ] ; then
    DBFDRIVER=1
 else
@@ -207,9 +208,9 @@
    DBFDRIVER=0
 fi
 
-if [ "$DBFDRIVER" -eq 1 ] ; then
+if [ $DBFDRIVER -eq 1 ] ; then
    # `wc -c` reports number of chars + 1 for the newline
-   if [ `echo "$COLPREFIX" | wc -c` -gt 7 ] ; then
+   if [ `echo "$COLPREFIX" | wc -c` -gt 10 ] ; then
        g.message -e "Cannot create unique names for columns. \
          Either use a shorter column prefix or switch to another DB \
 	 backend such as SQLite. DBF limits the length to 10 characters"
@@ -225,7 +226,7 @@
 DB_DATABASE="`v.db.connect -gl map=$VECTOR fs='|' layer=$GIS_OPT_LAYER | cut -d'|' -f4`"
 g.message -d message="database: [$DB_DATABASE]"
 
-#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
 TABLE=`v.db.connect map="$VECTOR" -gl fs='|' layer="$GIS_OPT_LAYER" | cut -d'|' -f2`
 if [ -z "$TABLE" ] ; then
    g.message -e 'There is no table connected to this map! Run v.db.connect or v.db.addtable first.'
@@ -233,11 +234,13 @@
 fi
 g.message -d message="table: [$TABLE]"
 
+# get key column
+KEYCOL=`v.db.connect map="$VECTOR" -gl fs='|' layer="$GIS_OPT_LAYER" | cut -d'|' -f3`
 
 BASECOLS="n min max range mean stddev variance cf_var sum"
 
 # do extended stats?
-if [ "$GIS_FLAG_E" -eq 1 ] ; then
+if [ $GIS_FLAG_E -eq 1 ] ; then
    # namespace is limited in DBF but the % value is important
    if [ "$DBFDRIVER" -eq 1 ] ; then
       PERCCOL="per$GIS_OPT_PERCENTILE"
@@ -249,18 +252,23 @@
    EXTRACOLS=""
 fi
 
-
 unset ADDCOLS
+unset COLNAMES
 for i in $BASECOLS $EXTRACOLS ; do
-  #check if column already present
-  if [ "$DBFDRIVER" -eq 1 ] ; then
+  # check if column already present
+  if [ $DBFDRIVER -eq 1 ] ; then
      CURRCOLUMN="`echo "${COLPREFIX}_$i" | cut -b1-10`"
   else
      CURRCOLUMN="${COLPREFIX}_$i"
   fi
+  if [ -n "$COLNAMES" ] ; then
+     COLNAMES="${COLNAMES},$CURRCOLUMN"
+  else
+     COLNAMES="$CURRCOLUMN"
+  fi
   g.message -d message="Looking for <$CURRCOLUMN> in <$VECTOR> ..."
-  v.info -c "$VECTOR" --quiet | sed -e 's+^+|+' -e 's+$+|+' | \
-    grep "|$CURRCOLUMN|" -q
+  v.info -c "$VECTORFULL" layer="$GIS_OPT_LAYER" --quiet | sed -e 's+^+|+' -e 's+$+|+' | \
+    grep "|$CURRCOLUMN|" --quiet
   if [ $? -eq 0 ] ; then
     if [ "$GIS_FLAG_C" -ne 1 ] ; then
        g.message -e "Cannot create column <$CURRCOLUMN> (already present). \
@@ -277,72 +285,86 @@
     else
        COLTYPE="DOUBLE PRECISION"
     fi
-    ADDCOLS="$ADDCOLS $CURRCOLUMN $COLTYPE"
+    ADDCOLS="${ADDCOLS}$CURRCOLUMN $COLTYPE"
   fi
 done
 
-g.message -v "Adding columns '$ADDCOLS'"
-v.db.addcol map="$VECTOR" columns="$ADDCOLS"
-if [ $? -ne 0 ] ; then
-   g.message -e "Cannot continue (problem adding columns)."
-   cleanup
-   exit 1
+if [ -n "$ADDCOLS" ] ; then
+   g.message -v "Adding columns <$ADDCOLS>"
+#   v.db.addcol map="$VECTORFULL" columns="$ADDCOLS" layer="$GIS_OPT_LAYER"
+#   if [ $? -ne 0 ] ; then
+#      g.message -e "Cannot continue (problem adding columns)."
+#      cleanup
+#      exit 1
+#   fi
+
+  # code borrowed from v.db.addcol
+  colnum=`echo "$ADDCOLS" | awk -F, '{print NF}'`
+
+  n=1
+  while [ "$n" -le "$colnum" ]
+  do
+    col=`echo "$ADDCOLS" | cut -d',' -f$n`
+    
+    echo "ALTER TABLE $TABLE ADD COLUMN $col" | db.execute database="${DB_DATABASE}" driver="${DB_SQLDRIVER}"
+    if [ $? -eq 1 ] ; then
+      g.message -e "Cannot continue (problem adding column)."
+      exit 1
+    fi
+    n=`expr $n + 1`
+  done
 fi
 
+# calculate statistics for zones:
+g.message -v "Processing data ..."
 
-#loop over cats and calculate statistics:
-g.message "Processing data ($NUMBER categories)..."
-CURRNUM=1
 # get rid of any earlier attempts
-rm -f "$SQLTMP"
+rm -f "$SQLTMP" "$STATSTMP" "$COLNAMETMP"
 
-# do extended stats?
-if [ $GIS_FLAG_E -eq 1 ] ; then
-   EXTSTAT="-e"
+# Processing statistics 
+unset $BASECOLS $EXTRACOLS
+if [ $GIS_FLAG_E -eq 1 ] ; then 
+    r.univar -t -e map="$RASTER" \
+    zones="${VECTOR}_$TMPNAME" percentile="$GIS_OPT_PERCENTILE" | \
+    cut -f1,3,5-8,10-13,15-18 -d'|' | sed 's+nan+NULL+g' > "$STATSTMP"
 else
-   EXTSTAT=""
+    r.univar -t map="$RASTER" zones="${VECTOR}_$TMPNAME" | \
+    cut -f1,3,5-8,10-13 -d'|' | sed 's+nan+NULL+g' > "$STATSTMP"
 fi
 
-for i in `cat "$TEMPFILE.cats"` ; do
-   g.message -v message="Processing category $i ($CURRNUM/$NUMBER)"
-   g.remove MASK --quiet 2>/dev/null
-   # r.mapcalc doesn't use the parser, so we set the env var version instead
-   GRASS_VERBOSE=0 \
-     r.mapcalc "MASK=if(${VECTOR}_${TMPNAME} == $i, 1, null())"
+if [ $GIS_FLAG_E -eq 1 ] && [ $DBFDRIVER -eq 1 ] ; then
+    eval "$PERCCOL=\$percentile_$GIS_OPT_PERCENTILE"
+fi
 
-   #n, min, max, range, mean, stddev, variance, coeff_var, sum
-   unset $BASECOLS $EXTRACOLS
-   eval `r.univar -g $EXTSTAT map="$RASTER" percentile="$GIS_OPT_PERCENTILE" | sed 's+nan+NULL+g'`
-   # How to test r.univar $? exit status? using -e creates the real possibility of out-of-memory errors
+# create array with new column names
+col1="`echo $COLNAMES | cut -f1 -d','`"
+col2="`echo $COLNAMES | cut -f2 -d','`"
+col3="`echo $COLNAMES | cut -f3 -d','`"
+col4="`echo $COLNAMES | cut -f4 -d','`"
+col5="`echo $COLNAMES | cut -f5 -d','`"
+col6="`echo $COLNAMES | cut -f6 -d','`"
+col7="`echo $COLNAMES | cut -f7 -d','`"
+col8="`echo $COLNAMES | cut -f8 -d','`"
+col9="`echo $COLNAMES | cut -f9 -d','`"
+if [ $GIS_FLAG_E -eq 1 ] ; then
+   col10="`echo $COLNAMES | cut -f10 -d','`"
+   col11="`echo $COLNAMES | cut -f11 -d','`"
+   col12="`echo $COLNAMES | cut -f12 -d','`"
+   col13="`echo $COLNAMES | cut -f13 -d','`"
+fi
 
-   cf_var="$coeff_var"
-   if [ $GIS_FLAG_E -eq 1 ] && [ $DBFDRIVER -eq 1 ] ; then
-      eval "$PERCCOL=\$percentile_$GIS_OPT_PERCENTILE"
-   fi
+# create SQL file for extended and normal statistics
+g.message -v "Creating SQL file ..."
+if [ $GIS_FLAG_E -eq 1 ] ; then
+    sed -e '1d' "$STATSTMP" | awk -F "|" '{printf "\nUPDATE '${TABLE}' SET '${col1}' = %i , '${col2}' = %.2f , '${col3}' = %.2f , '${col4}' = %.2f , '${col5}' = %.2f , '${col6}' = %.2f , '${col7}' = %.2f , '${col8}' = %.2f , '${col9}' = %.2f , '${col10}' = %.2f , '${col11}' = %.2f , '${col12}' = %2f , '${col13}' = %.2f WHERE '${KEYCOL}' = %i;", $2,$3,$4,$5,$6,$7,$8,$9,$10,$11,$12,$13,$14,$1}' > "$SQLTMP"
+else
+    sed -e '1d' "$STATSTMP" | awk -F "|" '{printf "\nUPDATE '${TABLE}' SET '${col1}' = %i , '${col2}' = %.2f , '${col3}' = %.2f , '${col4}' = %.2f , '${col5}' = %.2f , '${col6}' = %.2f , '${col7}' = %.2f , '${col8}' = %.2f , '${col9}' = %.2f WHERE '${KEYCOL}' = %i;", $2,$3,$4,$5,$6,$7,$8,$9,$10,$1}' > "$SQLTMP"   
+fi
 
-   for var in $BASECOLS $EXTRACOLS ; do
-      eval value=\${$var}
-
-      if [ $DBFDRIVER -eq 1 ] ; then
-	colname="`echo "${COLPREFIX}_$var" | cut -b1-10`"
-      else
-	colname="${COLPREFIX}_$var"
-      fi
-      g.message -d debug=3 message="UPDATE $TABLE SET ${colname}=${value} WHERE cat=$i;"
-      echo "UPDATE $TABLE SET ${colname}=${value} WHERE cat=$i;" >> "$SQLTMP"
-   done
-
-   CURRNUM=`expr $CURRNUM + 1`
-done
-
-
-g.message "Updating the database ..."
-db.execute input="$SQLTMP" database="$DB_DATABASE" driver="$DB_SQLDRIVER"
+g.message message="Updating the database ..."
+db.execute input="$SQLTMP" database="${DB_DATABASE}" driver="$DB_SQLDRIVER"
 EXITCODE=$?
 
-
-g.remove MASK --quiet 2>/dev/null
-
 cleanup
 
 if [ "$EXITCODE" -eq 0 ] ; then



More information about the grass-commit mailing list