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

GRASS GRASS user
1

Dear list,

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.

Any suggestions?

Regards

Chris

Express yourself instantly with MSN Messenger! MSN Messenger

2

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

3

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@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

4

On Tue, 2008-03-04 at 17:08 -0800, Dylan Beaudette wrote:

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.

Although Jose Luis Gomez replied to a question of mine about this in the
past, I never really gave it a try. Now I realise it's fast and does
what I want :wink:

I am still satisfied doing my attempt to script because I learnt 2-3
things more.

Cheers,

Dylan

5

Hey guys, wanted to echo this -- we are continuing development on starspan again -- it works with gdal 1.5, among other improvements, and if you have any feature requests (or if anyone can help us make this into a grass package :), please email me or Carlos Rueda!

--j

Nikos Alexandris wrote:

On Tue, 2008-03-04 at 17:08 -0800, Dylan Beaudette wrote:

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.

Although Jose Luis Gomez replied to a question of mine about this in the
past, I never really gave it a try. Now I realise it's fast and does
what I want :wink:

I am still satisfied doing my attempt to script because I learnt 2-3
things more.

Cheers,

Dylan

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--

Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn307@hotmail.com, Gchat: jgrn307

6

Jonathan,

I've seen mention of starspan previously and now I think it's time for me to learn more. Let me pose a problem to you to see if starspan would be of help. I am working on a modeling project with a couple of other people that involves the following:

(1) downloading and decoding gridded fields of numerical weather prediction (NWP) model output
(2) ingesting the decoded data into GRASS to calculate basin average precipitation & temperature (separate grids) for each subbasin
(3) writing out all calculated basin average values for each grid to separate files (one file contains one time step of data for all subbasins) for both temperature & precipitation
(4) for data management reasons, the files from (3) are written to a PostgreSQL database
(5) once all time steps of gridded precipitation and temperature field data are written to the PostgreSQL database, another process collects the data and generates individual ascii time series files for each subbasin for both temperature & precipitation
(6) once (5) is completed a hydrologic model runs using the temperature & precipitation time series as input and hydrologic forecasts are generated.

This is suppose to be a *real-time* process. The problem I am having is a matter of scale. What I did not say is that there are 12 different sets of NWP output covering a period of 168 hours at 6-hour time steps for both temperature & precipitation. So, this means I must process 12*2*(168/6) = 672 grids. Also, I need the mean areal values for 686 subbasins within the domain for each of the grids.

Steps (2-4) take about 20 seconds total for each of the grids… which is ~7.5 hours total
Step-5 also takes about 20 seconds for each time series file… which is ~7.5 hours total

So, about 15 hours total. Now, I can cut this time in half by running the processing of the temperature & precipitation grids and generating their separate time series files in parallel, rather than sequentially. So, I can get to about 7 hours fairly easily — what I am shooting for is to get the processing time from 7 hours to about 3 hours (or less).

I need to more efficiently generate the many basin average time series files from the numerous grids. Can starspan help by reducing the time to calculate the the basin average values faster?

I would also appreciate any/all suggestions on how to efficiently go from 'starspan generated basin average values' to my time series files. Realize, of course, the the individual grids are only a slice in time, so I have to track the grids and their resulting individual basin values (in time) to generate the time series files.

To compound the problem, very soon, I need to add model grids from an additional 21 models, bringing the total from 12 to 33!

Regards,
Tom

Jonathan Greenberg wrote:

Hey guys, wanted to echo this -- we are continuing development on starspan again -- it works with gdal 1.5, among other improvements, and if you have any feature requests (or if anyone can help us make this into a grass package :), please email me or Carlos Rueda!

--j

Nikos Alexandris wrote:

On Tue, 2008-03-04 at 17:08 -0800, Dylan Beaudette wrote:

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.

Although Jose Luis Gomez replied to a question of mine about this in the
past, I never really gave it a try. Now I realise it's fast and does
what I want :wink:

I am still satisfied doing my attempt to script because I learnt 2-3
things more.

Cheers,

Dylan

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--
Thomas E Adams
National Weather Service
Ohio River Forecast Center
1901 South State Route 134
Wilmington, OH 45177

EMAIL: thomas.adams@noaa.gov

VOICE: 937-383-0528
FAX: 937-383-0033

7

Hello Jonathan,

I just tried to see if I could grab the latest starspan via cvs. The cvs commands do not seem to be right. It says:

  cvs -d :pserver:anonymous@cvs1:/cvsroot/starspan login
  cvs -d :pserver:anonymous@cvs1:/cvsroot/starspan checkout modulename

Shouldn't cvs1 be a fully qualified domain name?

Bob

Bob Moskovitz
Seismic Hazard Evaluation Project
California Geological Survey
http://gmw.consrv.ca.gov/shmp

