8

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

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: enter image description here

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???
6
  • I don't understand what u is meant to be here. Each feature is the combination of the source feature and the feature you buffered. That is complicating the difference situation. plot(st_geometry(st_difference(buff, st_union(u)))) seems to extract the bit of the buffer above u but I'm not sure the way you construct u makes sense in the first place... – Spacedman Mar 11 '20 at 11:28
  • Hi @Spacedman, the u should just spatially merge my first geometry one and it's neighbors nbrs.buff and create single file. So, I just simply want to erase u from the buffer. Thanks, plot(st_geometry(st_difference(buff, st_union(u)))) makes correctly erase u from buff. Do you think I should rerun union tool u2 <- st_union(u) to have correct output? Or how would you simplify the union operation? – maycca Mar 11 '20 at 11:41
  • 1
    If you want to add a feature back to a spatial data frame, rbind: uu = rbind(nbrs.buff, one). – Spacedman Mar 11 '20 at 14:22
  • Have you checked if either u or buff have self-intersection problems? – Mox Mar 11 '20 at 17:07
  • @Mox, great idea!! I don't know how to check for self-intersection problems, but I am sure that they do, as the one vertices completely overlaps with vertices in nbrs.buff, as I one was to nbrs.buff geometry!!! when I create u = st_union(rbind(nbrs.buff, one)) (as suggested by @Spacedman) and then run st_difference() it creates correct output! I'll test st_dissolve and st_combine to see which one is more appropriate to run. – maycca Mar 12 '20 at 9:17
3

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)

enter image description here

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.