1

I have 111 Landsat8 scenes which I want to calculate several indices for. The algebra to calculate them is most of the time Index= (b1-b2)/(b1+b2). I have the single bands (bands 1-11 respectively) all stored in one folder, named LC08_L1TP_194023_20140225_20170425_01_T1_B1 etc.

I am looking for a good way to automatically calculate the indices.

I have made an approach in R already, which works fine for one file, but I don't know how to make a loop reading the correct files in the folder out of it.

red=raster("Path:/LC08_L1TP_194024_20171015_20171024_01_T1_B4.tif")
nearinfrared=raster("Path:/LC08_L1TP_194024_20171015_20171024_01_T1_B5.tif")

ndvi = (nearinfrared - red) / (nearinfrared + red)
writeRaster(ndvi, filename="Path/NDVI_Test",format="GTiff")    
  • 1
    Please poke around on the site, this is answered in so many different ways on a huge variety of data, including landsat, it is difficult to select a duplicate post. You can pull R specific tags using [R] in the search term. – Jeffrey Evans Jul 27 '18 at 15:17
0

Your question is too broad to post a specific answer, but the basic procedure is:

List files

library(raster)
green <- list.files(path = your_path, pattern = 'B3.tif$', full.names=T)
red <- list.files(path = your_path, pattern = 'B4.tif$', full.names=T)
nir <- list.files(path = your_path, pattern = 'B5.tif$', full.names=T)
# etc

Create functions for each index

norm_diff <- function(x,y){(y-x)/(y+x)}
other_indx <- funcion(x,y){((y-x)/(y+x+0.5))*1.5}

Iterate through rasters

for(i in seq_along(red)){
  ndvi <- overlay(x = raster(red[i]), y = raster(nir[i]), fun = norm_diff)
  gndvi <- overlay(x = raster(green[i]), y = raster(nir[i]), fun = norm_diff)
  savi <- overlay(x = raster(red[i]), y = raster(nir[i]), fun = other_indx)

  writeRaster(ndvi,filename = paste0(your_path, 'ndvi_',i,'.tif'))
  writeRaster(gndvi,filename = paste0(your_path, 'gndvi_',i,'.tif'))
  writeRaster(savi,filename = paste0(your_path, 'savi_',i,'.tif'))
  removeTmpFiles(h=0.1)
}

This procedure is really basic. If your files are sorted by date, this can cover all your needs. There are more packages or functions to optimize this, but they are more complex to understand.

  • 1
    wow. I tried to get Shelbys Solution to run for the last two days. I should have sticked with your Code though. Its more simple so cavemen like me cant put that many mistakes in the code. Thank you very much Aldo! – Felix Tischer Jul 29 '18 at 11:36
3

The first step I would take is to create a list of all of the prefixes for the images, which will help you to iterate through them all. Using the directory level that contains all of the folders which contain the individual bands:

#list of files to use for reference, band1 picked arbitrarily 
L8files = list.files("L8", pattern = "band1", full.names=TRUE)

#get prefixes 
getprefix = function(string){
 substr(string, 4, 51) 
} 
L8list = lapply(L8files, getprefix)

Once you have a list of every prefix, the code below will allow you to iterate through the list. For each image, you're going to want to stack all of the bands. A raster stack is a single file, but all the bands are still kept separate within the image.

This is optional, but you might want to create another folder to store all of the stacked raster images and then work out of that:

#create folder for processed images
dir.create("stacked")

Here's the code to stack them all:

stack= function(file){
 setwd("C:/directory_with_image_folders")
 prefix = substr(file, 1, 48)
 suffix = "tif"

 inband1 = raster(paste (prefix, paste("1", suffix, sep ="."), sep="")) #blue
 inband2 = raster(paste (prefix, paste("2", suffix, sep ="."), sep="")) #green
 inband3 = raster(paste (prefix, paste("3", suffix, sep ="."), sep="")) #red
 inband4 = raster(paste (prefix, paste("4", suffix, sep ="."), sep="")) #nir
 inband5 = raster(paste (prefix, paste("5", suffix, sep ="."), sep="")) #swir1
 inband6 = raster(paste (prefix, paste("6", suffix, sep ="."), sep="")) #swir2

 #stack bands
 inimage = stack(inband1, inband2, inband3, inband4, inband5, inband6)

 setwd("C:/stacked")
 sat = substr(file, 1, 4)
 date = substr(file, 18, 25)
 writeRaster(inimage, filename= paste(date, sat, sep="_"), format="GTiff", overwrite=TRUE)
}

for (i in L8list){
  stack(i)
}

Finally, you can iterate through the raster stack images to calculate the NDVI for each. Once again, I recommend a new directory for the NDVI images, but this is optional:

#create folder & list 
dir.create("ndvi")

L8list = list.files("stacked", pattern = "LC08", full.names=TRUE)

#ndvi function
ndvi=function(i){
 setwd("C:/stacked")
 stack = stack(i)
 ndvi = (stack[[4]] - stack[[3]]) / (stack[[4]] + stack[[3]])

 #save output
 setwd("C:/ndvi")
 writeRaster(ndvi, filename= paste(substr(i, 11, 23), "ndvi", sep="_"), format="GTiff", overwrite=TRUE)

And there you go! That should allow you to easily go through and create NDVI images out of each image you have. You can also go through and alter any of the naming conventions that I've used here.

  • I am confused by this line: L8files = list.files("L8", pattern = "band1", full.names=TRUE) , I made L8files = list.files(path = "C:/Users/Felix/Desktop/Bachelorarbeit/Daten/R_Test/Only_TIF", full.names=TRUE) out of it. – Felix Tischer Jul 28 '18 at 23:30
  • I am probably just to unexperienced with R to get your code going. Though i learned a lot trying to figure out whats wrong. Thank you very much for your solution Shelby. If I could i would mark both answers as a solution since i am combining them. Aldo´s stacking and loop works for me an your naming tricks give me the right output. – Felix Tischer Jul 29 '18 at 11:52

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.