如何高效筛选与Data_2相交的Data_1中SpatialPolygonsDataFrame多边形

我拥有两个大型SpatialPolygonsDataFrame对象Data_1和Data_2，各自包含数百万个多边形。我希望筛选出Data_1中至少与Data_2中一个多边形相交的多边形。

我知道可以使用raster::intersect()查找两个多边形的交集，因此最初尝试用双重for循环逐一对比Data_1和Data_2中的多边形，保存相交多边形的ID，代码如下：

Find_intersection = function(Data_1){

Save_poli = NA

for (i in 1:length(Data_1@polygons)) {

  Poli_1 = Data_1@polygons[[i]]@Polygons[[1]]@coords
  Poli_1 = as.data.frame(Poli_1)

  coord <- data.frame(Longitude = Poli_1$x, Latitude = Poli_1$y)
  coordinates(coord) <- ~Longitude + Latitude
  proj4string(coord) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

  Poli_1 <- Polygon(coord)
  Poli_1 <- Polygons(list(Poli_1), 1)
  Poli_1 <- SpatialPolygons(list(Poli_1))
  proj4string(Poli_1) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

    for (j in 1:length(Data_2@polygons)) {
  
      Poli_2 = Data_2@polygons[[j]]@Polygons[[1]]@coords
      Poli_2 = as.data.frame(Poli_2)
  
      coord <- data.frame(Longitude = Poli_2$x, Latitude = Poli_2$y)
      coordinates(coord) <- ~Longitude + Latitude
      proj4string(coord) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84+towgs84=0,0,0")
  
      Poli_2 <- Polygon(coord)
      Poli_2 <- Polygons(list(Poli_2), 1)
      Poli_2 <- SpatialPolygons(list(Poli_2))
      proj4string(Poli_2) = CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
  
      Intersection = raster::intersect(Poli_1, Poli_2)
  
      if(is.null(Intersection) == F) {
         Save_poli = c(Save_poli, Data_1@polygons[[i]]@ID)
         break
       }
      }
    }
  return(Save_poli)
 }

Save_poli = Find_intersection(Data_1)

但该方法效率极低，因Data_1和Data_2的数据量过大，运行速度过慢。我也查找过相关问题，但解决方案类似或不适用。

请问是否存在无需逐一检查所有组合，即可找出Data_1中与Data_2多边形相交的多边形的方法？

直接对整个SpatialPolygonsDataFrame对象使用空间操作函数，而非逐一遍历，是处理大规模空间数据的核心优化思路，以下是两种高效实现方式：

方法1：使用sf包（推荐，空间索引优化） #

sf包默认支持空间索引，能大幅减少不必要的相交计算，适合百万级多边形场景：

1. 先将SpatialPolygonsDataFrame转换为sf对象：

library(sf)
sf_data1 <- st_as_sf(Data_1)
sf_data2 <- st_as_sf(Data_2)

2. 利用st_intersects()获取相交关系，返回的是每个sf_data1要素对应的sf_data2要素索引，过滤出至少有一个相交的要素：

# 获取相交关系矩阵
intersect_idx <- st_intersects(sf_data1, sf_data2, sparse = TRUE)
# 筛选出至少有一个相交的Data_1要素
filtered_data1 <- sf_data1[lengths(intersect_idx) > 0, ]
# 如果需要转回SpatialPolygonsDataFrame格式
filtered_sp_data1 <- as(filtered_data1, "Spatial")

st_intersects()内部会自动使用空间索引，只检查边界可能重叠的多边形组合，避免了O(n*m)的暴力遍历。

方法2：使用raster或sp包的批量操作 #

如果你倾向于继续使用sp体系，可直接对整个对象调用raster::intersect()，内部已做优化：

library(raster)
# 直接计算两个对象的交集，结果包含所有相交的Data_1多边形（带属性）
intersected_result <- intersect(Data_1, Data_2)
# 提取原始Data_1中对应的ID
save_poli <- unique(intersected_result@data$ID) # 假设ID是Data_1中的标识列
# 筛选出Data_1中符合条件的多边形
filtered_data1 <- Data_1[Data_1@data$ID %in% save_poli, ]

这种方式避免了手动拆解多边形的冗余操作，利用包内部的优化逻辑处理大规模数据。

额外优化建议 #

• 确保两个对象的投影坐标系完全一致，避免计算过程中自动投影带来的性能损耗；
• 若数据量极大，可考虑先对Data_2做空间合并（比如用unionSpatialPolygons()），减少需要对比的多边形数量；
• 优先选择sf包，其空间运算效率远高于传统sp包，且API更简洁。

内容的提问来源于stack exchange，提问作者Daniel B.G