R intersect two list of multiple polygons-CodePudding

I have two objects Data_1 and Data_2 in the form of two Large SpatialPolygonsDataFrame each containing several million polygons.

I'm trying to filter Data_1 in order to leave only those polygons which intersect with at least one polygon in Data_2.

I know that I can find the intersection of two polygons using something like raster::intersect(), so my first brute force attempt was to make use of two for cycles to compare each polygon in Data_1 with each one in Data_2 and save the ID of those intersecting:

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)

But this approach is inefficient and the size of both Data_1 and Data_2 makes it too slow. I also found related questions (https://gis.stackexchange.com/questions/156660/loop-to-check-multiple-polygons-overlap-in-r) but the solutions are similar or don't apply.

Is there a way to know which polygons in Data_1 intersect with the polygons in Data_2 without checking every single one of the combinations?

Thanks!

CodePudding user response：

You can use terra::is.related for that like this

library(terra)
d1 <- vect(Data_1)
d2 <- vect(Data_2)

i <- is.related(d1, d2, "intersects") 
x <- d1[i,]

Illustrated with some example data:

library(terra)
p1 <- vect("POLYGON ((0 0, 8 0, 8 9, 0 9, 0 0))")
p2 <- vect("POLYGON ((5 6, 15 6, 15 15, 5 15, 5 6))")
p3 <- vect("POLYGON ((8 2, 9 2, 9 3, 8 3, 8 2))")
p4 <- vect("POLYGON ((2 6, 3 6, 3 8, 2 8, 2 6))")
p5 <- vect("POLYGON ((2 12, 3 12, 3 13, 2 13, 2 12))")
p6 <- vect("POLYGON ((10 4, 12 4, 12 7, 11 7, 11 6, 10 6, 10 4))")

d1 <- rbind(p3, p4, p5, p6)
d2 <- rbind(p1, p2)

i <- is.related(d1, d2, "intersects")
x <- d1[i,]

plot(d2, border="gray")
lines(d1, col="blue", lwd=4)
lines(x, col="red", lwd=2)