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 '20 at 8:00
  • Welcome to SE! Can you post your R code so far? – GISKid Jun 11 '20 at 23:40
  • @Spacedman: No, no colour images. Grey Satellite images.they are on the same grid and extent. I cropped them already. – catherine Jun 12 '20 at 9:21
  • @GISKid: Thanks. – catherine Jun 12 '20 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 '20 at 12:00

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.