I have a raster image (Sentinel-2 band 4) and a shapefile covering the same area than the image (same projection and extension). The shapefile has different polygons which are agricultural fields.

I want to create a new file where, for each polygon, the mean value of the pixels that are inside it is calculated.

Could you please help me to find an appropriate R code for that?

  • 2
    See ?raster::extract for something like 'extract(ras, poly, fun = "mean")'
    – mdsumner
    Commented Apr 18, 2017 at 13:17

3 Answers 3


The raster::extract function, when applied to polygons or with the buffer argument, returns a list object where each element in the list contains a vector of the raster values intersecting the polygon. If the input raster object was a stack or brick, containing multiple rasters, the list elements are a matrix rather than a vector.

Providing the fun argument to extract is just a way of aggregating or summarizing the values within the function, without having to the deal with the list object. This really only works with simple functions that operate on vectors, not matrices.

Here we can work through what extract returns and manipulate the results. First add the required packages and create some example data. We will create a raster (r) and some polygons (poly).

poly <- raster(nrow=10, ncol=10)
  poly[] <- runif(ncell(poly)) * 10
    poly <- rasterToPolygons(poly, fun=function(x){x > 9})
      r <- raster(nrow=100, ncol=100)
        r[] <- runif(ncell(r)) 
  plot(poly, add=TRUE, lwd=4) 

Now we can look at what extract returns when the fun argument is not provided.

( v <- extract(r, poly) )

You can see that it is a list with vectors of raster values corresponding to each polygon. Please note that the list elements are ordered, that is to say that the firs list element corresponds to the first polygon. Because of this, any summary of the list will stay ordered with the polygon data.

For illustration purposes, let's write our own function that calculates the proportion of values above a given threshold.

pct <- function(x, p=0.30) {
  if ( length(x[x >= p]) < 1 )  return(0) 
    if ( length(x[x >= p]) == length(x) ) return(1) 
     else return( length(x[x >= p]) / length(x) ) 

Now, we can apply this function to the list object using the lapply function. Since the object is a list it will return a list so, we simply wrap the call in unlist to coerce it into a vector.

unlist(lapply(v, pct))

Since this data is ordered with our polygons we can just assign it into a new column "pcts".

poly@data$pcts <- unlist(lapply(v, pct))
  spplot(poly, "pcts")

Here is a quick example of what happens when the list contains matrices resulting from passing extract a multiband object. In this case, any function passed to lapply would have to account for the different data structure eg., calling a specific column or operating on the entire matrix.

r.stack <- stack(r,r,r)
v <- extract(r.stack, poly) 
lapply(v, head) #display first 6 lines of each matrix

This is relevant because something like Sentinel data is multiband. To return the mean for each band and add it to the polygons data.frame you could apply something like the colMeans function and, using do.call, coerce to a matrix/data.frame.

as.data.frame(do.call(rbind, lapply(v, colMeans))) 
( poly@data <- data.frame(poly@data, do.call(rbind, lapply(v, colMeans))) )

You can also allow the extract function to to the recycling of the raster layer using the df=TRUE argument to achieve the same results.

extract(r.stack, poly, fun=mean, df=TRUE) 
  • Thanks for your detailed answer. Just came back to it to review some concepts and it helped me again. Just a small comment, does the approach colMeans creates the same result as using directly fun=mean in the extract function? In case yes, wouldn't it be more straightforward?
    – GCGM
    Commented Dec 19, 2017 at 9:15
  • 2
    @GCGM on a per raster basis yes but you cannot apply the mean, using the fun argument to a stack or brick. Commented Dec 19, 2017 at 11:59
  • Thanks, just to point out, in that type of situation (raster list) I am using lapply - e.g. - result<-lapply(rasterlist, FUN=function (rasterlist) {extract(rasterlist, polygons, fun=mean, df=TRUE)}
    – GCGM
    Commented Dec 19, 2017 at 12:08

To elaborate on the comment, here's a short sample code to answer the question:

#import raster file, we'll call it sentinel_2
#import polygon file, we'll call it agri_fields
library(raster) # use the raster library
agri_s2=extract(x=sentinel_2, y=agri_fields, fun=mean, df=TRUE)#

The issue I have with the output is that we lose any information about the original polygon file in this table. To solve that, you can add a new column to the extracted table with a polygon ID that you might have in your polygon file.

  • 1
    You do not "loose" any information about the original polygon data. Based on the lapply function used in the raster::extract to apply the fun argument, the data stays ordered so, the relationship between the output and the polygon data is one to one. Besides, with the df=TRUE option the ID column represents the row.names of the polygon objects @data slot data.frame, which provides a link ID that can be used in merge. Commented Dec 6, 2017 at 23:15
  • I don't understand the syntax @data slot data.frame. I also don't see rownames in my outputs. In my sample polygon dataset, I have rows 0 to 19. If I call rownames(polygonfile), I get NULL. In my sample extraction file with df=TRUE, I have rows 1 to 20. My familiariaty with R is not advanced enough to follow your sample code in your answer. Commented Dec 8, 2017 at 16:02
  • All sp class **DataFrame objects have a data slot that contains the attribute data of the object. The rownames of the dataframe correspond with the ID's in each feature class object. If you look at rownames(polygonfile@data) it returns the rownames. Many functions that track a unique ID will use the rownames as it is an element that can be relied on to exists and be unique. There has been an effort to create standard methods to make it so you do not need to call the data.frame through calling the slot but, it is good to be aware of it. Please look at my answer to see how to use a slot. Commented Dec 8, 2017 at 16:22

Following the answer suggested by @HMSP, great solution is now offered by terra package, and using terra native formats for rasters and vectors. First convert raster and vector to terra format, and then run extraction by terra::extract(raster, polygon). Note, that the output here is not a raster (unlike to QGIS or ArcGIS extraction) but a dataframe. If you set list = FALSE, all raster values will be added to a dataframe, with unique ID indicating individual features within the polygon shapefile.


# convert to native format
poly<- vect(grids)
rs  <- rast(CVsteepslope)

# extract

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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