I am using the following code to extract mean pixel values over an area represented by a polygon:

mean<-lapply(S2_stack, FUN=function (S2_stack) {data.frame(mean=extract(S2_stack, polygons2, fun=mean))})


> S2_stack
class      : RasterStack 
dimensions : 1454, 1595, 2319130, 4  (nrow, ncol, ncell, nlayers)
resolution : 10, 10  (x, y)
extent     : 744490, 760440, 4773400, 4787940  (xmin, xmax, ymin, ymax)
crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
names      : Subset_S2_MSIL2A_20180511T100029_MCI, Subset_S2_MSIL2A_20180710T100029_MCI, Subset_S2_MSIL2A_20180906T101021_MCI, Subset_S2_MSIL2A_20181115T101251_MCI 

and polygons2:

polygons2 <- readOGR('/path/Insitu_poly.shp')


> class(polygons2)
[1] "SpatialPolygonsDataFrame"
[1] "sp"

When I try to run my code I am getting the following error:

 Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ?extract? for signature ?"numeric", "SpatialPolygonsDataFrame"?

Not sure why. It seems extract needs a SpatialPolygonsDataFrame and as far as I know, it should run

Any idea?

  • @Aaron same error using mean<-lapply(S2_stack, FUN=function (S2_stack) {data.frame(mean=raster::extract(S2_stack, polygons2, fun=mean))}) Not sure how to proceed. Chekcint the documentation, I dont see a problem in the class of polygons2
    – GCGM
    Feb 26, 2020 at 15:20
  • 3
    Use simply extract(S2_stack, polygons2, fun=mean, df = T). Don't use lapply
    – aldo_tapia
    Feb 26, 2020 at 15:31
  • 1
    I think the problem is with the raster, not the polygons. I've never used lapply when extracting from a raster stack. Maybe try just doing mean <- extract(S2_stack, polygons2, fun=mean)
    – Nick
    Feb 26, 2020 at 15:40
  • 1
    Your stack has four bands. Do you want the average in each band separately? Or over all bands? lapply(1:4, FUN=function (i) {raster::extract(S2_stack[[i]], polygons2, fun=mean)}) maybe?
    – Spacedman
    Feb 26, 2020 at 16:32
  • Using the code without lapply worked. However, I do not get completley why as I am providing a list (of rasters) and a funciton to be applied. Please add the comment as an answer and I ll be glad to accept it.
    – GCGM
    Mar 2, 2020 at 10:47

2 Answers 2


Assuming that the correct way to proceed in this case is simply to do

mean <- extract(S2_stack, polygons2, fun=mean)

as said in the comments, however, if you want to pass your stack in lapply you have to unstack it first.

mean<-lapply(unstack(S2_stack), FUN=function (S2_stack) 
polygons2, fun=mean))})
  • 1
    You do not have to unstak the data first, simply define an iterator in lapply, lapply(1:nlayers(r), FUN=function (i) {raster::extract(r[[i]], p, fun=mean)}) But, using lapply is completely unnecessary and I would highly recommend using a different package entirely. The exact_extract function in the exactextractr package far outperforms raster::extract and benchmarks are quite notable for multiband data. The is also the idea of migrating to the terra package which will eventually replace raster anyway. Apr 16, 2021 at 19:09

I recommend you to use the exactextractr package instead of the extract function of the raster package which is extremely slow. Here is an example of how to use it:

if (!require("raster")) install.packages("raster")
if (!require("rgdal")) install.packages("rgdal")
if (!require("exactextractr")) install.packages("exactextractr")


I assume two variables loaded. img is a RasterStack with several layers and shp is a ESRI Shapefile containing polygons. This shapefile contains a field called class with the land cover. You can load this variables using: img <- raster::stack('file.tif') and shp <- rgdal::readOGR('file.shp').

# extract pixel values on a list of data.frames where
# length(ls.df) == number of polygons of the ESRI Shapefile
ls.df <- exact_extract(img, shp, include_cols="class", progress=TRUE)
# filter only full pixel (coverage_fraction==1)
for(i in 1:length(ls.df)){
ls.df[[i]] <- ls.df[[i]][ls.df[[i]]$coverage_fraction == 1,]
# create a single data.frame containing all the pixel values
df <- data.frame()
for(i in 1:length(ls.df)){
df <- rbind(df, ls.df[[i]])

The exact_extract function computes the coverave fraction of each pixel covered by the polygon. You remove this information using as follow:

# remove coverage_fraction column
df$coverage_fraction <- NULL

You can learn more about this package here: https://github.com/isciences/exactextractr

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.