Function in R to get difference between two sets of polygons comparable in speed to QGIS

Question

I have two very large polygon shapefiles for which I need to calculate the difference between the two shapefiles. One is a shapefile of forest polygons, the other is a shapefile of road polygons. I want to remove the area of the forest polygons that intersects with road polygons, while maintaining each polygon as a separate feature and its data.

After much trial and error, I found that erase using the sp format or using st_erase using the sf format performs the operation that I need.

st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))

However, the processing time is really long. With a small sample of my dataset, I have processing times several minutes long. (And R crashes with a tenth of the dataset). Meanwhile, the function Vector-->Geoprocessing tools--> Difference in QGIS takes 83 seconds with a tenth of the dataset.

erase(forest_sp, roads_sp)
#2.83488 minutes 
forest_sp-roads_sp
#2.8698 minutes
st_erase(forest, roads) 
#58.77sec

Is there a function in R performs this operation quickly? Or is R inferior to QGIS in processing times for this kind of spatial function?

Background

There are multiple questions on gis.stackexchange that deal with this type of problem. Some are with the sp format and suggest gDifference or erase (or its analog polygon1-polygon2):

Erase, gDifference

Clip polygon and retain data?

gDifference

Remove parts of line that fall within a polygon

Identify area within one shapefile not in another shapefile in R

Reverse clipping (erasing ) in R?

Other questions suggest solutions with the sf format: st_difference

Reverse clipping (erasing ) in R?

Using simple features (sf) in R, how do I erase polygons overlapping with another layer

ogr2ogr equivalent of QGIS Union

All state that they perform that same thing (remove intersecting portions) and the official tag definition on gis.stackexchange states that a difference operation removes sections. However, I have found that often they don't. For example, gDifference with byID=FALSE that merges all into one MULTIPOLYGON does remove overlapping sections but gDifference with byID=TRUE that keeps polygons separate includes overlapping sections. With sf, st_difference also keeping overlapping sections and ends up with 10x as many polygons as my original file.

However, erase and the helper function st_erase do work as expected.

Difference

A "difference" operation is a geoprocessing operation that completely removes one overlapping polygon and the intersection of that polygon from the other. It uses one polygon to "take a bite" out of another.

Erase

The action of removing data attributes or features through the processes of positive or negative selection and deletion.

Answer 1

The rmapshaper package offers a faster function for "erasing" or finding the difference between two polygon layers in R. Below is a comparison of the st_erase approach with sf and the ms_erase function in rmapshaper.

# Let's try an example using publicly available Census vector layers through `tigris`
if(!require(dplyr)){install.packages("dplyr")}
library(dplyr)
if(!require(sf)){install.packages("sf")}
library(sf)
if(!require(tigris)){install.packages("tigris")}
library(tigris)
if(!require(rmapshaper)){install.packages("rmapshaper")}
library(rmapshaper)

# First download two separate layers that overlap. We'll use a layer of Census 
# tracts and a layer of areal water features for Worcester County in Massachusetts.
# worcester_tracts has 172 obs and  10 variables; about 285KB in size.
worcester_tracts <- tracts(state = "MA", county = "Worcester", cb = TRUE) %>% 
  st_transform(., crs = 2805) %>% # use projected layers when doing geometric stuff
  st_make_valid()
# worcester_water has 3,270 obs and 9 variables; about 4.9MB in size. 
worcester_water <- area_water(state = "MA", county = "Worcester") %>% 
  st_transform(., crs = 2805) %>% 
  st_make_valid()

# Start with st_erase function. Note that I altered function slightly by removing 
# the st_combine function because it throws an error and isn't necessary. 
# Takes about 27 secs.
st_erase = function(x, y) st_difference(x, st_union(y))
system.time(erased_tracts1 <- st_erase(worcester_tracts, worcester_water))

# Now try it with ms_erase. Takes about 4.3 secs; over 6x faster. 
system.time(erased_tracts2 <- ms_erase(worcester_tracts, worcester_water))

The results are identical. You can speed things up even more by retaining as few variables, or attributes, as possible in the layers when doing these operations. Also, if your computer memory is being overwhelmed, consider breaking up the process into smaller geographic chunks (e.g., municipalities, counties, a grid), rather than working on the whole dataset at one time.

Answer 2

It looks like something has been updated in sf or rmapshaper as sf runs faster on my system now (see benchmark below).

Also if you only do the intersection on features that st_intersect() this can speed up the operation a little too. It can make a big difference if your x is much smaller than your y or vice versa.

### Install/load packages
# Install packages if you dont have them
if(!require(dplyr)){install.packages("dplyr")}
if(!require(sf)){install.packages("sf")}
if(!require(tigris)){install.packages("tigris")}
if(!require(rmapshaper)){install.packages("rmapshaper")}
if(!require(rbenchmark)){install.packages("rbenchmark")}

# load packages
library(dplyr)
library(sf)
library(tigris)
library(rmapshaper)
library(rbenchmark)


### Load Example layers
worcester_tracts <- tracts(state = "MA", county = "Worcester", cb = TRUE) %>% 
  st_transform(., crs = 2805) %>% # use projected layers when doing geometric stuff
  st_make_valid()
worcester_water <- area_water(state = "MA", county = "Worcester") %>% 
  st_transform(., crs = 2805) %>% 
  st_make_valid()


### User created functions to benchmark
# Getting difference of geoms of x and y that partially intersect
st_erase_big = function(x, y){
  x$ref <- 1:nrow(x)
  x_inter = x[apply(st_intersects(x,y), 1, any),] #x in y
  y_inter = y[apply(st_intersects(y,x), 1, any),] #y in x
  difference <- st_difference(x_inter, st_union(y_inter)) #x_inter not in y_inter
  nointer_x = x[!x$ref %in% x_inter$ref,] # x not intersecting
  dplyr::bind_rows(difference,nointer_x)
} 


st_erase = function(x, y) st_difference(x, st_union(y))


### Benchmarking
benchmark("st_erase" = {st_erase(worcester_tracts, worcester_water)},
          "ms_erase" = {ms_erase(worcester_tracts, worcester_water)},
          "st_erase_big" = {st_erase_big(worcester_tracts, worcester_water)},
          replications = 5,
          columns = c("function", "replications", "elapsed","relative"))

The fastest is st_erase_big by a whisker, closely followed by st_erase and ~50% slower is ms_erase

 function     replications elapsed relative
 ms_erase            5   30.31    1.452
 st_erase            5   21.05    1.008
 st_erase_big        5   20.88    1.000