lapply 在尝试堆叠栅格时创建无限循环答案

【问题标题】：lapply creates endless loop when attempting to stack rasterslapply 在尝试堆叠栅格时创建无限循环
【发布时间】：2020-07-08 01:25:29
【问题描述】：

我正在尝试从我的目录中生成由其他堆叠光栅图像组成的单个光栅图像：

ncfiles <- list.files("~/Desktop/Summer 2020/Tropomi/Aerosol Height", full.names = T, pattern = "*.nc")

当我关注这个example 并运行这个循环时：

bigstack <- stack()

test <- function(file) { for (i in 1: length(ncfiles)){
  GetMyImage <- tryCatch(        
    {
      fname <-(ncfiles[i])
      f <- nc_open(fname)
      print(fname)
    },
    error=function(e) {
      message('Caught Error')
      print(e)
    },
    warning=function(w) {
      message('Caught Warning')
      print(w)
    },
    finally = {
      message('All done')
    }
  )
  if(inherits(errorCondition("ERROR :", next)))
  {
    varx <- attributes(f$var) $names
    vary <- ncvar_get(f, varx)
    rm(f)
    proj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
    rbrick <- brick(vary, crs=proj)
    rm(vary)
    extent(rbrick) <- c(-180, -140, 10, 30)
    return(rbrick)
  }}}

然后尝试堆叠：

allyrs <- lapply(ncfiles, test)

我的循环不断重复运行。是什么导致我的循环无休止地运行？以及如何生成所需的堆叠光栅图像？感谢您的任何意见！

【问题讨论】：

造成问题的一个原因可能是您在 test 函数中运行了一个 for 循环。 lapply 已经将test 应用于ncfiles 中的所有文件，因此无需在test 本身内循环。尝试删除for 语句，并将ncfiles[i] 替换为file（test 函数的参数）
你打算用inherits(errorCondition(...))做什么？首先，inherits 接受两个参数，x=（对象）和what=（1 个或多个字符串，指示与对象类进行比较的类）；你错过了what。其次，errorCondition 的设计目的是始终将某些东西放在类 "error" 之外，这就是它的功能，因此直接将其与inherits 一起使用对我来说毫无意义，即使您确实为inherits 提供了第二个参数。第三，next 是一个保留字，旨在用于某些循环（例如，for、while），而不是作为函数的参数。
如果您的所有栅格都在同一范围内并且 ncfiles 是您的栅格文件列表，您可以通过 do.call(raster::stack, lapply(raster::raster, ncfiles)) 读取。

标签： r for-loop geospatial lapply raster

【解决方案1】：

我建议让事情更加模块化。首先是你想要一个文件的功能

open_one_file <- function(fname) {
    f <- nc_open(fname)
    varx <- attributes(f$var) $names
    vary <- ncvar_get(f, varx)
    proj <- "+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"
    rbrick <- brick(vary, crs=proj)
    extent(rbrick) <- c(-180, -140, 10, 30)
    return(rbrick)
}

现在获取文件名

ncfiles <- list.files(full.names = T, pattern = "\\.nc$")

用几个文件测试函数。例如

test <- open_one_file(ncfiles[3])

如果这适用于某些文件，请使用循环或lapply 来读取它们。显然有些失败，所以我添加了一个try 子句。

x <- lapply(ncfiles, function(i) try(open_one_file(i)))

如果你有失败的文件，使用循环可能更容易

x <- list()
for (i in 1:length(ncfiles)) {
    print(i); print(ncfiles[i]); flush.console()
    x[[i]] <- open_one_file(ncfiles[i]) 
 }

您可以像这样在列表x 中组合 RasterBricks

y <- stack(x)

正如@at80 指出的那样，您更愿意使用标准的 ncdf 文件

x <- lapply(ncfiles, brick)

或者

x <- lapply(ncfiles, function(f) brick(f, var="variable names"))

又来了

y <- stack(x)

【讨论】：

【解决方案2】：

我的问题是使循环过于复杂，并且包含反直觉的组件，就像之前在 cmets 中建议的那样。我将循环重新格式化为：


ncfiles <- list.files("~/Desktop/Summer 2020/Tropomi/Aerosol Height", full.names = T, pattern = "*.nc")

bigstack <- stack()

for (i in 1: length(ncfiles)){
      fname <-(ncfiles[i])
      f <- nc_open(fname)
  ah <- ncvar_get(f, varid = "DETAILED_RESULTS/aerosol_optical_thickness")
  lon <- ncvar_get(nc, varid = "PRODUCT/longitude")
  lat <- ncvar_get(nc, varid = "PRODUCT/latitude")
  nc_close(f)
   s1 <- data.frame(as.vector(lon), as.vector(lat), as.vector(ah))
   crsLatLon <- "+proj=longlat +datum=WGS84"
   ex <- extent(c(-180,180,-90,90))
   pmraster <- raster(ncol=360*10, nrow=180*10, crs=crsLatLon,ext=ex)
   pmraster <- rasterize(s1[,1:2], pmraster, s1[,3], fun=mean, na.rm=T)
   exHI <- extent(c(-180,-140,10,30))
   levelplot(crop(pmraster,exHI))
   bigstack <- stack(bigstack, pmraster)
   print("test")
}

stacking <- calc(bigstack, fun=mean, na.rm=T)
levelplot(crop(stacking, exHI))

【讨论】：