Geographic Information Systems Asked on May 15, 2021
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???
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)
Answered by César Arquero on May 15, 2021
Get help from others!
Recent Questions
Recent Answers
© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP