[R-sig-Geo] Calculating descriptive stats from many maps

[Robert 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:
> On Mon, 9 Feb 2009, Tomislav Hengl wrote:
>
>>
>> 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:
>>
>>> maps$Nsum <- rowSums(maps at data, na.rm=T, dims=1)
>>> maps$avg <- maps$Nsum/(length(names(meuse.grid at data))-1)
>>
>> You could also loop the sd, mean and quantile function over a range of
>> cells:
>>
>>> for(i in length(names(maps at data))) {
>>> maps at data$sd[i] <- sd(maps at data[i,])
>>> maps at data$mean[i] <- mean(maps at data[i,])
>>
>> ...
>>>
>>> }
>>
>> This could take a lot of time!
>
> 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
>
>>
>> 2) if your maps are rather large, try also using the SAGA function:
>>
>>> rsaga.get.usage(lib = "geostatistics_grid", module=5)
>>
>> 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
>>
>>> -----Original Message-----
>>> From: r-sig-geo-bounces at stat.math.ethz.ch
>>> [mailto:r-sig-geo-bounces at stat.math.ethz.ch] On Behalf
>>> Of Rainer M Krug
>>> Sent: Monday, February 09, 2009 4:58 PM
>>> To: R-sig-Geo at stat.math.ethz.ch
>>> Subject: [R-sig-Geo] Calculating descriptive stats from many maps
>>>
>>> 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
>>>
>>> --
>>> Rainer M. Krug, PhD (Conservation Ecology, SUN), MSc (Conservation
>>> Biology, UCT), Dipl. Phys. (Germany)
>>>
>>> Centre of Excellence for Invasion Biology
>>> Faculty of Science
>>> Natural Sciences Building
>>> Private Bag X1
>>> University of Stellenbosch
>>> Matieland 7602
>>> South Africa
>>>
>>> _______________________________________________
>>> R-sig-Geo mailing list
>>> R-sig-Geo at stat.math.ethz.ch
>>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at stat.math.ethz.ch
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
> --
> Roger Bivand
> Economic Geography Section, Department of Economics, Norwegian School of
> Economics and Business Administration, Helleveien 30, N-5045 Bergen,
> Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
> e-mail: Roger.Bivand at nhh.no
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at stat.math.ethz.ch
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>