Why is there a speed difference between sf_intersection in R and QGIS intersection

Question

I have two shape files that each have around 4600 polygons. I want to find the overlapping area of these polygons.

In R, I am using the following code:

intersect_pct <- st_intersection(ward_2016, ward_2020) %>% 
  mutate(intersect_area = st_area(.)) 

Where ward_2020 and ward_2016 are my shape files. Unfortunately, I never even get to the point where I can call mutate(intersect_area = st_area(.)) as the st_intersection takes so long to compute.

However, in QGIS it takes about 30 seconds to compute the intersection.

I tried using sf_join. This computes quickly, however this treats two polygons that share the same boarder as intersecting.

I am using an M1 MacBook Pro.

Answer 1

The speed difference is likely due to QGIS using GEOS backend and the current stable sf (v1.0 plus) using s2 library from Google; for more information consider official documentation https://r-spatial.github.io/sf/articles/sf7.html

The s2 library / package introduces a new dimension to intersections - the concept of a model. In the default behaviour (model = "semi-open" for polygons) each polygon contains only half of its boundaries - as consequence bordering polygons don't have any points in common ({sf} sometimes overrides this default; it is new and not fully settled yet).

This goes against the logic of GEOS (where polygons contain their boundaries in the entirety = bordering polygons share a line).

Good piece of news is that you can use the model argument to drive the behaviour - model = "closed" forces the old behavior, model = "open" should in theory instruct sf::st_join() to ignore case when only border is shared.

For an example of impact of model consider this example, built on the well known & much loved nc.shp from {sf}.

In the first case (with defaults) County Mecklenburg shares no points with the rest of NC; in case of a closed model it shares the internal boundaties / the polygon is open, because to the south from Charlotte is South Carolina, out of scope of the nc.shp file.

library(sf)

shape <- st_read(system.file("shape/nc.shp", package="sf")) %>%
  st_transform(4326)

# cnty mecklenburg - as in Charlotte of Mecklenburg-Strelitz
mecklenburg <- shape[68, ]
rest_of_nc <- shape[-68, ]

# default - nothing
st_intersection(rest_of_nc,
                mecklenburg) %>% 
  st_geometry() %>% 
  plot()

# closed model - all internal boundaries of cnty Mecklenburg
st_intersection(rest_of_nc,
                mecklenburg,
                model = "closed") %>% 
  st_geometry() %>% 
  plot()

Answer 2

While the polygon model may matter here, my reading is this a question regarding speeding up the intersect. One alternative is to use st_intersects instead of st_intersection and then filter or subset on the logical outcomes of that call. Chapter 4 of Geocomputation With R has a good working example of this.

From that chapter:

    sel_sgbp = st_intersects(x = nz_height, y = canterbury)
class(sel_sgbp)
#> [1] "sgbp" "list"
sel_sgbp
#> Sparse geometry binary predicate list of length 101, where the
#> predicate was `intersects'
#> first 10 elements:
#>  1: (empty)
#>  2: (empty)
#>  3: (empty)
#>  4: (empty)
#>  5: 1
#>  6: 1
#>  7: 1
#>  8: 1
#>  9: 1
#>  10: 1
sel_logical = lengths(sel_sgbp) > 0
canterbury_height2 = nz_height[sel_logical, ]

st_intersects has some performance advantages here and may be more useful for a meatier query as you've described.