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__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
.          .                    .
.          .                    .
.          .                    .
  • It is hard to help as you do not provide example data. I do not think you need a loop. See ?raster::extract Aug 22, 2018 at 12:36
  • I don´t see the benefit of providing sample data. I have hundreds of rasters in TIFF format. Each cell has a xy coordinate and an value, cell size is 30m. I thought that this is standard for raster. Also I dont know how to upload data to this question. Aug 22, 2018 at 12:45
  • I am not suggesting that you provide your data. I am asking to provide an example with data and code so that someone can help you. See the examples in the R help files. For example, start with data in ?raster::stack or ?raster::extract. There are plenty of examples on this site and on stackoverflow. And see: stackoverflow.com/questions/5963269/… and r-project.org/posting-guide.html Aug 22, 2018 at 13:17
  • 1
    @FelixTischer every minute you spend preparing a question here will pay back - if you provide a way for us to work with sample data then it saves us the bother of us all having to create sample data for us to work out the code and make sure it works. Its hard to give a concrete code solution from an abstract question.
    – Spacedman
    Aug 22, 2018 at 13:24
  • 1
    I am very sorry. Thanks for the links to how properly ask questions. I now uploaded some of my data. Aug 22, 2018 at 13:44

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

r <- raster(system.file("external/test.grd", package="raster"))
s <- stack(r, r*2, r/2)
#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

#[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'))

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

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)
#       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.

  • Hey Robert, thanks for your answer. Could you please add comments to each command you are using? I am very new to R and dont know most of the commands you used. I would like to understand the single steps of your code. Aug 22, 2018 at 14:53
  • I think I understood most of your code. But it ends with having the means. How do I get this into my desired output format i posted at the end of my question, table with the time the mean value and the satellite abbreviation? Aug 23, 2018 at 9:21

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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