[R-sig-Geo] Iterating through raster stack to map values +/- 1/2/3 standard deviations from mean

[Philip Monk]

I would like to create two stacks of raster images, each containing 12
images in total (one for each month over a one-year time period), then
process each image to:

*) display pixels which are 1/2/3 standard deviations from the mean (I
would amend and run the code several times to plot different values -
not all in one go)
*) add a mask to indicate the aoi
*) crop the image to the aoi
*) save each raster as an individual file.

The images are geoTIFFs created in and exported from Google Earth
Engine (GEE).  Each pixel is a calculated value of Land Surface
Temperature derived from Landsat 8 (Collection 2) scenes.  They are
identical in extent et al.

The mask is a shapefile exported from GEE.  it has the same extent but
a different crs.  The code below addresses that, but I have also used
QGIS to dissolve its geometry to create an outline (for the purposes
of indicating the aoi on the mapped data as its structure is not
relevant for this analysis).

At present I can process an individual image:

####

library(raster)
library(rgdal)

LST   <- raster("temperature.tif")
Mask  <- readOGR("mask.shp")

# For mask not modified in QGIS
crs(LST)
crs(Mask)
lst = projectRaster(LST, crs=crs(Mask))

plot(lst)
plot(Mask, add=TRUE)

masked <- mask(x = lst, mask = Mask, inverse=FALSE)
plot(masked)

cropped <- crop(masked, Mask)
plot(cropped)

# Calculate mean and sd
sdVal <- cellStats(cropped, "sd")
meanVal <- cellStats(cropped, "mean")

s <- cropped
s[which(s[] > (meanVal-sdVal)] <- NA # for example
plot(s)
plot(Mask, add=TRUE)

####

I can also create a stack:

####

#Load libraries
library(raster)

# List files
files <- list.files("...", pattern = "tif", full.names = TRUE)

# Create stack
siteStack <- raster::stack(files)

####

I cannot, however, work out how to iterate through the stack, and save
individual tifs.

One suggestion I've come across is to use a for loop:

####

# List all the files you want to work on
raster.files = dir(“/myfolder/my.raster.folder”)

for (i in 1:length(raster.files)){
                map = raster(raster.files[i]

# do some other stuff …

# Save results as csv or write.raster() to save maps
}

####

I've seen elsewhere that you should avoid for loops and use lapply().

Either way, though I can find plenty of examples that process the
stack using raster::calc() to create a single output file, I can't get
anything to work that sequentially processes each image in isolation.
Stackoverflow, R4DS, the usual websites and blogs, even my
University's internal R help Slack group have all drawn a blank.
(Which I find a bit surprising - this doesn't feel like that
complicated or an unusual thing to do).

I've never written anything like this in R, so any help would be very welcome.

Thanks for your time,

Philip
Part-time PhD student, Lancaster Environment Centre, Lancaster University.