流线型代码，多年的数据处理-应用，栅格，堆栈，平均，裁剪到AOI。

Question

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

function

• Stacks

• 使用lapply将文件放在一个集合文件夹中，例如，1997年数据

• 将文件列表放入砖块中--它们是NetCDF文件，因此，我使用了砖块-砖块-每一个

• 都是一个栅格，每一个

• ，平均从堆栈

• 将新的平均栅格种植到感兴趣区域(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)

非常感谢。

Answer 1

像这样的东西应该能起作用

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)
• If效率更高，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, paste0(year, ".tif"))
    print(fout); flush.console()
    # if (file.exists(fout)) next
    files <- list.files(path=fp, pattern="\\.nc$", full.names=TRUE)
    r <- rast(files)
    s <- mean(r)
    s <- crop(s, AOI, filename=fout) #, overwrite=TRUE)
}