I'm working with RasterStack of 300 simulations sims_trend. and I also have a SpatialPolygonsDataFrame limits. I'm trying to calculate quantiles for each polygon through the raster stack.

I know that raster::extract() is close to what I want. But, for what I know, it only does the calculations for the pixels inside the polygon for each layer. As an example I'll show the output of raster::extract() applied to my sims_trend.

b <- extract(sims_trend, limits, method="simple",
             fun=function(x, na.rm){quantile(x, probs=.975, na.rm=T)}, df=T)
> b[1,1:4]

 layer.1  layer.2  layer.3  layer.4 
21.00252 22.37075 26.89889 24.42825 

So the output is the quantile for each polygon of each layer, but what I wanted is the quantile for each polygon for all layers (basically one column)

1 Answer 1


Don't pass an argument to fun and the object will then be a list, with an element for each polygon, containing the pixel values. You can then pass a function tolapply that summarizes the values however you want. For a multiband (stack/brick) object, the simplest approach would be to collapse the data.frame into a vector and operate on that as representing "all" values intersecting a given polygon thus, removing the by layer structure of the results.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.