栅格重复随机采样R代码优化：1小时任务缩至15分钟

嘿，我来帮你搞定这个采样性能的问题！你的代码现在跑超1小时，咱们一步步把它提速到15分钟以内～

先搞懂原代码为啥慢 #

看你的代码和性能分析结果，核心问题有这几个：

• 每次调用sample_grid都要重新生成网格多边形（rasterToPolygons），这玩意儿创建Spatial对象的开销贼大，你重复1000次就等于做了1000次无用功
• 用sapply逐个处理每个多边形的采样，循环次数太多，而且Spatial对象的操作本身就慢
• replicate套sapply的嵌套循环，没有利用向量化或者并行计算的优势

针对性优化步骤 #

1. 先把所有网格的边界预计算好，别重复造轮子 #

原代码每次采样都重新生成网格多边形，完全没必要！咱们先把5到25所有网格大小的边界坐标一次性算好，存起来后面直接用，省掉大量Spatial对象创建的时间。

2. 用向量化操作替代逐个多边形采样 #

别再一个个给多边形采样了，咱们直接一次性生成所有需要的随机点坐标，然后批量提取栅格值，这比循环快N倍。

3. 跳过SpatialPolygons，直接用坐标计算 #

其实根本不用把栅格转成多边形！咱们直接算每个网格的x/y范围，在范围内生成随机坐标，然后直接提取值，跳过所有Spatial对象的繁琐操作，这是提速的关键。

4. 用上并行计算，把CPU核心都利用起来 #

1000次重复采样完全可以分给多个CPU核心同时跑，直接把时间砍到原来的1/4甚至更少（取决于你有多少核心）。

优化后的完整代码 #

library(raster)
library(dplyr)
library(parallel)

# 预计算所有网格大小的边界坐标（替代生成多边形）
precompute_grid_bounds <- function(r, grid_sizes) {
  lapply(grid_sizes, function(w) {
    ext <- extent(r)
    # 生成网格的x和y分割点
    x_seq <- seq(ext@xmin, ext@xmax, by = w)
    y_seq <- seq(ext@ymin, ext@ymax, by = w)
    # 生成所有网格的xmin/xmax、ymin/ymax
    expand.grid(xmin = x_seq[-length(x_seq)], ymin = y_seq[-length(y_seq)]) %>%
      mutate(xmax = xmin + w, ymax = ymin + w)
  }) %>% set_names(grid_sizes)
}

# 单个网格大小的快速采样函数
fast_sample_grid <- function(r, grid_bounds, n_reps = 1000) {
  n_cells <- nrow(grid_bounds)
  # 一次性生成所有重复的随机点坐标（每个网格1000个点）
  x_coords <- runif(n_cells * n_reps, grid_bounds$xmin, grid_bounds$xmax)
  y_coords <- runif(n_cells * n_reps, grid_bounds$ymin, grid_bounds$ymax)
  
  # 批量提取所有点的栅格值
  vals <- extract(r, cbind(x_coords, y_coords))
  
  # 整理成1000行（每次重复），求每行的均值（对应原代码的mean）
  matrix(vals, nrow = n_reps, ncol = n_cells) %>%
    rowMeans(na.rm = TRUE)
}

# 主流程
# 创建原栅格
r <- raster(ncol = 50, nrow = 50, xmn = 0, xmx = 50, ymn = 0, ymx = 50)
values(r) <- runif(ncell(r))
grid_sizes <- 5:25

# 预计算所有网格的边界
grid_bounds_list <- precompute_grid_bounds(r, grid_sizes)

# 启动并行计算（用CPU核心数-1，别占满）
cl <- makeCluster(detectCores() - 1)
# 把需要的变量和包传到每个核心
clusterExport(cl, c("r", "fast_sample_grid"))
clusterEvalQ(cl, library(raster))

# 并行处理每个网格大小的采样
results <- parLapply(cl, grid_bounds_list, function(bounds) {
  fast_sample_grid(r, bounds, n_reps = 1000)
})
stopCluster(cl)

# 整理成和原代码一样的输出格式（每行是一次重复，每列是一个网格大小）
results_matrix <- do.call(cbind, results)
colnames(results_matrix) <- grid_sizes

为啥这代码快？ #

• 预计算边界：只生成一次所有网格的边界，避免了1000次重复创建多边形的开销
• 向量化生成坐标：一次性生成所有随机点，替代了逐个多边形的循环采样，速度提升几个数量级
• 并行计算：把不同网格大小的任务分给多个CPU核心，直接把时间砍到原来的1/4左右
• 跳过Spatial对象：完全不用创建SpatialPolygons，避免了性能分析里看到的initialize、getClassDef这些耗时操作

额外小技巧 #

如果你的采样允许一点点误差（比如点刚好落在网格边缘也没关系），可以直接把原栅格聚合到目标网格大小，然后对每个聚合像元随机选一个原像元的值，代码会更简洁，速度还能再快：

# 快速替代方案（近似采样）
fast_approx_sample <- function(r, w, n_reps = 1000) {
  # 聚合栅格
  agg_r <- aggregate(r, fact = w)
  # 每个聚合像元对应的原像元索引
  idx <- Which(agg_r, cells = TRUE)
  # 每次重复随机选每个聚合像元里的一个原像元
  replicate(n_reps, mean(values(r)[sapply(idx, sample, size = 1)], na.rm = TRUE))
}

内容的提问来源于stack exchange，提问作者Ndr