Creating a Loop in R for averaging Raster values to create time series

Question

I have calculated several different indices for lakes stored as GEOTiff. I now want to write a loop with R to average the values for each lake, so only one value is left for each lake. And then create a table with the global means/medians and the date which is given in the filenames.

Input is raster data with following naming convention:

ndci__LC08__192023__20130515.TIFF

ndci__LC08__192023__20130803.TIFF . . . I uploaded five of my rasters here:

https://www.dropbox.com/s/fax1a16tvg0iohd/ndci_LC08_192023_201305151.tif?dl=0 https://www.dropbox.com/s/ghre1290pj0mwnl/ndci_LC08_192023_201404161.tif?dl=0 https://www.dropbox.com/s/cu1x3kxnd8koz2d/ndci_LC08_192023_201503181.tif?dl=0 https://www.dropbox.com/s/s2uqsj6f14b5u06/ndci_LC08_192023_201508091.tif?dl=0 https://www.dropbox.com/s/044b5vll1rsrd8a/ndci_LC08_192023_2013080321.tif?dl=0

So the output format I am looking for is:

Date       Mean/Median of Index Satellite
20130515   5.4                  LC08
20130803   8.7                  LC08
.          .                    .
.          .                    .
.          .                    .

Answer 1

Instead of your own data, you should try to use example data created by code or from R examples. That avoids unnecessary work, and makes the answers more useful --- they will be easier to use, also after those links go dead.

Example data:

library(raster)
r <- raster(system.file("external/test.grd", package="raster"))
s <- stack(r, r*2, r/2)
s
#class       : RasterStack 
#dimensions  : 115, 80, 9200, 3  (nrow, ncol, ncell, nlayers)
#resolution  : 40, 40  (x, y)
#extent      : 178400, 181600, 329400, 334000  (xmin, xmax, ymin, ymax)
#coord. ref. : +init=epsg:28992 +towgs84=565.237,50.0087,465.658,-0.406857,0.350733,-1.87035,4.0812 +proj=sterea +lat_0=52.15616055555555 +lon_0=5.38763888888889 +k=0.9999079 +x_0=155000 +y_0=463000 +ellps=bessel +units=m +no_defs 
#names       :  test.1,  test.2,  test.3 
#min values  : 128.434, 256.868,  64.217 
#max values  : 1805.78, 3611.56,  902.89

names(s)
#[1] "test.1" "test.2" "test.3"


xy <- matrix(c(179570, 179682, 179347, 179384, 179570, 331252, 330945, 330945, 331290, 331252), ncol=2)
lakes <- spPolygons(xy, xy+1000, df=data.frame(name=c('A', 'B'))
lakes

#class       : SpatialPolygons 
#features    : 2 
#extent      : 179347, 180682, 330945, 332290  (xmin, xmax, ymin, ymax)
#coord. ref. : NA 

plot(r)
plot(lakes, add=TRUE, col='blue')

So we have three rasters with values of interest, and two lakes for which we want to summarize the data. We can do:

y <- extract(s, lakes, fun=mean)
y
#       test.1    test.2   test.3
#[1,] 536.9039 1073.8077 268.4519
#[2,] 346.0194  692.0389 173.0097

apply(y, 1, mean)
#[1] 626.3878 403.6893

Two values, one for each lake.