[GRASS-SVN] r40385 - grass-addons/vector/v.rast.stats2

svn_grass at osgeo.org svn_grass at osgeo.org
Tue Jan 12 05:44:41 EST 2010


Author: dassau
Date: 2010-01-12 05:44:40 -0500 (Tue, 12 Jan 2010)
New Revision: 40385

Added:
   grass-addons/vector/v.rast.stats2/Makefile
   grass-addons/vector/v.rast.stats2/description.html
   grass-addons/vector/v.rast.stats2/v.rast.stats2
Log:
added new script

Added: grass-addons/vector/v.rast.stats2/Makefile
===================================================================
--- grass-addons/vector/v.rast.stats2/Makefile	                        (rev 0)
+++ grass-addons/vector/v.rast.stats2/Makefile	2010-01-12 10:44:40 UTC (rev 40385)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = v.rast.stats2
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script

Added: grass-addons/vector/v.rast.stats2/description.html
===================================================================
--- grass-addons/vector/v.rast.stats2/description.html	                        (rev 0)
+++ grass-addons/vector/v.rast.stats2/description.html	2010-01-12 10:44:40 UTC (rev 40385)
@@ -0,0 +1,78 @@
+<H2>DESCRIPTION</H2>
+
+<EM><b>v.rast.stats2</b></EM> - Calculates basic univariate statistics from
+a raster map only for the parts covered by the specified vector map.
+The vector map will be rasterized according to the raster map resolution.
+Then univariate statistics are calculated per vector category (cat) from
+the raster map and the results uploaded to the vector map attribute table.
+New columns are generated in the attribute table if not already present.
+<p>
+Nine columns are generated (n, min, max, range, mean, stddev, variance, 
+coeff_var, sum) according to the output of <em>r.univar2.zonal</em>.
+If the <B>-e</B> extended statistics flag is given the 1st quartile,
+median, 3rd quartile, and given percentile are also calculated.
+
+
+<H2>NOTES</H2>
+
+The module may take a long time to run if the raster region contains a large
+number of cells. In this case the <B>--verbose</B> flag may be used to track
+progress.
+<p>
+The script stops if a (prefixed) upload column is already present in the
+vector map attribute table, unless otherwise instructed with the <B>-c</B>
+continue flag. The column prefix will be separated from the statistic name
+with an underscore. For example with a prefix of "<tt>elev</tt>" the sum
+column will be named <tt>elev_sum</tt>.
+<P>
+If a DBF database is being used, note that column names are restricted by the
+DBF specification to 10 characters. Therefore it is advised to be economical
+in the use of the column prefix when using DBF as any additional characters
+will be chopped off.
+<p>
+If a MASK is present, it will be restored after the script finished.
+The script changes temporarily to the resolution of the given raster map.
+<P>
+<!-- r.univar limitation
+Large amounts of system memory can be used when the <B>-e</B> extended
+statistics flag is used with a very large region setting. If the region
+is too large the module should display memory allocation errors.
+Basic statistics can be calculated using any size input region.
+-->
+
+<H2>EXAMPLES</H2>
+
+Example to upload DEM statistics to vector field patches:
+
+<div class="code"><pre>
+# work on copy of original map:
+g.copy vect=fields,myfields
+# if needed, zoom to raster map:
+g.region rast=elevation.dem -p
+# calculate DEM statistics, upload to vector map table:
+v.rast.stats2 myfields raster=elevation.dem colprefix=dem
+# verify results:
+v.info -c myfields
+v.db.select myfields
+v.univar myfields column=dem_range type=centroid
+</pre></div>
+
+
+<H2>SEE ALSO</H2>
+
+<EM>
+<A HREF="r.univar.html">r.univar</A>,
+<A HREF="r.univar2.zonal.html">r.univar2.zonal</A>, 
+<A HREF="v.univar.html">v.univar</A>,
+<A HREF="v.rast.stats.html">v.rast.stats</A>,
+<a HREF="v.what.rast.html">v.what.rast</a>,
+<a href="v.what.vect.html">v.what.vect</a>
+</EM>
+
+<H2>AUTHOR</H2>
+
+Markus Neteler, CEA (<a href="http://www.eden-fp6project.net/">EDEN Project</a>)
+Otto Dassau
+
+<p>
+<i>Last changed: $Date: 2007-10-10 08:11:21 +0200 (Mi, 10. Okt 2007) $</i>

