Loop for cropping, cloud masking Sentinel 2 L2A images

Question

I have 6 stacks Sentinel 2 images (10m resolution) that I want to crop to my AOI and cloud masking. After which I want to calculate NDVI in R. The stacks are 3 from different months of 2016 and 3 from 2018.

I have read my AOI and different bands to create a stack as shown below

forest_shp <- readOGR(dsn= ".", layer = "forest_fichtelgebirge1")
 
band2_20160708 <- raster("T32UQA_20160708T101032_B02_10m.jp2")  
band3_20160708 <- raster("T32UQA_20160708T101032_B03_10m.jp2")  
band4_20160708 <- raster("T32UQA_20160708T101032_B04_10m.jp2")  
band8_20160708 <- raster("T32UQA_20160708T101032_B08_10m.jp2")  
all_bands_20160708 <- stack (band2_20160708, band3_20160708, band4_20160708, band8_20160708)

#cropping/masking

new_all_bands_2016070 <- mask(crop(forest_shp), all_bands_20160708)

I seek any looping ideas?

Answer 1

Assuming you have put your stacks into a list stacks_list, looping in R can be done either using the regular for structure:

for (stk in stacks_list) {
# Calculate NDVI here
}

Or with lapply from the apply family:

ndvis = lapply(stacks_list, FUN = function(s) {
# Calculate NDVI for stack s here
})

I would add that if you convert the stacks to a stars object (or read the tiff files as stars objects from the beginning), then you can use the st_apply() function. THis works similarly to any of the *apply family with the addition that it implements parallel processing. So you first load the parallel packages and define your cluster (parallel::makeCluster()), then you call st_apply with the CLUSTER parameter for parallel processing.

clust = makeCluster(num_cores)
ndvis = st_apply(stars_list, 
      CLUSTER = clust,
      MARGIN = ..., # Fill in here
      FUN = function(s) {
      # Calculate NDVI on stars object here
      })

---

Regarding cropping and masking

You will want something like this inside your looping structure (whichever approach you choose), where s is a single stack from the list of stacks stacks_list:

stk_cropped = mask(crop(s, forest_shp), forest_shp))

The crop() function returns a raster or rasterStack from s clipped to the extent of forest_shp. Then the mask() function sets all values outside of forest_shp to NULL. (In case the forest_shp is an irregular polygon)

HTH