[R-sig-Geo] Re: [GRASSLIST:5390] Script to 'sample' several grids with resultsaved to text output?

[Roger Bivand]

On Mon, 17 Jan 2005, Sander Oom wrote:

> Dear R-geo and GRASS users,
> 
> -Before trying to reinvent the wheel I would like to check whether 
> someone has tried this before.-
> 
> I have several rasters with environmental data. I also have point data, 
> representing animal locations (collected with GPS collars).
> 
> I would like to use GRASS to 'sample' (apologies for the ArcInfo AML 
> terminology, but it might help people to get what I need) the rasters 
> with the point data and save the results in a text file for further 
> analysis with R.
> 
> The resulting text file should contain a field/attribute of my choice 
> from the points data file (some sort of ID for each point) and a 
> field/attribute for each of the sampled rasters.
> 
> Preferably I would use an R script that calls GRASS, as both point data 
> collection and further analysis is done within R. However a GRASS script 
> working on input and output text files might also be sufficient.

Something like the following works with GRASS 5.4 and R 2.0.1, using the 
spearfish test data set for GRASS 5.3:

> rtmp <- tempfile()
> system(paste("s.out.ascii bugsites | r.what -f vegcover,geology >", 
+ rtmp))
> system(paste("tail -5", rtmp))
601228|4917635||5|mixed forest|5|limestone
603146|4917453||3|coniferous forest|5|limestone
605627|4916428||3|coniferous forest|5|limestone
607098|4917437||3|coniferous forest|4|sandstone
608471|4915733||3|coniferous forest|8|claysand
> input <- scan(rtmp, sep="|", what=list(x=numeric(0), y=numeric(0), 
+ NULL, NULL, vegcover=character(0), NULL, geology=character(0)))
Read 90 records 
> str(input)
List of 7
 $ x       : num [1:90] 590232 590430 590529 590546 590612 ...
 $ y       : num [1:90] 4915039 4915204 4914625 4915353 4915320 ...
 $         : NULL
 $         : NULL
 $ vegcover: chr [1:90] "coniferous forest" "coniferous forest" 
"coniferous forest" "coniferous forest" ...
 $         : NULL
 $ geology : chr [1:90] "igneous" "igneous" "igneous" "sandstone" ...
> input1 <- data.frame(input[c(1,2,5,7)])
> str(input1)
`data.frame':   90 obs. of  4 variables:
 $ x       : num  590232 590430 590529 590546 590612 ...
 $ y       : num  4915039 4915204 4914625 4915353 4915320 ...
 $ vegcover: Factor w/ 4 levels "coniferous forest",..: 1 1 1 1 1 4 1 1 1 1
 $ geology : Factor w/ 5 levels "claysand","igneous",..: 2 2 2 5 5 5 5 3 2 5

If I had figured out how to read the output of r.what into s.in.ascii, 
storing the data as a sites file, then using the R GRASS package and 
sites.get() would have been neater, but this does work.

Hope this helps,

Roger Bivand

> 
> Any help would be much appreciated,
> 
> Sander Oom.
> 
> PS: I am a total novice with GRASS.
> 
> 

-- 
Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Breiviksveien 40, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 93 93
e-mail: Roger.Bivand at nhh.no