基于坐标加速从Shapefile提取生物地理区域的方法

问题背景 #

此前基于sf包逐行匹配坐标与生物地理区域的方法，在处理大规模坐标集时速度过慢，尝试用Terra包转栅格提取未成功，现提供三种高效替代方案。

方案1：使用sf的st_join向量化匹配 #

直接用空间连接替代逐行循环，sf会自动优化空间查询，速度提升显著，同时自动处理不在区域内的坐标（返回NA）。

代码实现 #

# 加载包
library(sf)

# 读取并转换生物地理区域Shapefile
bioregions <- read_sf('PATH/BiogeoRegions2016.shp')
bioregions <- st_transform(bioregions, crs = 3035) 

# 生成大规模点集（示例）
pt <- 
  data.frame(x = rep(c(14,28, 3, 4),1000), 
             y = rep(c(41, 59, 42, 42),1000), 
             id = 1:4000) |>
  st_as_sf(coords = c('x', 'y')) |>
  st_set_crs(4326) |>
  st_transform(3035)

# 空间连接匹配区域
result <- st_join(pt, bioregions[, c("short_name")], join = st_within, left = TRUE)
# 注：left=TRUE保留所有点，不在区域内的short_name为NA

方案2：使用terra包的矢量匹配 #

无需转栅格，直接用terra的矢量操作处理，速度比sf更快，适合超大规模数据。

代码实现 #

# 加载包
library(terra)

# 读取生物地理区域Shapefile并转换投影
bioregions <- vect('PATH/BiogeoRegions2016.shp')
bioregions <- project(bioregions, "EPSG:3035")

# 生成大规模点集并转换投影
pt_df <- data.frame(x = rep(c(14,28, 3, 4),1000), 
                    y = rep(c(41, 59, 42, 42),1000), 
                    id = 1:4000)
pt <- vect(pt_df, geom = c("x", "y"), crs = "EPSG:4326")
pt <- project(pt, "EPSG:3035")

# 提取匹配的区域信息
result <- extract(bioregions, pt, "short_name")
# 合并原数据与结果
result_df <- cbind(pt_df, result[, "short_name"])
# 不在区域内的short_name为NA

方案3：显式构建空间索引（sf进阶优化） #

如果数据量极端庞大，显式为区域图层构建空间索引，进一步加速查询。

代码实现 #

library(sf)

bioregions <- read_sf('PATH/BiogeoRegions2016.shp')
bioregions <- st_transform(bioregions, crs = 3035) 

# 显式构建空间索引
bioregions <- st_make_valid(bioregions) # 确保图层拓扑有效
st_geometry(bioregions) <- st_precision(st_geometry(bioregions), 1000) # 降低精度提升速度
st_crs(bioregions)$epsg <- 3035
bioregions <- st_set_geometry(bioregions, st_geometry(bioregions))

pt <- 
  data.frame(x = rep(c(14,28, 3, 4),1000), 
             y = rep(c(41, 59, 42, 42),1000), 
             id = 1:4000) |>
  st_as_sf(coords = c('x', 'y')) |>
  st_set_crs(4326) |>
  st_transform(3035)

result <- st_join(pt, bioregions[, c("short_name")], join = st_within, left = TRUE)

内容的提问来源于stack exchange，提问作者J. Lan