1) Commenting on scripts is extremely useful, especially if you want to look at it later and understand what’s going on. Test your understanding of the script for pre-processing a single L8 scene shown in the previous section. Copy the following code into the script window of R-Studio. For each line, try to give a brief explanation of what is being done in the next line of code! (Explanation of the first line is given)
# load the raster library in the current session
# (needed for rasterstack and writeRaster)
library(raster)
#
product <- "/media/sf_exchange/landsatdata/LC081930232017060201T1-SC20180613160412.tar.gz"
#
productname <- substr(product, 1, nchar(product) - 7)
#
untar(product, exdir = productname)
#
productfiles <- list.files(productname, full.names = TRUE)
#
bands <- c(grep('band1', productfiles, value=TRUE),
grep('band2', productfiles, value=TRUE),
grep('band3', productfiles, value=TRUE),
grep('band4', productfiles, value=TRUE),
grep('band5', productfiles, value=TRUE),
grep('band6', productfiles, value=TRUE),
grep('band7', productfiles, value=TRUE)
)
#
rasterstack <- stack( bands )
#
writeRaster(rasterstack,
paste0(productname, ".tif"),
format = "GTiff",
overwrite = TRUE
)
#
unlink(productname, recursive=TRUE)
#
makeOVR <- paste0("gdaladdo -ro ", productname, ".tif 2 4 8 16")
system( makeOVR )
Done? Compare it with the solution:
# load the raster library in the current session
# (needed for rasterstack and writeRaster)
library(raster)
# define file path of a single L8 dataset
product <- "/media/sf_exchange/landsatdata/LC081930232017060201T1-SC20180613160412.tar.gz"
# get the export folder for decompression by deleting suffix (.tar.gz)
productname <- substr(product, 1, nchar(product) - 7)
# decompress and save to export folder
untar(product, exdir = productname)
# create a vector of all names of files in decompressed folder
productfiles <- list.files(productname, full.names = TRUE)
# search for individual bands inside file vector
bands <- c(grep('band1', productfiles, value=TRUE),
grep('band2', productfiles, value=TRUE),
grep('band3', productfiles, value=TRUE),
grep('band4', productfiles, value=TRUE),
grep('band5', productfiles, value=TRUE),
grep('band6', productfiles, value=TRUE),
grep('band7', productfiles, value=TRUE)
)
# create a rasterstack of the band selection
rasterstack <- stack( bands )
# write the rasterstack to hard drive
writeRaster(rasterstack,
paste0(productname, ".tif"),
format = "GTiff",
overwrite = TRUE
)
# delete the uncompressed folder in order to save disc space
unlink(productname, recursive=TRUE)
# create pyramid layers for faster visualization
makeOVR <- paste0("gdaladdo -ro ", productname, ".tif 2 4 8 16")
system( makeOVR )
2) Download any Landsat-8 Level-1 dataset using the USGS EarthExplorer. Modify the code given in the previous task so that you create and store a datastack of all the spectral channels available in this Level-1 dataset!
library(raster)
product <- "/media/sf_exchange/landsatdata/LC081930232017060201T1-SC20180613160412.tar.gz"
productname <- substr(product, 1, nchar(product) - 7)
untar(product, exdir = productname)
productfiles <- list.files(productname, full.names = TRUE)
bands <- c(grep('band1', productfiles, value=TRUE),
grep('band2', productfiles, value=TRUE),
grep('band3', productfiles, value=TRUE),
grep('band4', productfiles, value=TRUE),
grep('band5', productfiles, value=TRUE),
grep('band6', productfiles, value=TRUE),
grep('band7', productfiles, value=TRUE),
grep('band8', productfiles, value=TRUE),
grep('band9', productfiles, value=TRUE),
grep('band10', productfiles, value=TRUE),
grep('band11', productfiles, value=TRUE)
)
rasterstack <- stack( bands )
writeRaster(rasterstack,
paste0(productname, ".tif"),
format = "GTiff",
overwrite = TRUE
)
unlink(productname, recursive=TRUE)
makeOVR <- paste0("gdaladdo -ro ", productname, ".tif 2 4 8 16")
system( makeOVR )