Save lappy results without overwritting - R

Question

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)

VI<-as.list(Index1,Index2,Index3)

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:

shp<-readOGR("myshapefile.shp")

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?

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) {
data.frame(
    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")