[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