R spatial: Erase one polygon from another: correct use of `st_difference`?

Question

I am trying to erase one shapefile from the other, similar to ArcGIS Erase tool:

This worked pretty well with rgeos::gDifference(spdf1, spdf2), which needs as input geometry of SpatialPolygonDataFrame.

Now, I want to switch from SPDF files to sf or sfc objects. But, st_difference does not keep a not overlapping geometry. Instead, it keeps buffer minus one :

Instead, I want to keep part encircled in blue:

How to correctly Erase features using sf objects? I tried tools st_difference, st_sym_difference or using features st_geometry(sf) for operation, but I get to have the same, for me unwanted results. Also, I found a suggestion of st_erase function:

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

But using it on my data erased <- st_erase(buff, u) generates error:

`Error in CPL_geos_op2(op, st_geometry(x), st_geometry(y))   Evaluation error: TopologyException: Input geom 1 is invalid: Self-intersection at or near point`

Here is my example:

# Load data
shp = system.file("shape/nc.shp", package="sf")

my.sf <- st_read(shp, quiet = TRUE)

# Convert crs to projected system to make buffer
my.sf.web<- st_transform(my.sf, 3857)

# Subset the data to create two independent shps
i = 10

# Split datasets in two files
one  = my.sf.web[i, ]
left = my.sf.web[-i,]

# Create buffer 
buff = st_buffer(one, 35000 ) # distance


# CHeck which polygons overlaps with my buffer
out.overlap = st_overlaps(buff, left)


# Subset the polygons that overlaps with the buffer
nbrs.buff <- left[st_overlaps(buff,left)[[1]],]


u <- st_union(st_geometry(nbrs.buff), st_geometry(one), by_feature = FALSE)
#u <- st_union(st_combine(st_geometry(nbrs.buff)), 
 #             st_combine(st_geometry(one)))
int.buff.one = st_difference(buff, u)   # NOT WORKING HERE???

Answer 1

Your code worked just fine for me, so I think that make it throws that error because you have some topology error in your shape. I've just adapted your code adding a check and repair part in order to avoid error (st_is_valid and st_make_valid).

Here you can download the data and code I've used.

library(sf)
library(mapview)

# Load data
my.sf <- st_read("shape/shape.gpkg", quiet = TRUE)

# Check for errors
st_is_valid(my.sf, reason = TRUE)

# Repare if needed
my.sf <- st_make_valid(my.sf)

# Subset the data to create two independent shps
i = 2
one  = my.sf[i, ]
left = my.sf[-i,]

# Create buffer 
buff = st_buffer(one, 1000) # distance

# CHeck which polygons overlaps with my buffer
out.overlap = st_overlaps(buff, left)

# Subset the polygons that overlaps with the buffer
nbrs.buff <- left[st_overlaps(buff,left)[[1]],]

u <- st_union(st_geometry(nbrs.buff), st_geometry(one), by_feature = FALSE)

int.buff.one = st_difference(buff, u) 

# RESULTS
mapview(my.sf)+mapview(int.buff.one, alpha=0.5, color="red",lwd=3)