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.


#load data

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$")


outfiles = file.path("C:/Users/cathe/Documents/Output",

for (i in seq_along(ALL_FILES)) { 
     r = crop(stack(ALL_FILES[i]), shape_data)
     writeRaster(r, filename=outfiles[i],
              datatype= "INT1U",

mean <- mean(r)


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

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

  • You say "Images" - are they colour images? What's the mean and SD of a set of colours? And are they all on the same grid and extent, or will they need warping to conform? Have you tried mean(r1,r2,r3,r4) ?
    – Spacedman
    Jun 11, 2020 at 8:00
  • Welcome to SE! Can you post your R code so far?
    – GISHuman
    Jun 11, 2020 at 23:40
  • @Spacedman: No, no colour images. Grey Satellite images.they are on the same grid and extent. I cropped them already. Jun 12, 2020 at 9:21
  • @GISKid: Thanks. Jun 12, 2020 at 9:22
  • How could you tell if it is correct? Make a small example that you can check by hand. Check that the 77th pixel in mean is the mean of the 77th pixel of the input rasters? Create some rasters with known values in a folder as a test set.
    – Spacedman
    Jun 12, 2020 at 12:00

1 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:

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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