大型栅格的并行处理（拆分为多个部分，运行函数，然后重新组合）

Question

总体目标是在堆叠栅格上运行相关分析，其中我使用移动窗口并将焦点像素堆栈与移动窗口中每个其他像素的堆栈相关联。最终结果是一个栅格，它是该像素与周围像素的平均相关性。但我这样做的规模如此之大，以至于处理时间非常长。我需要分成几个部分并在每个部分上运行此函数，然后重新组合。由于我正在做一个移动窗口，我还想将这些部分重叠在移动窗口的范围内，所以我对栅格的大部分没有边缘效应。

我在下面编写了示例堆叠栅格和移动窗口相关函数的代码。我需要有关如何分离光栅然后并行运行多个部分的帮助。

加载包

library(raster)
library(dplyr)

第一步是制作一个随机的 100x100x5 堆叠栅格。我创建的函数使用数组，因此也将栅格更改为数组。

# Set the number of rows and columns
num_rows <- 100
num_cols <- 100

# Create an empty raster
empty_raster <- raster(nrows = num_rows, ncols = num_cols)

# Define the number of layers (stack)
num_layers <- 5

# Create a stacked raster with random values for each layer
stacked_raster <- stack(lapply(1:num_layers, function(i) setValues(empty_raster, rpois(ncell(empty_raster), lambda = 10))))

#make data into an array
data <- raster::as.array(stacked_raster)

设置一些功能参数

y <- dim(data)[3] #number of years of data for correlation



m <- 5 #number of nearest cells you want in each direction 

#create matrix with label for each cell
cell <- matrix(1:10000, 100, 100, byrow=F)

#create distance matrices
lon <- matrix(rep(1:100, 100), nrow=100, ncol=100, byrow = T)
lat <- matrix(rep(1:100, 100), nrow=100, ncol=100, byrow = F)

为 forloop 和函数创建做好准备

#get ready for for loop

#create matrix to hold average correlation value for each cell
avg_cor_values <- matrix(NA, 100, 100, byrow=F)

#create function to do the necessary stuff
corr_function <- function(i){
  focal <- which(cell==i, arr.ind=T)
  
  #get min and max lat and lon for desired grid
  min_lon <- focal[1] - m
  #set min lon to 1 if calculation is negative
  if(min_lon < 1){
    min_lon <- 1
  }
  
  max_lon <- focal[1] + m
  #set max lon to nrow(lon) if calculation is more than nrow(lon)
  if(max_lon > nrow(lon)){
    max_lon <- nrow(lon)
  }
  
  min_lat <- focal[2] - m
  #set min lat to 1 if calculation is negative
  if(min_lat < 1){
    min_lat <- 1
  }
  
  max_lat <- focal[2] + m
  #set max lat to ncol(lat) if calculation is more than ncol(lat)
  if(max_lat > ncol(lat)){
    max_lat <- ncol(lat)
  }
  
  #create dataframe with position of each cell that needs correlation calculation
  focal_grid <- expand.grid(lon=min_lon:max_lon, lat=min_lat:max_lat)
  
  #determine which row in focal_grid matches the focal cell location
  drop_row1 <- which(focal_grid[1:nrow(focal_grid),1]==focal[1], arr.ind=T)  
  drop_row2 <- which(focal_grid[1:nrow(focal_grid),2]==focal[2], arr.ind=T)
  #remove the focal cell from the grid
  focal_grid <- focal_grid[-drop_row1[which(drop_row1 %in% drop_row2)],]
  
  #calculate correlation between focal cell and all cells within focal grid
  cor_results <- NA
  for(n in 1:nrow(focal_grid)){
    cor_results[n] <- cor(data[focal[1], focal[2], 1:y], data[focal_grid[n,1], focal_grid[n,2], 1:y], use = "pairwise.complete.obs") 
  }
  
  avg_cor_values[focal[1], focal[2]] <- mean(cor_results, na.rm=TRUE)
}

在栅格上运行该函数

corr_function_vec <- Vectorize(corr_function)

#run code for all cells in grid
Sys.time()
product <- corr_function_vec(1:max(cell))
Sys.time()

corr_product <- matrix(product, 100, 100, byrow=F)
corr_raster <- raster(corr_product)
plot(corr_raster)

Answer 1

如果您在类 Unix 系统（Linux 或 macOS）上运行，则以下内容可以使用您的函数并在可用内核上加速：

#run serial code for all cells in grid
system.time({
  product <- corr_function_vec(1:max(cell))
})

#run parallel code for all cells in grid
system.time({
  cores = parallel::detectCores()
  product2 <- parallel::mclapply(1:max(cell), corr_function, mc.cores = cores)
  product2 <- do.call(c, product2)
})
all.equal(product, product2)

Vectorize

函数已经在幕后使用了

mapply

，因此这只是将其替换为

parallel

包中的mclapply。您可以跳过

mc.cores

参数，因为我无论如何都提供默认值（为了清楚起见）。

如果您想移植到 Windows，futures 包与 furrr 一起工作也可以，但在 Windows 上，您将在每个并行实例中复制

cell

数组（因为 Unix

fork

位于

mclapply

不适用于共享内存并行计算）。

如果您有一个集群并希望将

cell

数组拆分到节点上，还有更多选项，但由于需要重叠，这里会变得更加复杂。