[GIS] conditional st_join and st_intersection in R (sf)

postgisrsf

I just started working with the R sf package and am trying to rewrite a simple procedure from PostGIS to extract shared borders between polygons. In PostGIS, I usually do a "conditional" self-join that filters the rows before doing the intersection (e.g. to avoid intersecting a lineshape with itself): 

 WITH table1 AS 
 (SELECT 
    id,
    ST_Boundary(polygon) as geom 
  FROM table0)
 SELECT 
   t1.id, 
   t2.id, 
   ST_Intersection(t1.geom, t2.geom) as geom
   FROM table1 t1  
     LEFT JOIN table1 t2 
     ON ST_Intersects(t1.geom, t2.geom) 
     AND t1.id<t2.id

In R, I can get the same result by doing the intersection first and subsetting the results after:

library(sf)
library(raster)

## Load example data (raster package)
fr0 <- getData('GADM', country="FRA", level=0)
fr1 <- getData('GADM', country="FRA", level=1)

## Get border intersections between units
fr1 <- st_boundary(st_as_sf(fr1))
fr1.is <- st_intersection(fr1, fr1)

## Remove self-intersections
fr1.is <- fr1.is[fr1.is$ID_1<fr1.is$ID_1.1,]

Result:

Result

This works fine, but it would be much more efficient to do a conditional join first and then do the intersection, as is usually done in PostGIS. Unfortunately, my R skills are not that great yet. Any suggestions on how I can do this?

Best Answer

I can only suggest you loop over each feature and compute the intersection with all the features not considered so far:

nr = nrow(fr1)
is1 <- do.call(rbind,
       lapply(1:(nr-1),
          function(i){
           st_intersection(fr1[i,], fr1[-(1:i),])
           }))

> dim(is1)
[1] 43 27

which is the same as your code.

Whether this method is faster or not I don't know - would need to benchmark it under a variety of conditions.

Related Question