[GRASS-SVN] r48057 -
grass/branches/releasebranch_6_4/scripts/v.rast.stats
svn_grass at osgeo.org
svn_grass at osgeo.org
Fri Sep 2 04:52:46 EDT 2011
Author: mmetz
Date: 2011-09-02 01:52:46 -0700 (Fri, 02 Sep 2011)
New Revision: 48057
Modified:
grass/branches/releasebranch_6_4/scripts/v.rast.stats/v.rast.stats
Log:
v.rast.stats: update to grass-addons version
Modified: grass/branches/releasebranch_6_4/scripts/v.rast.stats/v.rast.stats
===================================================================
--- grass/branches/releasebranch_6_4/scripts/v.rast.stats/v.rast.stats 2011-09-01 20:20:43 UTC (rev 48056)
+++ grass/branches/releasebranch_6_4/scripts/v.rast.stats/v.rast.stats 2011-09-02 08:52:46 UTC (rev 48057)
@@ -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