Quantile of a polygon over a raster stack in R

Question

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)

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.