[R-sig-Geo] How to efficiently clip large grids/raster with polygons?

[Ariel Ortiz-Bobea]

Hello Robert,
Sorry for the long delay. This is so useful!
Thanks a lot for the example, this should help many others in the future.
Thanks again!
Ariel

---
Ariel Ortiz-Bobea
PhD Candidate, Agricultural & Resource Economics
University of Maryland - College Park

________________________________________
From: Robert Hijmans [via R-sig-geo] [ml-node+s2731867n6807080h6 at n2.nabble.com]
Sent: Sunday, September 18, 2011 11:13 PM
To: Ariel Ortiz Bobea
Subject: Re: How to efficiently clip large grids/raster with polygons?

Ariel, if you have 23 stacks, you could, in principle, combine these into a
single stack and then extract the values. But because your stack are already
very large (90,000 layers each), you might be better off considering an
alternative option: to determine the weights once, and then re-use them
while looping over the 23 stacks, as illustrated (without the loop) below.


library(raster)
r <- raster(ncol=36, nrow=18)
r[] <- 1:ncell(r)
cds1 <- rbind(c(-180,-20), c(-160,5), c(-60, 0), c(-160,-60), c(-180,-20))
cds2 <- rbind(c(80,0), c(100,60), c(120,0), c(120,-55), c(80,0))
polys <- SpatialPolygons(list(Polygons(list(Polygon(cds1)), 1),
Polygons(list(Polygon(cds2)), 2)))

# extract the cellnumbers and weights
v <- extract(r, polys, weights=TRUE, cellnumbers=TRUE)
head(v[[1]])

# example stack
s <- stack(r,r*2,r*3,r*4,r*5)

# normal approach
system.time( a <- extract(s, polys, weights=TRUE, fun=mean) )

# reuse of weights approach (assuming there are no NA values)
fun <- function(x) apply((s[v[[x]][, 'cell']] * v[[x]][,'weight']) / sum(
v[[x]][,'weight']), 2, sum, na.rm=TRUE)

# as above, but here taking care of possible NA values (the same cells for
all layers)
fun2 <- function(x) {
d <- na.omit(cbind(s[v[[x]][, 'cell']], v[[x]][,'weight']))
 w <- d[, ncol(d)]
d <- d[, -ncol(d)]
apply(d*w/sum(w), 2, sum, na.rm=TRUE)
}

# much faster than the example above
system.time( b <- sapply(1:length(v), fun ))
system.time( d <- sapply(1:length(v), fun2 ))


 a
b
d

# results are the same (for this example!)
t(b) == a



Probably a good idea to test this first for a few counties.

Best, Robert

On Sun, Sep 18, 2011 at 11:35 AM, Ariel Ortiz-Bobea <
[hidden email]</user/SendEmail.jtp?type=node&node=6807080&i=0>> wrote:

> Thanks Robert,
>
> It works great. It took the same time to clip the RasterStack than the
> single RasterLayer. I read your "Introduction to the ’raster’ package
> (version 1.9-11)" and I found very useful too to understand the different
> objects. Thanks for developing this resource!
>
> I wonder about a way to combine several variables (23) and time periods
> (~90,000) to efficiently do the polygon clipping on a large
> RasterBrick/Stacks. I feel that I need to choose between:
>
> 1- performing the clipping on 23 RasterStacks (one for each variable, e.g.
> temperature) in which each band corresponds to a time period (90,000 of
> them)
>
> 2- performing the clipping on 90,000 RasterStacks (one for each time
> period,
> e.g. 1980/1/1 09:00) in which each band corresponds to a variable (23 of
> them).
>
> ... and option 1 would be expectedly faster. But, is there a way to combine
> both options to do the clipping "ALL at once" of all the time periods and
> variables? I read about the "setZ" function for time-tagging layers and the
> ability to perform functions for each raster cell over time... but I'm not
> sure I can make use of it in this context (stacking 2 time-period layers
> doubles my band number... basically, I feel it does not lead to an
> additional dimension in the structure of the object).
>
> Any thoughts would be greatly appreciated!
>
> Regards,
>
> Ariel
>
>
>
> -----
> 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-efficiently-clip-large-grids-raster-with-polygons-tp6804923p6806329.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=6807080&i=1>
> 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=6807080&i=2>
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-efficiently-clip-large-grids-raster-with-polygons-tp6804923p6807080.html
To unsubscribe from How to efficiently clip large grids/raster with polygons?, click here<http://r-sig-geo.2731867.n2.nabble.com/template/NamlServlet.jtp?macro=unsubscribe_by_code&node=6804923&code=YW9ydGl6Ym9iZWFAYXJlYy51bWQuZWR1fDY4MDQ5MjN8LTYxMzExNzA2NQ==>.


-----
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-efficiently-clip-large-grids-raster-with-polygons-tp6804923p6926442.html
Sent from the R-sig-geo mailing list archive at Nabble.com.