Added: grass-addons/vector/v.rast.stats2/v.rast.stats2
===================================================================
--- grass-addons/vector/v.rast.stats2/v.rast.stats2	                        (rev 0)
+++ grass-addons/vector/v.rast.stats2/v.rast.stats2	2010-01-12 10:44:40 UTC (rev 40385)
@@ -0,0 +1,337 @@
+#!/bin/sh
+
+############################################################################
+#
+# MODULE:	v.rast.stats2
+# 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, 2010 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
+#%  description: Calculates univariate statistics from a GRASS raster map based on vector polygons and uploads statistics to new attribute columns. Needs r.univar2.zonal from GRASS addons repository.
+#%  keywords: vector, raster, statistics
+#%End
+#%flag
+#%  key: c
+#%  description: Continue if upload column(s) already exist
+#%END
+#%flag
+#%  key: e
+#%  description: Calculate extended statistics
+#%END
+#%option
+#% key: vector
+#% type: string
+#% key_desc: name
+#% gisprompt: old,vector,vector
+#% description: Name of vector polygon map
+#% required : yes
+#%End
+#%option
+#% key: layer
+#% type: integer
+#% description: Layer to which the table to be changed is connected
+#% answer: 1
+#% required : no
+#%end
+#%option
+#% key: keycol
+#% type: string
+#% description: Key column of attribute table used to join statistics 
+#% answer: cat
+#% required : no
+#%end
+#%option
+#% key: raster
+#% type: string
+#% key_desc: name
+#% gisprompt: old,cell,raster
+#% description: Name of raster map to calculate statistics from
+#% required : yes
+#%END
+#%option
+#% key: colprefix
+#% type: string
+#% description: Column prefix for new attribute columns
+#% required : yes
+#%end
+#%option
+#% key: percentile
+#% type: integer
+#% description: Percentile to calculate (requires extended statistics flag)
+#% options: 0-100
+#% answer: 90
+#% required : no
+#%end
+
+if  [ -z "$GISBASE" ] ; then
+    echo "You must be in GRASS GIS to run this program." >&2
+    exit 1
+fi
+
+if [ "$1" != "@ARGS_PARSED@" ] ; then
+    exec g.parser "$0" "$@"
+fi
+
+PROG=`basename $0`
+
+#### check if we have awk
+if [ ! -x "`which awk`" ] ; then
+    g.message "awk required, please install awk or gawk first" 
+    exit 1
+fi
+
+# setting environment, so that awk works properly in all languages
+unset LC_ALL
+LC_NUMERIC=C
+export LC_NUMERIC
+
+#### setup temporary file
+TMP="`g.tempfile pid=$$`"
+if [ $? -ne 0 ] || [ -z "$TMP" ] ; then
+    g.message -e "Unable to create temporary files" 
+    exit 1
+fi
+
+### variables for temp files
+SQLTMP="$TMP.sql"
+STATSTMP="$TMP.csv"
+COLNAMETMP="$TMP.col"
+
+### we need a random name
+TMPNAME=`basename "$TMP"`
+
+cleanup()
+{
+   #restore settings:
+   g.region region="$TMPNAME"
+   g.remove region="$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..."
+      g.rename "${TMPNAME}_origmask",MASK --quiet
+   fi
+   rm -f "$TMP" "$TMPNAME" "$TMP.cats" "$SQLTMP"
+}
+
+# what to do in case of user break:
+exitprocedure()
+{
+   g.message -e 'User break!'
+   cleanup
+   exit 1
+}
+# shell check for user break (signal list: trap -l)
+trap "exitprocedure" 2 3 15
+
+RASTER="$GIS_OPT_RASTER"
+COLPREFIX="$GIS_OPT_COLPREFIX"
+KEYCOL="$GIS_OPT_KEYCOL"
+
+### setup enviro vars ###
+eval `g.gisenv`
+: ${GISBASE?} ${GISDBASE?} ${LOCATION_NAME?} ${MAPSET?}
+
+# does map exist in CURRENT mapset?
+eval `g.findfile element=vector file="$GIS_OPT_VECTOR" mapset="$MAPSET"`
+VECT_MAPSET=`echo "$GIS_OPT_VECTOR" | grep '@' | cut -f2 -d'@'`
+if [ -z "$VECT_MAPSET" ] ; then
+   VECT_MAPSET="$MAPSET"
+fi
+if [ ! "$file" ] || [ "$VECT_MAPSET" != "$MAPSET" ] ; then
+   g.message -e "Vector map <$GIS_OPT_VECTOR> not found in current mapset"
+   exit 1
+else
+   VECTOR=`echo "$GIS_OPT_VECTOR" | cut -f1 -d'@'`
+fi
+
+#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
+MASKFOUND=0
+eval `g.findfile element=cell file=MASK`
+if [ "$file" ] ; then
+   g.message "Raster MASK found, temporarily disabled"
+   g.rename MASK,"${TMPNAME}_origmask" --quiet
+   MASKFOUND=1
+fi
+
+#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."
+   cleanup
+   exit 1
+fi
+EWRES=`r.info -s "$RASTER" | grep ewres | cut -d'=' -f2`
+
+#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
+
+# create raster from vector map for r.univar.zonal 
+v.to.rast in="$VECTOR" out="${VECTOR}_${TMPNAME}" use=attr \
+col="$KEYCOL" 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:
+r.category "${VECTOR}_${TMPNAME}" fs=';' --quiet | cut -d';' -f1 > "$TMP.cats"
+#echo "List of categories found: $CATSLIST"
+NUMBER=`cat "$TMP.cats" | wc -l | awk '{print $1}'`
+if [ $NUMBER -lt 1 ] ; then
+   g.message -e "No categories found in raster map"
+   cleanup
+   exit 1
+fi
+
+#check if DBF driver used, in this case cut to 10 chars col names:
+DBFDRIVER=0
+v.db.connect -g "$VECTOR" fs=";" | grep -w "$GIS_OPT_LAYER" | cut -d';' -f5 | grep -i dbf --quiet
+if [ $? -eq 0 ] ; then
+   DBFDRIVER=1
+else
+   # in case a table is missing, we'll trap a crash later...
+   DBFDRIVER=0
+fi
+# we need this for non-DBF driver:
+DB_SQLDRIVER=`v.db.connect -g "$VECTOR" fs=";" | cut -d';' -f5`
+DB_DATABASE="`v.db.connect -g "$VECTOR" fs=";" | cut -d';' -f4`"
+
+#Find out which table is linked to the vector map on the given layer
+TABLE=`v.db.connect map="$VECTOR" -g fs=";" | grep -w "$GIS_OPT_LAYER" | awk -F ";" '{print $2}'`
+if [ -z "$TABLE" ] ; then
+   g.message -e 'There is no table connected to this map! Run v.db.connect or v.db.addtable first.'
+   exit 1
+fi
+
+BASECOLS="n min max range mean stddev variance cf_var sum"
+
+# do extended stats?
+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"
+   else
+      PERCCOL="percentile_$GIS_OPT_PERCENTILE"
+   fi
+   EXTRACOLS="first_quartile median third_quartile $PERCCOL"
+else
+   EXTRACOLS=""
+fi
+
+unset ADDCOLS
+for i in $BASECOLS $EXTRACOLS ; do
+  #check if column already present
+  if [ $DBFDRIVER -eq 1 ] ; then
+     CURRCOLUMN="`echo "${COLPREFIX}_${i}" | cut -b1-10`"
+  else
+     CURRCOLUMN="${COLPREFIX}_${i}"
+  fi
+  v.info -c $VECTOR --quiet | sed 's+^+|+g' | sed 's+$+|+g' | \
+    grep "|$CURRCOLUMN|" --quiet
+  if [ $? -eq 0 ] ; then
+    if [ $GIS_FLAG_C -ne 1 ] ; then
+       g.message -e "Cannot create column <$CURRCOLUMN> (already present). \
+         Use -c flag to update values in this column."
+       cleanup
+       exit 1
+    fi
+  else
+    if [ -n "$ADDCOLS" ] ; then
+       ADDCOLS="$ADDCOLS,"
+    fi
+    if [ "$i" = "n" ] ; then
+       COLTYPE="INTEGER"
+    else
+       COLTYPE="DOUBLE PRECISION"
+    fi
+    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
+fi
+
+# calculate statistics for zones:
+g.message -v "Processing data ..."
+
+# get rid of any earlier attempts
+rm -f "$SQLTMP" "$STATSTMP" "$COLNAMETMP"
+
+# Processing statistics 
+unset $BASECOLS $EXTRACOLS
+if [ $GIS_FLAG_E -eq 1 ] ; then 
+    r.univar.zonal -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
+    r.univar.zonal -t map="$RASTER" zones="${VECTOR}_${TMPNAME}" | \
+    cut -f1,3,5-8,10-13 -d'|' | sed 's+nan+NULL+g' > "$STATSTMP"
+fi
+
+if [ $GIS_FLAG_E -eq 1 ] && [ $DBFDRIVER -eq 1 ] ; then
+    eval "$PERCCOL=\$percentile_$GIS_OPT_PERCENTILE"
+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
+    echo "$colname" >> "$COLNAMETMP" 
+done
+
+# create array with new column names
+col=( `cat "$COLNAMETMP" | tr '\n' ' '` )
+
+# 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 '${VECTOR}' SET '${col[0]}' = %i , '${col[1]}' = %.2f , '${col[2]}' = %.2f , '${col[3]}' = %.2f , '${col[4]}' = %.2f , '${col[5]}' = %.2f , '${col[6]}' = %.2f , '${col[7]}' = %.2f , '${col[8]}' = %.2f , '${col[9]}' = %.2f , '${col[10]}' = %.2f , '${col[11]}' = %2f , '${col[12]}' = %.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 '${VECTOR}' SET '${col[0]}' = %i , '${col[1]}' = %.2f , '${col[2]}' = %.2f , '${col[3]}' = %.2f , '${col[4]}' = %.2f , '${col[5]}' = %.2f , '${col[6]}' = %.2f , '${col[7]}' = %.2f , '${col[8]}' = %.2f WHERE '${KEYCOL}' = %i;", $2,$3,$4,$5,$6,$7,$8,$9,$10,$1}' > "$SQLTMP"   
+fi
+
+g.message -v "Updating the database ..."
+db.execute input="$SQLTMP" database="${DB_DATABASE}" driver="$DB_SQLDRIVER"
+EXITCODE=$?
+
+cleanup
+
+if [ "$EXITCODE" -eq 0 ] ; then
+  g.message "Statistics calculated from raster map <$RASTER> and uploaded to \
+    attribute table of vector map <$VECTOR>."
+  g.message "Done."
+fi
+
+exit $EXITCODE


Property changes on: grass-addons/vector/v.rast.stats2/v.rast.stats2
___________________________________________________________________
Added: svn:executable
   + *



More information about the grass-commit mailing list