Open Nowosad opened 1 month ago
Won't it be an issue related to GEOS? I get the same results in terra
and geos
, but strangely the first example in geos
looks different (pink color is inside, not outside).
library("geos")
p = as_geos_geometry("POLYGON ((50 50, 50 150, 150 150, 150 50, 50 50))")
p_b = geos_buffer(p, -20, params = geos_buffer_params(single_sided = TRUE))
plot(p_b, col = "lightpink", border = "maroon", lwd = 6)
plot(p, border = "deepskyblue3", add = TRUE, lwd = 6)
singleSide: logical; if ‘TRUE’, single-sided buffers are returned for
linear geometries, in which case negative ‘dist’ values give
buffers on the right-hand side, positive on the left; see
details
So it should'nt do anything for polygons, really ;-) but if you take the polygon as a linestring you get a buffer on the right hand side, right?
p = st_as_sfc("POLYGON ((50 50, 50 150, 150 150, 150 50, 50 50))")
p_b = st_buffer(st_cast(p, "LINESTRING"), -20, singleSide = TRUE)
plot(p_b, col = "lightpink", border = "maroon", lwd = 6)
plot(p, border = "deepskyblue3", add = TRUE, lwd = 6)
@edzer see this:
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.7.3, PROJ 9.2.1; sf_use_s2() is TRUE
nc = st_read(system.file("shape/nc.shp", package = "sf"))
nc_g = st_geometry(nc)
nc_g = st_transform(nc_g, "EPSG:9311")
nc_g1 = nc_g[1, ] |> st_cast("POLYGON") |> st_cast("LINESTRING")
nc_g1_b = st_buffer(nc_g1, -2000, singleSide = TRUE)
plot(nc_g1, col = "deepskyblue3", lwd = 6)
plot(nc_g1_b, col = "lightpink", add = TRUE, border = "maroon", lwd = 6)
It seems that the algorithm fails if the buffer self-intersects; this works:
library(sf)
# Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
nc = st_read(system.file("shape/nc.shp", package = "sf"))
# Reading layer `nc' from data source
# `/home/edzer/R/x86_64-pc-linux-gnu-library/4.4/sf/shape/nc.shp'
# using driver `ESRI Shapefile'
# Simple feature collection with 100 features and 14 fields
# Geometry type: MULTIPOLYGON
# Dimension: XY
# Bounding box: xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
# Geodetic CRS: NAD27
nc_g = st_geometry(nc)
nc_g = st_transform(nc_g, "EPSG:9311")
nc_g1 = nc_g[1, ] |> st_cast("POLYGON") |> st_cast("LINESTRING")
x = nc_g1
x[[1]] = x[[1]][-(25:27),] # remove last 3 vertices
nc_g1[[1]] = structure(x[[1]], class = class(nc_g1[[1]]))
nc_g1_b = st_buffer(nc_g1, -2000, singleSide = TRUE)
plot(nc_g1, col = "deepskyblue3", lwd = 6)
plot(nc_g1_b, col = "lightpink", add = TRUE, border = "maroon", lwd = 6)
Definitely an issue of GEOS, not sf
.
I was working on preparing an issue for GEOS, and then I found out that the geos package seems to return expected result.
@edzer what could make a difference here? Do you think this is still a GEOS problem?
library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.8.5, PROJ 9.3.1; sf_use_s2() is TRUE
nc = read_sf(system.file("shape/nc.shp", package = "sf"))
nc_g = st_geometry(nc)
nc_g = st_transform(nc_g, "EPSG:9311")
nc_g1 = nc_g[1, ] |> st_cast("POLYGON") |> st_cast("LINESTRING")
# sf
nc_g1_b = st_buffer(nc_g1, -2000, singleSide = TRUE)
plot(nc_g1, col = "deepskyblue3", lwd = 6)
plot(nc_g1_b, col = "lightpink", add = TRUE, border = "maroon", lwd = 6)
# geos
library(geos)
p = as_geos_geometry(nc_g1)
p_b = geos_buffer(p, -2000, params = geos_buffer_params(single_sided = TRUE))
plot(p_b, col = "lightpink", border = "maroon", lwd = 6)
plot(p, col = "deepskyblue3", add = TRUE, lwd = 6)
A possible cause outside the package code may be that geos
uses another version of GEOS
(3.11.1) than sf
(3.12.1). The 3.12.0 release notes mention work on single sided buffers.
Hi @edzer -- the
st_buffer()
withsingleSide = TRUE
works as expected on an example from https://postgis.net/docs/ST_Buffer.html, but differently on actual data. Is this an expected behavior?