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 <[EMAIL PROTECTED]> 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
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user
