[GRASS-user] How to list all unique pixel values (of a raster) covered by polygons (areas) for each polygon separately(?)

[nikos.alexandris at felis.uni-freiburg.de]

Here it goes... :

It's my first attempt to write a functional script
(attached below), so please, if you have the time, correct
my mistakes.

So far:

1. I managed to "grab" pixel values from rasters covered by
vector areas and export them in .txt (or .csv) file(s).


2. If I get the "upload to vector table" task working
(point 3.), I would like to make the script more flexible
so it can automatically run the process for more raster(s).

I mean: grab pix values from all files named after pattern
"MOD200????" -- e.g. MOD2006NIR, MOD2006R, MOD2006Blue,
MOD2007NIR, etc. --

3. I am stuck with how to get these on a table (other than
the existing one for my source vector polygons)


I want to create a new table where I upload the extracted
pixel values.


[ I think I got the meaning of layers (tables connected to
on e vector) but still, I can't:

a. update in a new "layer2", a column "areacat" with the
"cat" values from "layer1" (?)

b. (if "a" successful then put in the loop a v.db.update
command) 

About "a" following command: > v.to.db map=$AREASVECT
type=boundary layer=2 qlayer=1 option=query
column='areacat' qcolumn='cat'

gives:

Reading data from the map...
 100%
Querying database...
 100%
Updating database...
 100%
1 categories read from map
0 records selected from table
0 categories read from map exist in selection from table
0 categories read from map don't exist in selection from
table
0 records updated/inserted
0 update/insert errors ] 

*****************************

My (first!) script:

#!/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 "Script to list all unique pixel values in raster(s)
covered by vector areas of interest"

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 2

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

# Source for vector area(s) of interest
AREASVECT=$1

# Raster(s) from which pixel values will be extracted
RASTERSRC1=$2
RASTERSRC2=$3

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

echo "Creating new table to upload pixel values... "

v.db.addtable map=$AREASVECT layer=2 'columns=areacat
integer, pixval integer'
v.db.update map=$AREASVECT layer=2 column=areacat
qcolumn="cat"

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

# March region to AREASRAST
echo "Prepare region for analysis..."
g.remove MASK
g.region vect=$AREASVECT -p && echo "
                              ...matched to "$AREASVECT""


# 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:"
# where "i" the "cat" value derived from VECTOR1
for i in $CATi; do
	echo "\nProcessing Area [$i] out of [$CATMAX]"
	# Step 2(a). Masking areas of interest
	echo "*Masking..."
	r.mapcalc MASK="if(AREASRAST==$i, 1, null())"
	# 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
Pixel_Values_"$RASTERSRC1"_$i.csv)"
	r.stats -nxg input=$RASTERSRC1,$RASTERSRC2 fs="|"
output=Pixel_Values_$i.csv
	g.remove MASK && echo "otherwise r.mapcalc will work on
it!"

done

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

*******************************************8

Is there maybe the interest to integrate this to
v.rast.stats ?


On Sun, 20 Jan 2008 21:43:01 -0800 (PST)
 Hamish <hamish_b at yahoo.com> wrote:
> G. Allegri wrote:
> > About point 4), you can do it simply adding
> (v.db.addcol) at the
> > beginning of your job, then insert something like this
> into the loop:
> > 
> > echo "UPDATE (your_table) SET
> pixval=(what_you_want_to_upload) WHERE
> > FID=$i" | db.execute
> 
> 
> fyi there is v.db.update which may be easier to use than
> raw
> db.execute.
> 
> 
> 
> Hamish