Calculation of mean and standard deviation of raster images in R

Question

I have a set of raster images (GeoTIFF, Landsat 1984-2018) which I cropped with my AOI using a shapefile.

Now I want to stack them and calculate the mean and standard deviation of every pixel of the stacked images (in space and time). How can I stack them and calculate the mean and standard deviation? For the outputfile I need one image for the mean and one for the standard deviation.

library(raster)
library(rgdal)

#load data

setwd("C:/Users/cathe/Documents/GEOTIFF_offsetcorrected")
dbase = "C:/Users/cathe/Documents/GEOTIFF_offsetcorrected"

#polygon with crop-extend

shape_data <- readOGR("C:/Users/cathe/OneDrive/Documents/ArcGIS/Projects/Cologne/WGS_1984_UTM_Zone_32N/Cologne.shp", stringsAsFactors=FALSE)

#load tif files

ALL_FILES <- list.files(path = dbase, pattern = "*\\.tif$|*.TIF$")

#output

outfiles = file.path("C:/Users/cathe/Documents/Output",
                     paste0(basename(tools::file_path_sans_ext(ALL_FILES)),
                            ".tif"))
#crop

for (i in seq_along(ALL_FILES)) { 
     r = crop(stack(ALL_FILES[i]), shape_data)
     writeRaster(r, filename=outfiles[i],
              bylayer=FALSE,
              format="GTiff",
              datatype= "INT1U",
              options="COMPRESS=ZIP",
              overwrite=TRUE)
}

mean <- mean(r)

plot(mean)

writeRaster(x = mean, filename = "mean1.tif", driver = "GeoTiff", overwrite=TRUE)

I get a result, but I am not sure if it is correct?

Answer 1

I think you are looking for calc function from library raster. If you want a mean raster from your raster stack you can put fun = mean, or if you want standard deviation you can put fun = sd, etc.

For example, your rasters are r1, r2, r3:

library(raster)
stacked <- stack(r1, r1, r3) # make a raster stack

# calculate a mean raster
meanR <- calc(stacked, fun = mean)

# make a new raster as standard deviation of stacked rasters
sdR <- calc(stacked, fun = sd)

# to calculate raster range
raR <- calc(stacked, fun = range) # this will  give you a list of two rasters, one for lower range and one for upper range. To get the final range raster run the following code.

rangeR <- abs(raR[[1]] - raR[[2]])

Now you can save the output using writeRaster function.