Ariel,

Whether you do stack or addLayer, each file needs to be inspected by gdal
to get its georeference etc.
You can shortcut the process if you already know this, because all you
files have the same parameters. For example with a (do not try this at
home, all kinds of things can go wrong, but it might work in this case)
hack like this:

f <- list.files()
r <- raster(f[1])
s <- stack(r)
s@layers <- sapply(f, function(x) { r@file@name=f; r } )

Also, you may want to group your data in files that have multiple layers,
e.g. GTiff or ncdf. That would make it faster to create Raster* objects
from these files, and data access should also be much faster; although I am
not sure what the optimal number of layers would be (90000 may be too many
layers per file for optimal data access; but you might try ncdf for that).

ss = writeRaster(s, 'test.nc')

Robert



On Tue, Nov 1, 2011 at 4:46 PM, Ariel Ortiz-Bobea
<aortizbobea@arec.umd.edu>wrote:

> Hello everyone,
>
> I have about 90,000 individual grid files in GRIB format I'm importing and
> converting to raster format in R.
>
> Some of the calculations I need to do (clipping with polygons... thanks
> Robert J. Hijmans for valuable advice on this!) are best performed on many
> raster layers at once, e.g. over a RasterStack.
>
> I have tried "stack()" and "addLayer()" commands. "stack" seems to read all
> objects in the command line slowing the process down more than
> proportionately to object size (see code below) which I suppose explains
> how
> "addLayer" came about. However, the time it takes "addLayer" to add an
> additional layer is roughly proportional to the size of the RasterStack it
> is being added to... so when I have very "thick" stacks/bricks adding an
> extra layer takes some time.
>
> Would there be an alternative way to combine many multi-band raster layers
> more quickly (in R or outside of R) from many individual files and still
> manage to keep things orderly, i.e. being able to identify specific layers
> afterwards?
>
> Any thoughts would be greatly appreciated.
>
> Ariel
>
> #---------------
> # setup
>        library(raster)
>        r <- raster(ncol=139, nrow=97)
>        r[] <- 1:ncell(r)
>
>        # stack command
>        s1<-r
>        t1<- c()
>        while (nlayers(s1)<=200) {
>                print(nlayers(s1))
>                timer1 <- system.time( {
>                        s1 <- stack(s1,r)
>                        })
>                        t1 <- rbind(t1,timer1[3] )
>                }
>
>        # addLayer command
>        s2<-r
>        t2<- c()
>        while (nlayers(s2)<=200) {
>                print(nlayers(s2))
>                timer2 <- system.time( {
>                        s2 <- addLayer(s2,r)
>                        })
>                        t2 <- rbind(t2,timer2[3] )
>                }
>
>        # visualize
>        plot(1:200,t1, type="l", xlab="# of layers in object", ylab="time
> for
> adding extra layer (s)")
>        lines(1:200,t2,col="blue")
>
>        plot(1:200,t2, type="l", col="blue", xlab="# of layers in object",
> ylab="time for adding extra layer (s)")
>
>
>
> -----
> Ariel Ortiz-Bobea
> PhD Candidate in Agricultural & Resource Economics
> University of Maryland - College Park
> --
> View this message in context:
> http://r-sig-geo.2731867.n2.nabble.com/How-to-quickly-add-many-raster-layers-to-a-thick-RasterStack-tp6953662p6953662.html
> Sent from the R-sig-geo mailing list archive at Nabble.com.
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]

