在 R 中拆分光栅迭代或循环

Question

R 编程语言（对此不熟悉）

我正在尝试遍历 splitRaster 输出的多个平铺栅格。在循环期间，我想对每个栅格执行一些处理。

但是下面的代码会报错。

library(ForestTools)
library(raster)
library(sp)
library(rgdal)
library(SpaDES)
rm(list = ls())
tmpdir <- file.path(tempdir(), "splitRaster")
lin <- function(x){x * 0.1 + 0.6}
inCHM <- raster("input raster path and name.tif")
split <- splitRaster(inCHM, 5, 5, c(0.05, 0.05), tmpdir)
files <- list.files(path=tmpdir, pattern="*.grd", full.names=FALSE, recursive=FALSE)
file.names <- dir(tmpdir, pattern ="*.grd")

for(file.names in files ){
name <- file.names
ttops <- vwf(name, winFun = lin, minHeight = 5)
writeOGR(ttops, "output folder", name, driver = "ESRI Shapefile")
}

这是错误

[1] "Xrastername_tile1.grd"
Error in CRS(x) : 
  PROJ4 argument-value pairs must begin with +: Xrastername_tile1.grd

---

更多问题（2020 年 7 月 24 日），

我已经删除了用于故障排除的循环，而只是选择将在循环中使用的 splitRasters 输出之一，即 files[[3]]

当我运行以下代码时，错误是一样的；

library(ForestTools)
library(raster)
library(sp)
library(rgdal)
library(SpaDES)
rm(list = ls())
# set temp directory
tmpdir <- "C:\\R-Test\\Temp_Output"
# get raster
r <- raster("C:\\Lidar\\grid_treeheight_max_1m_nofill.tif")
# define projection
projection(r) <- "+proj=utm +zone=50 +south +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
# split raster brick
y <- splitRaster(r, 8, 8, c(0.05, 0.05), tmpdir)
# Get the complete file locations with full.names = T
files <- list.files(path=tmpdir, pattern="*.grd", full.names=FALSE, recursive=FALSE)
tmpfile <- paste(tmpdir, "\\", files[[3]], sep="")
lin <- function(x){x * 0.06 + 0.6}
ttops <- vwf(tmpfile, winFun = lin, minHeight = 5)

这是错误

Error in CRS(x) : 
  PROJ4 argument-value pairs must begin with +: D:\R-Test\Temp_Output\Xgrid_treeheight_max_1m_nofill_tile11.grd

当我使用上述代码中的一个 splitRaster 输出（文件 [[3]]）运行以下代码时，它运行时没有错误并且我能够绘制 ttops。

rm(list = ls())
# set temp directory
tmpdir <- "D:\\R-Test\\Temp_Output"
# get raster
r <- raster("D:\\R-Test\\Temp_Output\\Xgrid_treeheight_max_1m_nofill_tile11.grd")
lin <- function(x){x * 0.06 + 0.6}
ttops <- vwf(r, winFun = lin, minHeight = 5)

为什么会出现 PROJ4 错误？

这似乎是导致循环失败的错误？

Answer 1

我认为问题在于您尝试使用文件名而不是光栅对象来提供 vwf 函数。我还建议在循环中使用lapply 而不是for。这是一个应该可以工作的代码

library(raster)
library(ForestTools)
library(rgdal)

# Get the complete file locations with full.names = T
files <- list.files(path=tmpdir, pattern="*.grd", full.names=T, recursive=FALSE)

# Loop over each item of the list, i.e., each raster
lapply(files, function(x){
  # Load the image as raster
  image <- raster(x)
  # Calculate vwf (I added a dummy function for winFun)
  ttops <- vwf(image, winFun = function(x){x * 0.06 + 0.5}, minHeight = 5)
  # Write the file with the name of each raster
  writeOGR(ttops, "output_dir", names(x), driver = "ESRI Shapefile")
})