优化多时段栅格按多边形提取均值的方法

我有覆盖美国本土的多个逐小时栅格图层，以及一个多边形矢量图层，目前使用terra::extract提取每个多边形的栅格均值。该方法在多边形和栅格数量较少时运行正常，但扩展到数月数据及1000+个多边形时，速度无法满足需求，急需更优化的实现方案。

我已构思三种优化方向，但尚未实现：

• 并行处理，但不清楚如何避免序列化问题；
• 不再为每个栅格单独调用extract(..., fun = mean)，而是通过extract(mrast[[1]], mpoly, cells = TRUE, exact = TRUE)获取每个多边形覆盖的栅格单元及精确覆盖比例，利用所有栅格单元一致的特性，对各时段栅格提取对应单元值后加权计算均值；
• 结合上述两种方法，这应该是最快的方案。

最小可复现示例代码 #

# setwd("C:/temp")
library(terra)
library(R.utils)

murl <- "https://mtarchive.geol.iastate.edu/2023/06/27/mrms/ncep/MultiSensor_QPE_01H_Pass2/"
mfiles <- c("MultiSensor_QPE_01H_Pass2_00.00_20230627-070000.grib2.gz", 
           "MultiSensor_QPE_01H_Pass2_00.00_20230627-080000.grib2.gz")

## 下载、解压并读取栅格文件到R
lapply(seq_along(mfiles), \(i) download.file(paste0(murl, mfiles[i]), mfiles[i])) 
lapply(mfiles, \(f) gunzip(f, gsub(".gz", "", f), remove = TRUE, overwrite=TRUE))

mrast <- lapply(list.files(pattern = ".grib2$"), \(f) rast(f)) 
lapply(seq_along(mrast), \(i) names(mrast[[i]]) <<- time(mrast[[i]]))

## 创建多边形矢量文件
mpoly <- vect(dptply, "polygon") ## dptply数据见下文
crs(mpoly) <- "EPSG:3857"
mpoly <- project(mpoly, "+proj=longlat +datum=WGS84")

## 提取均值（此部分需优化）
startTime = Sys.time()
mavg <- lapply(mrast, \(r) {op <- extract(r, mpoly, 
                                         weights=TRUE, fun=mean, na.rm = TRUE)
                            message(time(r))
                            return(op)})
endTime = Sys.time()
print(endTime - startTime)

# 运行输出：
# 2023-06-27 07:00:00
# 2023-06-27 08:00:00
# Time difference of 10.87896 secs

预期输出示例 #

## 仅作演示
cbind(mavg[[1]][1], do.call(cbind, lapply(mavg, "[", 2)))
#>    ID 2023-06-27 07:00:00 2023-06-27 08:00:00
#> 1   1            0.000000           10.405161
#> 2   2            0.000000            8.961424
#> 3   3            0.000000            6.300000
#> 4   4            0.000000            5.902914
#> 5   5            0.000000            5.462630
#> 6   6            0.000000            4.100000
#> 7   7            0.000000            3.566511
#> 8   8            0.000000           13.275161
#> 9   9            2.600000           11.200000
#> 10 10            4.633752           16.301403

多边形数据 #

dptply <- structure(c(1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 
                      4, 4, 4, 4, 4, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 8, 
                      8, 8, 8, 8, 8, 9, 9, 9, 9, 9, 10, 10, 10, 10, 10, 1, 1, 1, 1, 
                      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                      1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 
                      1, 1, 1, 1, 1, 1, 1, -7945530.9815, -7938378.0621, -7938605.1389, 
                      -7945190.3663, -7945530.9815, -7938264.5237, -7936674.9861, -7930543.9123, 
                      -7932360.5267, -7936674.9861, -7938264.5237, -7938264.5237, -7934698.4679, 
                      -7933886.647, -7933878.7652, -7934666.9409, -7934698.4679, -7933735.8469, 
                      -7932622.0601, -7932634.9736, -7933745.532, -7933735.8469, -7934646.4483, 
                      -7933925.1099, -7933904.9326, -7934631.3153, -7934646.4483, -7933551.8299, 
                      -7932850.6689, -7932865.8018, -7933496.3424, -7933551.8299, -7931645.765, 
                      -7931670.3955, -7932384.6797, -7932470.8864, -7931645.765, -7940928.6868, 
                      -7941271.2217, -7941328.3108, -7948978.2562, -7949434.9693, -7940928.6868, 
                      -7998760.8958, -7999220.4225, -7998966.0077, -7998636.765, -7998760.8958, 
                      -8001395.2248, -8001358.6877, -8003002.8551, -8003149.0033, -8001395.2248, 
                      5262480.6428, 5262635.0714, 5251676.934, 5251831.1854, 5262480.6428, 
                      5262789.5025, 5268041.6687, 5261862.9539, 5252293.955, 5248438.2361, 
                      5248900.8391, 5262789.5025, 5269309.6397, 5269309.6397, 5269052.1707, 
                      5269041.443, 5269309.6397, 5269446.29, 5269441.8957, 5267934.7856, 
                      5267925.9984, 5269446.29, 5271218.9923, 5271198.3903, 5270580.3524, 
                      5270607.8199, 5271218.9923, 5272036.2399, 5272049.9757, 5271507.4245, 
                      5271239.5943, 5272036.2399, 5272266.345, 5270589.7139, 5270556.1843, 
                      5272216.0417, 5272266.345, 5246883.428, 5237972.8, 5237895.3531, 
                      5238282.5936, 5247503.6111, 5246883.428, 5227512.1703, 5227597.3282, 
                      5228266.5923, 5228286.8738, 5227512.1703, 5232252.738, 5229875.3018, 
                      5229528.6423, 5232104.1307, 5232252.738, 0, 0, 0, 0, 0, 0, 0, 
                      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                      0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
                      0, 0, 0, 0), 
                    dim = c(53L, 5L), 
                    dimnames = list(NULL, c("geom", "part", "x", "y", "hole")))

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