[GRASS-user] Re: Calculating standard error of many maps?

[Rainer M Krug]

On Tue, Feb 10, 2009 at 12:06 PM, Glynn Clements
<glynn at gclements.plus.com> wrote:
>
> Moritz Lennert wrote:
>
>> >>> Just one clarification: I would like to calculate these descriptive
>> >>> stats for each cell, to obtain "variability maps".
>> >>> Rainer
>> >>> On Mon, Feb 9, 2009 at 5:57 PM, Rainer M Krug <r.m.krug at gmail.com
>> >>> <mailto:r.m.krug at gmail.com>> wrote:
>> >>>> Hi
>> >>>>
>> >>>> I have 25000 maps, generated by simulation predictions, and I would
>> >>>> like to calculate some descriptive stats, like mean, standard
>> >>>> deviation, median, quartiles.
>> >>>>
>> >>>> Is there an easy way of doing this in GRASS?
>> >>
>> >>
>> >> r.series ?
>> >>
>> >> You will probably have to use xargs to collate that many map names...
>> >
>> > Or you can use g.mlist if the names have a common prefix, suffix, etc.
>> > This assumes you want a single map that displays the result of all the
>> > input maps.
>> >
>> > E.g.
>> >
>> > r.series input="`g.mlist rast pattern='prefix_*' sep=,`"
>> > output=prefix_stddev method=stddev
>>
>> Yes, but if I'm not mistaken, this won't work with 25000 maps as it
>> would create a too long command line.
>
> 25000 maps will probably also exceed the resource limit for the number of
> open descriptors. It's common for user processes to be limited to 1024
> descriptors; you need to modify /etc/security/limits.conf to increase the
> limit.

I asked the same question on the R mailing list, and got the following
promising response:
###
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')

####

I am trying it at the moment with 1002 maps and, apart from the fact
that it takes some time (which I am sure it would in GRASS as well),
it seems to be working.

I will be looking at R.series again later, when I run the simulations
again. But I assume that there will be serious problems with 25000
maps - I might have to do it with a sample of maps.

Rainer


>
> --
> Glynn Clements <glynn at gclements.plus.com>
>



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