[GRASS-user] please help: v.rast.stats

[Nikos Alexandris]

On Tue, 2008-03-04 at 12:57 +0000, christian Brandt wrote:
> Dear list,

Hi Christian
> 
> I am still trying to get statistic values of raster cells overlaid by
> a polygon theme using the v.rast.stats command.
> 
> Executing this commands doesn't result in the real statistics within
> one cell but in any virtual figures. Also the statistic attributes
> joined to the vector layer doesn't seem to match the cell content.

Not sure but:

did you check the resolution (g.region -p) ?

> 
> Any suggestions?

Once I needed to do something similar and I got the ideas on how to go
about it from the list of course... 

Here is what I finally did to grab stats of raster maps inside a vector
(polygon) map. Hope it in some way useful!

(I know I am not the best scripter... ! You probably need to adjust
everything to fit your needs. Pay attention to the parameters of
"r.stats" command... or maybe you need the r.univar )


#!/bin/sh
# List pixel values covered by areas of interest for each area
separately
# First attempt, Nikos Alexandris

if  [ -z "$GISBASE" ] ; then
    echo "You must be in GRASS GIS to run this program." >&2
    exit 1
fi

PROG="$0"

echo $0

if [ $# -lt 1 -o "$1" = "-h" -o "$1" = "--help" -o "$1" = "-help" ] ; 
then
 echo "\n* Please provide vector areas of interest and four (4) raster
sources!"
 echo "Usage:"
 echo "        $PROG [Vector areas of interest] [Raster source 1]
[Raster source 2] [Raster source 3] [Raster source 4]\n"
 exit 1
fi

echo "Script to list all unique top quality pixel values in raster(s)
covered by vector areas of interest"
echo "Requires top quality pixels MASK!"

echo "\nD i d    y o u    c h e c k    t h e    r e s o l u t i o n ?"
g.region -p

sleep 3

echo "\nExplanations:"
echo "a. Repeated value(s) listed as many times as present"
echo "b. Areas derived from vector polygons"
echo "c. Assuming first column is "cat""
sleep 1

# Source for vector area(s) of interest
AREASVECT=$1 
echo "\nVector source is: $1"

# Raster(s) from which pixel values will be extracted
RASTERSRC1=$2
echo "\nRaster source 1 is: $2"

RASTERSRC2=$3
echo "\nRaster source 2 is: $3"

RASTERSRC3=$4
echo "\nRaster source 3 is: $4"

RASTERSRC4=$5
echo "\nRaster source 4 is: $5"
sleep 2

echo "\nThis is year... ?"
read YEAR

# Define the number of areas based on "cat" column of AREASVECT
CATi="`db.select table=$AREASVECT | grep -v cat | cut -d"|" -f1 `"
for i in $CATi; do
	CATMAX=$(echo $i);
done
echo "Table [$AREASVECT] contains $CATMAX areas."
sleep 3

# Step 1. Rasterize vector polygons and assign cat values (one unique
value) to each raster area
echo "\nRasterising vector polygons"
v.to.rast input=$AREASVECT output=AREASRAST use=cat type=area layer=1
--overwrite
sleep 3

# Step 2. MASK-ing based on raster area(s) of interest and grabbing
pixel values through loop 
echo "\nLooping over all areas of interest to grab pixel values:"
rm -rf ALL_PV_$YEAR.*.csv

# where "i" the "cat" value derived from VECTOR1
for i in $CATi; do
	echo "\n Processing Area [$i] out of [$CATMAX]"
	# Step 2(a). Masking areas of interest
	echo "* Masking..."
	r.mapcalc MASK="if(AREASRAST==$i, 1, null())"
	echo "(Keep copy of MASK"$i")"
	g.copy rast=MASK,MASK.$AREASVECT.$YEAR.$i
	# Step 2(b). Extract all unique pixel values for above created MASK
along with their x,y and grid coordinates
	echo "* Grabbing pixel values... (exporting in: ALL_PV_"$i".csv)"
	r.stats -nxg input=$RASTERSRC1,$RASTERSRC2,$RASTERSRC3,$RASTERSRC4
fs="," output=$AREASVECT.ALL_PV_$YEAR.$i.csv
# If a nicely formatted output is desired, pipe the output into a
command which can create columnar output.
# For example, the command:
# r.stats input=a,b,c | pr -3 | cat -s
	cat $AREASVECT.ALL_PV_$YEAR.$i.csv >> $AREASVECT.ALL_PV_$YEAR.csv
	cat $AREASVECT.ALL_PV_$YEAR.$i.csv | cut -d "," -f5,6,7,8 >>
$AREASVECT.ALL_PV_$YEAR.NoXY.csv
	g.remove MASK --quiet && echo "MASK"$i" removed... proceeding to next
one!\n"
done

g.message "Cleaning temporarily created files..."
g.message "\n R E N A M E  output csv file before using script again!!!"
g.remove rast=AREASRAST --quiet

echo "\nDone!"
echo "Order of output: grid coordinates, pixel column and row, pixel
values $RASTERSRC1, $RASTERSRC2, $RASTERSRC3, $RASTERSRC4"

# Step 3. Add extracted pixel values to table of "VECTOR1"
# loop through all of the shapefiles in the directory and load them

> 
> Regards
> 
> Chris

Cheers,

Nikos