# Zonal stats percentile values using R?

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?

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.

``````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)
``````

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")
``````
• Thank you! Can you explain what you're doing in the first section of this script when setting the p and r objects? Your script works, but I'm getting NULL results when reading in my own files. – happymappy Dec 30 '19 at 11:46
• @happymappy please see my clarification. This was just creating some example data but, functionally the data need to be a `SpatialPolygonsDataFrame` object, representing the zonal data, and a `RasterLayer` class object representing the raster value data. And, obviously they must share the same projection and align. – Jeffrey Evans Dec 30 '19 at 17:02