How to Find the Average Raster Value of an area defined by a shapefile using R?

Question

I have a set of raster images which represent a specific month over the years, and I want to make a timeline of the average values of an area using a shapefile.

How do I extract the values from the rasters and import them in R in a manner that I can use them?

Answer 1

Here is example code. It is fairly straight forward to adapt this code to work in a loop for processing all of your rasters. If your rasters share a common extent and resolution you can create a raster stack and loop through the bands in the stack. To create a vector containing all rasters in a directory, in a specific format, you can use "list.files" and then pass this vector to stack.

Example:

rlist=list.files(getwd(), pattern="img$", full.names=TRUE) 

r <- stack(rlist)   


    # Add required libraries
    require(raster)
    require(sp)
    require(rgdal)

    # Set working directory, raster, in and out shapefiles
    setwd("C:/test")
    inshp="MyPolys"
    outshp="PolyMeans"
    rdata <- "Year2012.img"

    # Read polygon feature class shapefile
    sdata <- readOGR(dsn=getwd(), layer=inshp)

    # Read raster
    r <- raster(rdata)

    # Extract raster values to list object
    r.vals <- extract(r, sdata)

    # Use list apply to calculate mean for each polygon
    r.mean <- lapply(r.vals, FUN=mean)

    # Join mean values to polygon data
    sdata@data <- data.frame(sdata@data, m2012=r.mean)

    # Write results
    writeOGR(sdata, getwd(), outshp, driver="ESRI Shapefile", check_exists=TRUE, 
         overwrite_layer=TRUE)

Answer 2

Read the shapefile into a SpatialPolygonsDataFrame (readOGR function from package rgdal)

Read the raster into a Raster object (raster function from package raster)

Use extract(raster, spdf) to get the grid cells under each polygon. Then run mean on them.

Repeat over your set of raster images...

Answer 3

Both answer above are great, I wanted to point out that the raster::extract() function does allow you to call functions within it.

r1 <- raster
p1 <- polygon feature

p1$mean_val <- raster::extract(x = r1, y = p1, FUN=mean) 

Makes for a nice one liner for adding your summary values back to your polygon feature.