Here is my situation,

I have 10 multi-spectral images. For each one I am calculating three different vegetation indices. Those 3 raster objects share resolution, projection, extent and dimension (they are produced with raster calculations).

I am using this code to list them (not stack since I will use lapply later)


The next step in my methodology is to upload a shp representing the field boundaries to derive the mean, min and max value inside each polygon FOR EACH raster of the list.

The code I am using for this is:


Statistics<-function (VI) {

    meancalc<-extract(VI, shp, fun=mean)
    mincalc<-extract(VI, shp, fun=min)
    maxcalc<-extract(VI, shp, fun=max)

    shp@data<-data.frame(shp@data, mean=meancalc)
    shp@data<-data.frame(shp@data, min=meancalc)
    shp@data<-data.frame(shp@data, max=meancalc)

   writeOGR(shp, getwd(), layer="mylayer", driver="ESRI Shapefile")

I run the function using output<-lapply(VI, FUN=Statistics)

The output is the shp with three new attributes (mean, min and max). However, it is obviously overwriting the calculations each time lapply takes a new raster as input for the function.

What I want is the shapefile with individual attributes for each raster input. Something like: Index1mean, Index1min, Index1max + Index2mean, Index2min, Index2max + Index3mean, Index3min, Index3max

I know a solution would be to link the output name of the shp to the input raster name so that three different shp are created for each iteration. However, I would prefer to keep everything in the same shp

Any idea how to implement this on the code?

1 Answer 1


You may better create the entire dataframe of data out of the shapefile. And then, add it to the shapefile. Something like:

VI <- as.list(Index1,Index2,Index3)
shp <- readOGR("myshapefile.shp")

# output of the function is only dataframe for each Index
Statistics <- function (VI) {
    mean = extract(VI, shp, fun=mean)
    min = extract(VI, shp, fun=min)
    max = extract(VI, shp, fun=max)
output <- lapply(VI, FUN=Statistics)
# Combine lists into dataframe with multiple columns
out.df <- do.call("cbind", output)
# Change names of the df
names(out.df) <- paste(
  rep(c("Index1","Index2","Index3"), each = 3),
  rep(c("mean", "min", "max"), 3),
  sep = ".")
# Add in the shapefile
shp@data <- data.frame(shp@data, out.df)
writeOGR(shp, getwd(), layer="mylayer", driver="ESRI Shapefile")
  • Thanks for your solution. It works perfectly. Just something that did not: when doing names(out.df)<-paste(....... the output is Index1.mean, Index1.mean, Index1.mean, Index2.min, Index2.min, index2.min, Index3.max, Index3.max, Index3.max. I think the rep is not working in the way we want. I have solved that by using names(out.df)<-c(...) containing the names in the order I want. I would prefer a solution like yours but so far it is a minor problem that is easy to solve. Thanks again for your smart solution
    – GCGM
    Oct 24, 2017 at 7:57
  • 1
    Normally, if you use each for the Index but times for the "mean, min, max" in the rep function, this should output in the order you want. Oct 25, 2017 at 10:11
  • I would point out that it is very inefficient to have 3 calls to extract. Extracting values from a raster is a bottleneck so, doing it multiple times can significantly slow things down. A more efficient way would be to extract the values to an object and then operate on the object using something like lapply. Something like: x=extract(VI, shp); data.frame( min=unlist(lapply(x, min)),max=unlist(lapply(x, max))) Sep 16, 2019 at 22:54

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.