6

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)
2
  • Do you need to keep the attributes for the individual polygons in Points_Buffer or can you dissolve/merge all of the polygons before running the clip?
    – Aaron
    Aug 30 '20 at 1:46
  • @Aaron I need to keep the attributes for the individual polygons
    – Chris
    Aug 30 '20 at 22:36
0

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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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