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

[Dylan Beaudette]

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

Also-- see the application 'starspan'. works with GRASS, based on GDAL and 
GEOS.

Cheers,

Dylan

>
> 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
>
> _______________________________________________
> grass-user mailing list
> grass-user at lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user



-- 
Dylan Beaudette
Soil Resource Laboratory
http://casoilresource.lawr.ucdavis.edu/
University of California at Davis
530.754.7341