conditional st_join and st_intersection in R (sf)

Question

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:

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?

Answer 1

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.

Answer 2

You can split the difference and use your postGIS query but call it from R using sf::st_read. Within st_read you can add a connection to your postGIS db and run the query from there. This gives you the flexibility of using your preferred SQL without having to leave R. Here's a sample of what I've done for a project grabbing LEHD data for a city.

if(!require(pacman)){install.packages("pacman"); library(pacman)}
p_load(here, RPostgreSQL, sf)


user <- "jamgreen"
host <- "pgsql102.rc.pdx.edu"
pw <- scan(here::here("batteries.pgpss"), what = "")
dbname <- "bike_lanes"

con <- dbConnect("PostgreSQL", host = host, user = user, dbname = dbname, 
             password = pw, port = 5433)

minn_corridor <- st_read(dsn = con, query = "select a.geoid10 as geoid, a.c000 as     c000,
a.cns07 as cns07, a.cns08 as cns08, a.cns12 as cns12, a.cns14 as cns14, 
a.cns15     as cns15, a.cns16 as cns16, a.cns17 as cns17, 
a.cns18 as cns18, a.cns19 as cns19,     b.name as Name, 
b.buildstart as buildstart, b.buildend as buildend, b.group as corridor_group, 
b.type as grouptype,  a.year as year, a.geometry as geom
FROM minneapolis_lehd a, minneapolis_corridors b
WHERE ST_Intersects(ST_Buffer(b.geom, 20), a.geometry);")

This is a relatively simple call, but st_read will run a SQL query of arbitrary length as long as it is a single call is my understanding.