R中sf简单要素对象列表调用mapply运行报错问题求助
需要对接收两个sf(简单要素)对象作为输入的自定义函数调用mapply实现并行空间计算,相关可复现代码与数据已公开发布在开源代码仓库。
当前业务逻辑:将点生成的圆形缓冲区shapefile与共71万条记录的加利福尼亚州人口普查区块shapefile执行空间相交计算,统计圆形对每个普查区块的面积重叠占比,仅保留重叠率≥50%的区块。由于该计算必须使用精度更高但运算效率较低的st_intersection命令,全量计算开销极高,因此计划在HPC集群上将任务拆分为县级粒度的并行子任务,缩短整体运行时间。
county_func <- function(x, y) { # (x = circles, y = centroids/points) y_id <- y %>% # keep points only as an index for later st_drop_geometry() x_int <- st_intersects(x, blocks) # intersect circles with blocks to build index of overlap ints_holder <- data.frame() for(i in 1:nrow(x_int)){ # for each circle blocks_int <- subset(blocks, as.numeric(rownames(blocks)) %in% x_int[[i]]) # subset blocks that intersect blocks_int$hud_id <- y_id[i, 1] # add point id to these blocks for merge later ints_holder <- rbind(ints_holder, blocks_int) } x_blocks <- st_intersection(x, ints_holder) %>% # intersection between circles and blocks they intersect mutate(intersect_area = st_area(.)) %>% # calculate intersect area dplyr::select(GEOID10, intersect_area) # drop irrelevant data return(x_blocks) }
单条命令分步运行、传入两个独立sf对象作为输入时,上述函数可正常执行。但适配parallel库的mcmapply/mapply并行函数时出现兼容性问题:将圆形shapefile按县分组拆分为sf对象列表后,即使仅传入列表中单个县的子集运行函数,也会触发报错。
Error in UseMethod("st_drop_geometry") : no applicable method for 'st_drop_geometry' applied to an object of class "character"
该报错在移除输入对象所有字符型变量后仍会复现。已知apply类函数存在强制将传入对象转换为矩阵格式的机制,可能与sf对象存在兼容问题,但暂未找到除parallel库mcmapply系列外、适配HPC集群的其他并行方案。
R version 4.1.0 (2021-05-18) Platform: x86_64-apple-darwin17.0 (64-bit) Running under: macOS Catalina 10.15.7 Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8 attached base packages: [1] parallel stats graphics grDevices utils datasets methods base other attached packages: [1] dplyr_1.0.9 sf_1.0-7 loaded via a namespace (and not attached): [1] Rcpp_1.0.8.3 rstudioapi_0.13 magrittr_2.0.3 units_0.7-2 tidyselect_1.1.2 R6_2.5.1 [7] rlang_1.0.2 fansi_1.0.3 s2_1.0.7 stringr_1.4.0 wk_0.5.0 tools_4.1.0 [13] grid_4.1.0 KernSmooth_2.23-20 utf8_1.2.2 cli_3.3.0 e1071_1.7-9 DBI_1.1.1 [19] ellipsis_0.3.2 class_7.3-19 assertthat_0.2.1 tibble_3.1.7 lifecycle_1.0.1 crayon_1.5.1 [25] purrr_0.3.4 vctrs_0.4.1 glue_1.6.2 stringi_1.7.6 proxy_0.4-26 compiler_4.1.0 [31] pillar_1.7.0 generics_0.1.2 classInt_0.4-3 pkgconfig_2.0.3 > sf::sf_extSoftVersion() GEOS GDAL proj.4 GDAL_with_GEOS USE_PROJ_H PROJ "3.9.1" "3.4.0" "8.1.1" "true" "true" "8.1.1"
报错核心原因:mapply/mcmapply默认开启SIMPLIFY=TRUE参数,会自动将传入的列表参数简化为矩阵/原子向量,直接破坏sf对象的S3类结构,把几何列拆成普通字符向量传入函数,最终触发st_drop_geometry的方法匹配错误,和输入数据是否包含字符列无关。
具体修复与优化方案如下:
- 调用
mapply/mcmapply时必须显式设置SIMPLIFY = FALSE,禁止函数自动简化输入输出结构,保留sf对象的完整类属性。 - 修正原函数的全局变量依赖问题:原函数直接调用全局环境的
blocks对象,并行计算时子进程无法稳定访问主进程全局变量,必须把当前县对应的普查区块子集作为显式参数传入函数。 - 优化数据拆分逻辑:提前将圆形缓冲区、点、普查区块三类数据统一按县拆分为长度一致的三个列表,每个子任务仅加载单县的小体量数据,大幅降低内存占用与跨进程数据传输开销。
- 优化函数内部运算效率:将循环内逐次
rbind累积绑定结果的逻辑改为列表存储、最终一次性绑定,同时给st_intersects加prepared = TRUE参数启用预编译几何索引,单任务运算速度可提升3~10倍。
修正后的函数与并行调用示例:
# 修复后的函数,移除全局变量依赖,优化运算效率 county_func_fixed <- function(x, y, county_blocks) { y_id <- st_drop_geometry(y) # 启用预编译几何索引加速空间匹配 x_int <- st_intersects(x, county_blocks, prepared = TRUE) # 用列表暂存循环结果,替代逐次rbind ints_list <- vector("list", length(x_int)) for(i in seq_along(x_int)){ # 直接用空间索引取子集,避免匹配行名的额外开销 blocks_int <- county_blocks[x_int[[i]], ] blocks_int$hud_id <- y_id[i, 1] ints_list[[i]] <- blocks_int } ints_holder <- do.call(rbind, ints_list) x_blocks <- st_intersection(x, ints_holder) %>% mutate(intersect_area = st_area(.)) %>% dplyr::select(GEOID10, intersect_area, hud_id) return(x_blocks) } # 并行调用示例 library(parallel) results <- mcmapply( FUN = county_func_fixed, x = county_circles_list, # 按县拆分的圆形缓冲区sf列表 y = county_points_list, # 按县拆分的对应点sf列表 county_blocks = county_blocks_list, # 按县拆分的普查区块sf列表 SIMPLIFY = FALSE, # 核心参数,禁止自动简化破坏sf结构 mc.cores = detectCores() - 1 # 预留1个核心避免系统卡顿 ) # 合并所有县的计算结果 final_result <- do.call(rbind, results)
内容的提问来源于stack exchange,提问作者mhankinson

