st_buffer ignores nQuadSegs when sf_use_s2(TRUE)

[juansebastianl]

I'm not sure if this is an issue or expected behavior (sort of like the st_crop issue) but it took me a long time to figure out what was going on, so even if it is expected behavior I think it should be documented.

If I have the default sf_use_s2(TRUE) then st_buffer behaves completely differently than expected. For example running

tibble(id = 1, lon = 0, lat = 0) %>% 
  st_as_sf(coords = c("lon", "lat")) %>% 
  st_set_crs(4326) %>% 
  st_buffer(111000, nQuadSegs = 1, endCapStyle = "SQUARE") %>% 
  ggplot() + 
     geom_sf()

produces a circle of radius 1 degree no matter the value of nQuadSegs or endCapStyle, If I run this after running sf_use_s2(FALSE) I would get a square of radius larger than the whole earth.

[edzer]

Yes, agreed this needs improved documentation; see e.g. #1801 for a successful attempt to clear up a similar problem for st_simplify().

[juansebastianl]

Okay, it seems nQuadSegs is completely ignored because the s2 library doesn't reproduce the same ability to buffer as geos, right? Specifically this issue https://github.com/r-spatial/s2/issues/42 . There is code in the st_buffer.sfc function which is commented out which outputs a warning. Uncommenting that would be helpful, but it doesn't seem to make sense to me why this function should use s2 as a default when s2 can't actually produce the same results.

Would it be possible to have it use geos and have a warning and a flag along the lines of warning: st_buffer uses geos instead of s2 for calculation by default, use s2 = T to use s2 and then a warning when using s2 along the lines of warning: st_buffer ignores all arguments when using s2 = T .

When the issue in s2 is resolved then it can be reverted to the current default behavior.

[edzer]

but it doesn't seem to make sense to me why this function should use s2 as a default when s2 can't actually produce the same results.

because GEOS can only create round buffers on the equator, everywhere else they are ellipses, getting worse the further you go away from the equator. s2 creates round shapes on the sphere, regrettably not with smooth borders.

When the issue in s2 is resolved then it can be reverted to the current default behavior.

I haven't received any signals that this will happen any time soon. You can now switch off using s2geometry by sf_use_s2(FALSE). I see the imperfection but let me put it this way: how would you create a map with the ocean zones not further than 200 mile from the nearest coast? with GEOS?

[juansebastianl]

Yeah, it's definitely a trade-off. It's just frustrating from an user's perspective because even without the extra options this is still a breaking change. Before s2 the buffer size was in degree cells, now by default it's in meters. I can see why for new code/new users this is definitely better behavior, I just wish there was a way to make it a non-breaking change.

[juansebastianl]

I still think the meters vs degree change isn't documented but I couldn't think of a reasonable way to show a warning or add it to the docs.

[edzer]

Units can still be in degrees:

library(sf)
library(tidyverse)
tibble(id = 1, lon = 0, lat = 50) %>% 
  st_as_sf(coords = c("lon", "lat")) %>% 
  st_set_crs(4326) %>% 
  #st_buffer(111000, nQuadSegs = 1, endCapStyle = "SQUARE") %>% 
  st_buffer(units::set_units(1, degree)) %>%
  ggplot() + 
     geom_sf()

[juansebastianl]

I agree, it's just a breaking change. I can't think of how to make it non-breaking, though, so i think my PR is all that can be done. The only other thing I can think of it so add an output "buffering with units: " to the function.

[bart1]

Just for information I just encountered the same issue while i was trying to create simpler polygons:

require(sf)
#> Loading required package: sf
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1
require(units)
#> Loading required package: units
#> udunits database from /usr/share/xml/udunits/udunits2.xml
pt<-st_sfc(st_point(1:2), crs=st_crs(4326))
a<-st_buffer(pt, set_units(1,'m'), nQuadSegs = 30)
b<-st_buffer(pt, set_units(1,'m'), nQuadSegs = 3)
all.equal(a,b)
#> [1] TRUE

^{Created on 2021-09-29 by the reprex package (v2.0.1)}

[edzer]

so i think my PR is all that can be done.

I agree, but you didn't finish it.

[bart1]

I think that is good enough. For future reference as in my case the shape was not very critical st_simplify(..., dTolerance=...) managed the desired effect of having a simpler polygon