Efficient clipping of polygons by other polygons using R

Question

I have a dataset with overlapping polygons. I would like the clip the overlapping polygons (which I need to retain as overlapping) to a boundary (another polygon). I am using sf and st_intersection but due to the size of my data (60,000 polygons) this is not a quick process (I stopped it after an hour). For comparison I've tried using QGIS (I also stopped that after an hour with it still processing) and ArcMap (completed the task in 3 minutes). Is there a way I can complete the task using R in a reasonable time frame?

The code below recreates my data (I am after Points_Buffer being clipped to Boundary not Boundary_Clip):

library(sf)

download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/England_gor_2011_clipped.zip", 
              destfile = "England_gor_2011_clipped.zip" , mode='wb')
unzip("England_gor_2011_clipped.zip", exdir = ".")
file.remove("England_gor_2011_clipped.zip")

Boundary <- st_read("England_gor_2011_clipped.shp") %>% st_transform(., crs=27700) 

Boundary_Crop <- st_crop(Boundary,st_bbox(c(xmin = 401000, xmax = 500000, ymax = 172000, ymin = 74000), crs = st_crs(27700)))
                      
Points_Buffer <- st_sample(Boundary_Crop, 60000, exact = TRUE) %>% st_buffer(., 2000)

Points_Buffer_Clip <- st_intersection(Points_Buffer, Boundary)

Answer 1

So I have a solution that uses st_intersects to identify the polygons that intersect the boundary. I then run st_intersection on these polygons only before joining that output with the polygons that didn't need to be clipped.

library(sf)

# Original data

download.file("https://borders.ukdataservice.ac.uk/ukborders/easy_download/prebuilt/shape/England_gor_2011_clipped.zip", 
              destfile = "England_gor_2011_clipped.zip" , mode='wb')
unzip("England_gor_2011_clipped.zip", exdir = ".")
file.remove("England_gor_2011_clipped.zip")

Boundary <- st_read("England_gor_2011_clipped.shp") %>% st_transform(., crs=27700) 

Boundary_Crop <- st_crop(Boundary,st_bbox(c(xmin = 401000, xmax = 500000, ymax = 172000, ymin = 74000), crs = st_crs(27700)))
                      
Points_Buffer <- st_sample(Boundary_Crop, 60000, exact = TRUE) %>% st_buffer(., 2000)

# Solution

Boundary_Line <- st_cast(Boundary,"MULTILINESTRING") %>% st_cast(.,"LINESTRING") 

Buffer_Line_Intersect <- st_intersects(Boundary_Line,Points_Buffer) %>% unlist()

Points_Buffer_Does_Intersect <- Points_Buffer %>% .[Buffer_Line_Intersect ,]

Points_Buffer_Does_Not_Intersect <- Points_Buffer %>% .[-Buffer_Line_Intersect ,]

Boundary_Buffer_Crop <- st_crop(Boundary,st_bbox(Points_Buffer), crs = st_crs(27700))

Points_Buffer_Clip <- st_intersection(Points_Buffer_Does_Intersect, Boundary_Buffer_Crop)

# If your input polygons contain fields you may want to 'dissolve' the multiple single polygons the above process may produce to create multipart polygons using something along the lines of:
# Points_Buffer_Clip <- Points_Buffer_Clip %>% group_by(n) %>% summarise(n) %>% st_cast("MULTIPOLYGON") 

# If you have polygons with fields then you will have to rbind instead of this
Points_Buffer_Complete <- c(Points_Buffer_Does_Not_Intersect, Points_Buffer_Clip)

This whole process takes about 10 minutes for me. So it's longer than ArcMap, but not too long for my purposes.