1

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?

3

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")
| improve this answer | |
  • 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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