[STATSGRASS] raster cross-stats

[Markus Neteler]

On Thu, Dec 13, 2001 at 02:09:23PM +0100, Agustin Lobo wrote:
> 
> On Thu, 13 Dec 2001, Roger Bivand wrote:
> 
> > Hi Agus!
> > 
> > Great to see you "here" too! 

Hi Agus

as Roger says: nice to see you again!
 
> Well, I actually wrote some (probably primitive) scripts
> linking Splus and Grass ten years ago! I did pass the scripts
> to Markus, and maybe they were at least useful to initially fuel the
> current (and better) link with R.

Doh! They must be somewhere. I am sorry, but that time (maybe two
years ago?) I wasn't aware that potential and my current interests.
Unfortunately I can't find the, do you have a copy?
 
> > I'm enclosing below bits of correspondence
> > between Markus, Helena and me of about two weeks ago, together with a
> > draft shell script to use R functions for GRASS raster layers. My current
> > feeling is that with R being easy to install and maintain, shell scripts
> > using R (here for the leics dataset):
> 
> > > library(GRASS)
> > > G <- gmeta()
> > > work <- rast.get(G, c("landcov", "topo"), c(T,F))
> > > res <- by(work$topo, work$landcov, summary)
> 
> etc.
> 
> Yes, I was aware of this option, which is very similar to 
> my current solution: I import the necessary rasters to R
> (using more general tools for raw binary images, as I 
> use a number of different GIS and RS packages) and there
> I compute the statistics.
> 
> But there are 2 problems with this approach:
> 
> 1. We duplicate the information, which is really bad in the
> case of images, as they can be huge (hyperspectral, multitemporal,
> high resolution ...).
> 
> 2. We import large objects to R, which is not precisely good and fast
> at dealing with large objects.
> 
> In fact, an ideal link would let R work on the Grass files, but this
> is probably asking for too much.

Roger, did you ever look at the "libgrass" library from
Frank Warmerdam? We recently updated it to connect
MapServer directly to GRASS (see it working here:
http://grass.itc.it/start.html
a primitive version only to proof that it is running).

Using this library 
http://grass.itc.it/related_projects.html
-> libgrass: Standalone GRASS 5 I/O Library to directly read/write GRASS
Databases. 

you may exchange data directly. See 
gdal/frmts/grass
GDAL how to use the library if you are interested.

> Therefore, a good solution would be performing this type of
> "low level" stats directely at the GIS level, get a table
> (of moderate size, my tables would have ~ 20000 rows) and
> port it to R. If we have, for each category, mean, size, var,
> median, mad ..., we can use R for further analysis. 
> 
> I do something similar with contingency tables: I use grass
> r.stat to creat a 3 col table (categ.1_raster_1, categ.1_raster_2
> nb_of_cells) and then a simple R function to compute
> a contingency table. This is much faster and simpler than importing
> the 2 rasters to R.
> 
> Opinions?
> 
> Agus

Markus