I have several hundred Landsat8 scenes with 12 bands each. I have stored them all in one folder. And now i try to stack them all together in R for easier batch processing of indices. This is the code im using:
#Create List of files (whole path is used)
L8files = list.files(path =
"C:/Users/Felix/Desktop/Bachelorarbeit/Daten/R_Test/Only_TIF",
full.names=TRUE)
#Write Names of Files in List
getprefix = function(string){
substr(string, 61, 113)
}
L8list = lapply(L8files, getprefix)
#Defining function to stack
stack_raster= function(file){
setwd("C:/Users/Felix/Desktop/Bachelorarbeit/Daten/R_Test/Only_TIF")
prefix = substr(file, 1, 106)
suffix = "tif"
inband1 = raster(paste (prefix, paste("B1", suffix, sep ="."), sep=""))#coastal aerosol
inband2 = raster(paste (prefix, paste("B2", suffix, sep ="."), sep="")) #blue
inband3 = raster(paste (prefix, paste("B3", suffix, sep ="."), sep="")) #green
inband4 = raster(paste (prefix, paste("B4", suffix, sep ="."), sep="")) #red
inband5 = raster(paste (prefix, paste("B5", suffix, sep ="."), sep="")) #NIR
inband6 = raster(paste (prefix, paste("B6", suffix, sep ="."), sep="")) #swir1
inband7 = raster(paste (prefix, paste("B7", suffix, sep ="."), sep="")) #swir2
inband8 = raster(paste (prefix, paste("B8", suffix, sep ="."), sep="")) #PAN
inband9 = raster(paste (prefix, paste("B9", suffix, sep ="."), sep="")) #cirrus
inband10 = raster(paste (prefix, paste("B10", suffix, sep ="."), sep="")) #TIRS1
inband11 = raster(paste (prefix, paste("B11", suffix, sep ="."), sep="")) #TIRS2
inbandQBA = raster(paste (prefix, paste("QBA", suffix, sep ="."), sep=""))#Pre-Collection Quality Assessment
#stack bands
inimage = stack(inband1, inband2, inband3, inband4, inband5, inband6,
inband7, inband8,inband9, inband10, inband11, inbandQBA)
setwd("C:/Users/Felix/Desktop/Bachelorarbeit/Daten/R_Test/stacked")
sat = substr(file, 67, 70)
date = substr(file, 84, 91)
writeRaster(inimage, filename = paste(date, sat, sep="_"),
format="GTiff", overwrite=TRUE)
}
#Loop
for (i in L8list){
stack_raster(i)
}
I think the stack function doesn´t know which file to use, but i also dont know how to change that. Below you see what Debug shows me. I cant really read anything out of it but that the file doesnt exist.
function (x, band = 1, objecttype = "RasterLayer", native = FALSE,
silent = TRUE, offset = NULL, ncdf = FALSE, ...)
{
x <- trim(x)
if (x == "" | x == ".") {
stop("provide a valid filename")
}
start <- tolower(substr(x, 1, 3))
if (!start %in% c("htt", "ftp")) {
y <- NULL
try(y <- normalizePath(x, mustWork = TRUE), silent = TRUE)
if (!is.null(y)) {
x <- y
}
}
fileext <- toupper(extension(x))
if (fileext %in% c(".GRD", ".GRI")) {
grifile <- .setFileExtensionValues(x, "raster")
grdfile <- .setFileExtensionHeader(x, "raster")
if (file.exists(grdfile) & file.exists(grifile)) {
return(.rasterFromRasterFile(grdfile, band = band,
objecttype, ...))
}
}
if (!file.exists(x)) {
if (extension(x) == "") {
grifile <- .setFileExtensionValues(x, "raster")
grdfile <- .setFileExtensionHeader(x, "raster")
if (file.exists(grdfile) & file.exists(grifile)) {
return(.rasterFromRasterFile(grdfile, band = band,
objecttype, ...))
}
else {
}
}
}
if ((fileext %in% c(".HE5", ".NC", ".NCF", ".NC4", ".CDF",
".NCDF", ".NETCDF")) | (isTRUE(ncdf))) {
return(.rasterObjectFromCDF(x, type = objecttype, band = band,
...))
}
if (fileext == ".GRD") {
if (.isNetCDF(x)) {
return(.rasterObjectFromCDF(x, type = objecttype,
band = band, ...))
}
}
if (fileext == ".BIG" | fileext == ".BRD") {
return(.rasterFromRasterFile(x, band = band, objecttype,
driver = "big.matrix", ...))
}
if (!is.null(offset)) {
return(.rasterFromASCIIFile(x, offset, ...))
}
if (fileext %in% c(".BIN")) {
r <- .rasterFromNSIDCFile(x)
if (!is.null(r))
return(r)
}
if (!native) {
if (!.requireRgdal(FALSE)) {
native <- TRUE
}
}
if (native) {
if (fileext == ".ASC") {
return(.rasterFromASCIIFile(x, ...))
}
if (fileext %in% c(".BIL", ".BIP", ".BSQ")) {
return(.rasterFromGenericFile(x, type = objecttype,
...))
}
if (fileext %in% c(".RST", ".RDC")) {
return(.rasterFromIDRISIFile(x, ...))
}
if (fileext %in% c(".DOC", ".IMG")) {
return(.rasterFromIDRISIFile(x, old = TRUE, ...))
}
if (fileext %in% c(".SGRD", ".SDAT")) {
return(.rasterFromSAGAFile(x, ...))
}
}
if (fileext == ".DOC") {
if (file.exists(extension(x, ".img"))) {
return(.rasterFromIDRISIFile(x, old = TRUE, ...))
}
}
if (fileext %in% c(".SGRD", ".SDAT")) {
r <- .rasterFromSAGAFile(x, ...)
if (r@file@toptobottom | r@data@gain != 1) {
return(r)
}
}
if (!.requireRgdal(FALSE)) {
stop("Cannot create RasterLayer object from this file; perhaps you need
to install rgdal first")
}
test <- try(r <- .rasterFromGDAL(x, band = band, objecttype,
...), silent = silent)
if (class(test) == "try-error") {
if (!file.exists(x)) {
stop("Cannot create a RasterLayer object from this file. (file does not
exist)")
}
stop("Cannot create a RasterLayer object from this file.")
}
else {
return(r)
}
}
There is an easier trick that should work when each landsat image (with their bands) is in a separate folder. Then the following can be applied:
list <- list.files(path='C:/ENTERPATHHERE', full.names=TRUE)
inimage <- stack(list)
The above code works for me. If you apply it and it still doesn't for you you might want to check whether all the locations of your rasters are written correctly.
Hopefully I was able to help!