The st_buffer() with singleSide = TRUE is inconsistent

[Nowosad]

Hi @edzer -- the st_buffer() with singleSide = 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?

library(sf)
#> Linking to GEOS 3.12.1, GDAL 3.7.3, PROJ 9.2.1; sf_use_s2() is TRUE
# works as expected
p = st_as_sfc("POLYGON ((50 50, 50 150, 150 150, 150 50, 50 50))")
p_b = st_buffer(p, -20, singleSide = TRUE)
plot(p_b, col = "lightpink", border = "maroon", lwd = 6)
plot(p, border = "deepskyblue3", add = TRUE, lwd = 6)

# ?
nc = st_read(system.file("shape/nc.shp", package = "sf"))
#> Reading layer `nc' from data source 
#>   `/home/jn/R/x86_64-redhat-linux-gnu-library/4.3/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")
nc_g1_b = st_buffer(nc_g1, -2000, singleSide = TRUE)

plot(nc_g1, border = "deepskyblue3", lwd = 6)
plot(nc_g1_b, col = "lightpink", add = TRUE, border = "maroon", lwd = 6)

[kadyb]

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)

[edzer]

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)

[Nowosad]

@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)

[edzer]

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.

[Nowosad]

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)

[edzer]

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.