2

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.

3
  • 2
    What version of sf are you using, and what CRS is used in your objects? I believe QGIS uses GEOS, while sf (from version 1.0 above, and for unprojected coordinate systems, ie. 4326 and the like) uses s2 from Google. This may, and may not, be relevant... Commented Nov 24, 2021 at 15:07
  • 1
    @JindraLacko The latest stable version of sf and the CRS is EPSG:4326 - WGS 84
    – entropy
    Commented Nov 24, 2021 at 15:49
  • @JindraLacko What you said is indeed very relevant for this case. GEOS is much faster and I recommend using st_transfrom() to reproject your layers to a suitable projection and then switch-off s2 (sf_use_s2) before running any sf function in the family: st_intersection, st_difference, etc. Commented Nov 29, 2023 at 8:27

2 Answers 2

4

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()

borders of county mecklenburg

1
  • Thanks for the response, but I can't even test this with my data as when i use st_intersection it seems to just run forever
    – entropy
    Commented Nov 24, 2021 at 17:36
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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.