Zonal stats percentile values using R?

Question

I have a 1-m continuous DEM and a vector grid with ~160,000 250m cells. For each cell in the grid overlaying my DEM, I want to retrieve the 50th percentile pixel value. I'm curious if there's a way to retrieve this percentile value using zonal stats in R?

Answer 1

Here are two ways to do this. The raster extract function is the common approach. However, the zonal.stats function uses velox under the hood, which speeds things up considerably. You will not see much difference here but, with large problems it is a huge speed gain.

Add required libraries

library(raster)
library(sp)
library(spatialEco)

Here we create some example data. Please note that I am creating a raster that contains "zones" or polygons and then coercing it to a SpatialPolygonsDataFrame object. This is similar to what the raster::rasterToPolygon function does. Depending on your data, when operating on large rasters, an efficient approach is to have your zonal data represented as polygon vector data.

p <- raster(nrow=10, ncol=10)
  p[] <- runif(ncell(p)) * 10
    p <- as(p, "SpatialPolygonsDataFrame")  
      p <- p[p$layer > 9,]    
        r <- raster(nrow=100, ncol=100)
          r[] <- runif(ncell(r)) 
    plot(r)
      plot(p, add=TRUE, lwd=4) 

This defines a simple wrapper function for quantile that allows us to control the functions arguments while passing it to another function.

pct <- function(x, p=0.5, na.rm = TRUE) { quantile(x, p, na.rm = na.rm) }

This method extracts the raster data, and applies our pct function, for each zonal polygon using the extract function in the raster package.

raster::extract(r, p, fun=pct) 

This method extracts the raster data, and applies our pct function, for each zonal polygon using zonal.stats in the spatialEco package, using velox for the raster value extraction.

spatialEco::zonal.stats(x = p, y = r, stats = "pct")