0

我有数百个 Landsat8 场景，每个场景有 12 个波段。我已将它们全部存储在一个文件夹中。现在我尝试在 R 中将它们全部堆叠在一起，以便更轻松地批量处理索引。这是我使用的代码：

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

我认为堆栈函数不知道要使用哪个文件，但我也不知道如何更改它。下面你会看到 Debug 向我展示的内容。我无法真正从中读取任何内容，但该文件不存在。

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)
  }
}
4

1 回答 1

2

当每个陆地卫星图像（及其波段）位于单独的文件夹中时，应该有一个更简单的技巧。然后可以应用以下内容：

list <- list.files(path='C:/ENTERPATHHERE', full.names=TRUE)
inimage <- stack(list)

上面的代码对我有用。如果您应用它但它仍然不适合您，您可能需要检查您的栅格的所有位置是否都正确写入。

希望我能提供帮助！

于 2018-07-28T20:04:27.890 回答