CONFIDENTIALITY NOTICE: This communication is intended only for the use of the individual or entity to which it is addressed. This message contains information from the State of California, California Geological Survey, which may be privileged, confidential and exempt from disclosure under applicable law, including the Electronic Communications Privacy Act. If the reader of this communication is not the intended recipient, you are hereby notified that any dissemination, distribution, or copying of this communication is strictly prohibited.

-----Original Message-----
From: grass-user-bounces@lists.osgeo.org
[mailto:grass-user-bounces@lists.osgeo.org]On Behalf Of Jonathan
Greenberg
Sent: Tuesday, March 04, 2008 8:18 PM
To: Nikos Alexandris
Cc: grass-user@lists.osgeo.org
Subject: Re: [GRASS-user] please help: v.rast.stats

Hey guys, wanted to echo this -- we are continuing development on
starspan again -- it works with gdal 1.5, among other
improvements, and
if you have any feature requests (or if anyone can help us make this
into a grass package :), please email me or Carlos Rueda!

--j

Nikos Alexandris wrote:
> On Tue, 2008-03-04 at 17:08 -0800, Dylan Beaudette wrote:
>> 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.
>
> Although Jose Luis Gomez replied to a question of mine
about this in the
> past, I never really gave it a try. Now I realise it's fast and does
> what I want :wink:
>
> I am still satisfied doing my attempt to script because I learnt 2-3
> things more.
>
>> Cheers,
>>
>> Dylan
>
>
> _______________________________________________
> grass-user mailing list
> grass-user@lists.osgeo.org
> http://lists.osgeo.org/mailman/listinfo/grass-user
>

--

Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn307@hotmail.com, Gchat: jgrn307
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

8

Honestly, I have no idea if our cvs works -- just grab it from our website: http://starspan.casil.ucdavis.edu/

--j

Moskovitz, Bob wrote:

Hello Jonathan,

I just tried to see if I could grab the latest starspan via cvs. The cvs commands do not seem to be right. It says:

  cvs -d :pserver:anonymous@cvs1:/cvsroot/starspan login
  cvs -d :pserver:anonymous@cvs1:/cvsroot/starspan checkout modulename

Shouldn't cvs1 be a fully qualified domain name?

Bob

Bob Moskovitz
Seismic Hazard Evaluation Project
California Geological Survey
http://gmw.consrv.ca.gov/shmp

CONFIDENTIALITY NOTICE: This communication is intended only for the use of the individual or entity to which it is addressed. This message contains information from the State of California, California Geological Survey, which may be privileged, confidential and exempt from disclosure under applicable law, including the Electronic Communications Privacy Act. If the reader of this communication is not the intended recipient, you are hereby notified that any dissemination, distribution, or copying of this communication is strictly prohibited.

-----Original Message-----
From: grass-user-bounces@lists.osgeo.org
[mailto:grass-user-bounces@lists.osgeo.org]On Behalf Of Jonathan
Greenberg
Sent: Tuesday, March 04, 2008 8:18 PM
To: Nikos Alexandris
Cc: grass-user@lists.osgeo.org
Subject: Re: [GRASS-user] please help: v.rast.stats

Hey guys, wanted to echo this -- we are continuing development on starspan again -- it works with gdal 1.5, among other improvements, and if you have any feature requests (or if anyone can help us make this into a grass package :), please email me or Carlos Rueda!

--j

Nikos Alexandris wrote:

On Tue, 2008-03-04 at 17:08 -0800, Dylan Beaudette wrote:

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.

Although Jose Luis Gomez replied to a question of mine

about this in the

past, I never really gave it a try. Now I realise it's fast and does
what I want :wink:

I am still satisfied doing my attempt to script because I learnt 2-3
things more.

Cheers,

Dylan

_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--

Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn307@hotmail.com, Gchat: jgrn307
_______________________________________________
grass-user mailing list
grass-user@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/grass-user

--

Jonathan A. Greenberg, PhD
Postdoctoral Scholar
Center for Spatial Technologies and Remote Sensing (CSTARS)
University of California, Davis
One Shields Avenue
The Barn, Room 250N
Davis, CA 95616
Cell: 415-794-5043
AIM: jgrn307, MSN: jgrn307@hotmail.com, Gchat: jgrn307

9

Hey guys, wanted to echo this -- we are continuing development on
starspan again -- it works with gdal 1.5, among other improvements, and
if you have any feature requests (or if anyone can help us make this
into a grass package :), please email me or Carlos Rueda!

--j

Dear Jonathan,

how about calculating basic stats (area size per class --given that a
"class" column exists--) on a single polygon vector?

Nikos
--
View this message in context: http://www.nabble.com/please-help%3A-v.rast.stats-tp15826629p15950816.html
Sent from the Grass - Users mailing list archive at Nabble.com.