CODEWITHSUNDEEP
rpolygondistanceshapefilegeoplot
MapQuestioner
MapQuestioner

Reputation: 1

Incorrect Result From st_difference

I am trying to generate a set of rings around a set of points in the sf package, and have run into an issue where for some reason st_difference returns the whole first geometry despite the two inputs overlapping. I want a ring of all points between 19 and 20 units away from the center point.

Here is my sample code

library(sf)
library(sp) # SpatialPoints
set.seed(42)

x_coords <<- runif(n = 10, min = -100, max = 100)
y_coords <<- runif(n = 10, min = -100, max = 100)

PhenomPoints <<- SpatialPoints(coords = cbind(x_coords,y_coords)) 
PhenomPoints_SF <<- st_as_sf(PhenomPoints) 
buffer_zones <<- st_buffer(PhenomPoints_SF, dist=20)
buffer_zones_2 <<- st_buffer(PhenomPoints_SF, dist=19)
ring <- st_difference(buffer_zones,buffer_zones_2)

plot(buffer_zones, col = "Grey")
plot(buffer_zones_2, col = "Blue", add=TRUE)
plot(ring, col = "Yellow", add=TRUE)

However, the ring seems to be basically identical to the buffer_zones. I tried using st_intersects, and it generated an intersection, but when I then used that with st_differences, it also returned the original. To see if the issue was the blank CRS, I tried it with a random CRS value and it still didn't work.

Upvotes: 0

Views: 126

Answers (1)

margusl
margusl

Reputation: 17504

For pairwise difference you probably need to iterate over geometry pairs yourself, one option for this is calling st_difference() through mapply().

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE

set.seed(42)
x_coords <- runif(n = 10, min = -100, max = 100)
y_coords <- runif(n = 10, min = -100, max = 100)

points_sfc <-
  st_multipoint(cbind(x_coords, y_coords)) |>
  st_sfc() |>
  st_cast("POINT")

buffer_zones   <- st_buffer(points_sfc, dist=20)
buffer_zones_2 <- st_buffer(points_sfc, dist=15)

ring <- 
  mapply(st_difference, buffer_zones, buffer_zones_2, SIMPLIFY = FALSE) |>
  st_sfc()

plot(points_sfc, xlim = c(-120, 120), ylim = c(-120, 120))
plot(ring, col = "Yellow", add= TRUE)

Make sure your arguments are sfc / geometrty columns. If your buffers are sf objects, first extract geometries:

mapply(st_difference, st_geometry(buffer_zones), st_geometry(buffer_zones_2), SIMPLIFY = FALSE)

Upvotes: 1

Related Questions