简化代码，需要处理多年的数据 - lapply、raster、stack、mean、crop to AOI

编程问答 2022-05-03

问题描述

我在 R 中有一些代码，它执行以下操作：

  1. 使用 lapply 将文件导入设置文件夹中，例如1997年数据
  2. 将文件列表变成砖块 - 它们是 NetCDF 文件，所以我使用了砖块功能
  3. 将砖块按每一年的月份堆叠成一个栅格堆栈。
  4. 从堆栈中计算平均值
  5. 将新的平均栅格裁剪到感兴趣区域 (AOI)。

我有一个工作代码，见下文，但它很笨重，我觉得在一个循环中可以更好然后遍历每年的文件夹（我有 1997 年到 2018 年的数据） .任何人都可以帮助将其简化为我可以通过更改文件路径运行的简单循环代码吗？我以前使用过一些循环，但不是从头开始的。

# Packages:
library(raster)
library(parallel) # Check cores in PC
library(lubridate) # needed for lapply
library(dplyr) # ""
library(sf) # For clipping data
library(rgdal)

# ChlA
  
# Set file paths for input and outputs:
usingfp <- "/filepath/GIS/ChlA/1997/"
the_dir_ex <- "Data/CHL/1997"

# List all NETCDF files in folder:
CHL_1997 <- list.files(path = usingfp,pattern = "\\.nc$",full.names = TRUE,recursive = FALSE)

# Make file list into brick
CHL_1997_brick <- lapply(CHL_1997,FUN = brick,the_dir = the_dir_ex)
# Stack bricks
s <- stack(CHL_1997_brick)

# Calculate mean from stack
mean <- calc(s,fun = mean,na.rm = T)
plot(mean)

# Load vector boundary to "crop" to
AOI <- readOGR("/filepath/AOI/AOI.shp")
plot(AOI,main = "Shapefile imported into R - crop extent",axes = TRUE,border = "blue",add = T)

# crop the raster using the vector extent
CHL_1997_mean <- crop(mean,AOI)
plot(CHL_1997_mean,main = "Cropped mean CHL - 1997")

# add shapefile on top of the existing raster
plot(AOI,add = TRUE)

非常感谢。

解决方法

这样的东西应该可以工作

library(raster)

AOI <- shapefile("/filepath/AOI/AOI.shp")
path <- "/filepath/GIS/ChlA/"
years <- 1997:2018

for (yr in years) {
    fp <- file.path(path,yr)
    fout <- file.path(fp,paste0(year,".tif"))
    print(fout); flush.console()
    # if (file.exists(fout)) next
    files <- list.files(path=fp,pattern="\\.nc$",full.names=TRUE)
    b <- lapply(files,brick)
    s <- stack(b)
    s <- mean(s)
    s <- crop(s,AOI,filename=fout) #,overwrite=TRUE)
}

注意事项：

  • mean(s) 比 calc(s,mean)
    • 更有效率
  • 如果 AOI 相对较小，则首先执行会更有效 使用 crop，然后使用 mean（然后使用 writeRaster）

您也可以像这样使用 terra：

library(terra)
AOI <- vect("/filepath/AOI/AOI.shp")
path <- "/filepath/GIS/ChlA/"
years <- 1997:2018

for (yr in years) {
    fp <- file.path(path,year)
    fout <- file.path(fp,full.names=TRUE)
    r <- rast(files)
    s <- mean(r)
    s <- crop(s,overwrite=TRUE)
}
loopsnetcdfrrr-rasterspatial