Skip to content

Calculating descriptive stats from many maps

5 messages · Rainer M Krug, Tomislav Hengl, Roger Bivand +1 more

Original
Rainer M Krug
# 
Hi

I have 25000 maps, generated by simulation predictions, covering the
same area, and would like to calculate some descriptive stats, like
mean, standard deviation, median, quartiles of all cells, to create a
"variability map".

Is there an easy way of doing this in R?

Thanks,

Rainer
Tomislav Hengl
# 
Dear Rainer,

This is of course possible in R, and can be done in several ways:

1) for example, you can derive the average value using the rowSums function:
You could also loop the sd, mean and quantile function over a range of cells:
...
This could take a lot of time!

2) if your maps are rather large, try also using the SAGA function:
SAGA CMD 2.0.3
library path:   C:/Progra~1/saga_vc/modules
library name:   geostatistics_grid
module name :   Statistics for Grids

This is probably the fastest method you can use.

HTH

T. Hengl
Roger Bivand
#
On Mon, 9 Feb 2009, Tomislav Hengl wrote:

            

      
      
      
    

        
Tom, Rainer,

Yes, using sapply(slot(maps, "data"), summary) or similar, you get the 
map-wise statistics. But have I misunderstood, or are the statistics in 
question cell-wise? This would involve stacking subset areas for all 25' 
maps, wouldn't it? Brutally, a loop in readGDAL() from rgdal with offset= 
and region.dim= shifted? Is there a canned way to do this in the R-forge 
raster package (by the way, regularly one of the R-forge packages showing 
most developer activity)?

Roger

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

        

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Robert J. Hijmans
    

    

      
      #
    

  


  

      
Dear Rainer,

This is how can you can do it with the raster package

# install.packages("raster", repos="http://R-Forge.R-project.org")
require(raster)

# Try it for a few files first..
n <- 10

# create a list (or vector) of file names, e.g. :
fn <- list()
for (i in 1:n) { fn[i] <- paste('myfile', i, '.tif', sep='') }

# make a RasterStack
s <- stack(fn)

r1 <- mCalc(s, fun=mean)
r2 <- mCalc(s, fun=sd)

#r can be plotted, coerced to sp objects, etc.
plot(r1)

# or saved to file
r1 <- setFilename(r1, 'cellmeans.tif')
r1 <- writeRaster(r1, format='GTiff')


Robert

        
On Tue, Feb 10, 2009 at 1:05 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

            

      
      
      
    

        

      
      
    

  
  


  


      

    

    

      
      

        


  

        

      Rainer M Krug
    

    

      
      #
    

  


  

      
Thanks a lot to all of you.

You are right Roger, I need cell-wise statistics

I like the idea of the "raster" package, and I will try it out just now.

Concerning SAGA: I'll look into that if "raster" does not work (or is to slow).

I'll report back

Rainer

        
On Tue, Feb 10, 2009 at 4:07 AM, Robert Hijmans <r.hijmans at gmail.com> wrote: