[R-sig-Geo] How to quickly add many raster layers to a thick RasterStack?

[Ariel Ortiz-Bobea]

Thanks Robert,

This works fairly well and the direct access is much faster. I just need to manually specify the extent and projection.

Each file has several bands and I should specify which band I'll be importing (otherwise I only import the first band) this way:

r <- raster(f[1], band=2) # for example

Thanks!

AOB

---
Ariel Ortiz-Bobea
PhD Candidate, Agricultural & Resource Economics
University of Maryland - College Park
________________________________________
From: Robert Hijmans [via R-sig-geo] [ml-node+s2731867n6971221h77 at n2.nabble.com]
Sent: Monday, November 07, 2011 12:14 PM
To: Ariel Ortiz Bobea
Subject: Re: How to quickly add many raster layers to a thick RasterStack?

correction, it should be:

s at layers <- sapply(f, function(x) { r at file@name=x; r } )

On Sun, Nov 6, 2011 at 11:04 PM, Robert J. Hijmans <[hidden email]</user/SendEmail.jtp?type=node&node=6971221&i=0>>wrote:

> 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 at layers <- sapply(f, function(x) { r at 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 <
> [hidden email]</user/SendEmail.jtp?type=node&node=6971221&i=1>> 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
>> [hidden email]</user/SendEmail.jtp?type=node&node=6971221&i=2>
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]</user/SendEmail.jtp?type=node&node=6971221&i=3>
https://stat.ethz.ch/mailman/listinfo/r-sig-geo


________________________________
If you reply to this email, your message will be added to the discussion below:
http://r-sig-geo.2731867.n2.nabble.com/How-to-quickly-add-many-raster-layers-to-a-thick-RasterStack-tp6953662p6971221.html
To unsubscribe from How to quickly add many raster layers to a thick RasterStack?, click here<http://r-sig-geo.2731867.n2.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=6953662&code=YW9ydGl6Ym9iZWFAYXJlYy51bWQuZWR1fDY5NTM2NjJ8LTYxMzExNzA2NQ==>.


-----
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-tp6953662p6982713.html
Sent from the R-sig-geo mailing list archive at Nabble.com.