Workaround needed for sf/S2 bug for long/lat with st_buffer

[finnlindgren]

When sf::st_buffer() is given long/lat CRS coordinates and S2 is active (the default), negative dist argument seems to return an empty geometry, likely due to misinterpreting the units (which are supposed to be interpreted in metres for that case). The sf documentation just says that if an empty geometry is returned for negative dist, turn off S2. This is not a fully workable solution here, since that's up to the package user and affect the units interpretation.

Potential workaround: detect long/lat CRS and define a predictable behaviour of some kind, e.g. by extracting the raw coordinates, and working in raw degrees, or by locally turning off S2 to get predictable behaviour.

Potential workaround 2: see if st_concave_hull() is available, and work out an algorithm that takes advantage of it.

[enoch26]

Reprex example

library(sf)
#> Linking to GEOS 3.11.1, GDAL 3.6.2, PROJ 9.1.1; sf_use_s2() is TRUE
df <- data.frame(x = c(0, 1, 0, 1), y = c(0, 0, 1, 1))
df_sf <- st_as_sf(df, coords = c("x", "y"), crs = 4326) %>% st_bbox() %>% 
  st_as_sfc()
st_buffer(df_sf, dist = -0.25)
#> Geometry set for 1 feature  (with 1 geometry empty)
#> Geometry type: POLYGON
#> Dimension:     XY
#> Bounding box:  xmin: NA ymin: NA xmax: NA ymax: NA
#> Geodetic CRS:  WGS 84
#> POLYGON EMPTY
reprex::reprex()
#> ℹ Non-interactive session, setting `html_preview = FALSE`.
#> CLIPR_ALLOW has not been set, so clipr will not run interactively
#> Error in switch(where, expr = stringify_expression(x_expr), clipboard = ingest_clipboard(), : EXPR must be a length 1 vector

[finnlindgren]

Thanks! I got the s2 buffer to work, using the trick of inverting the polygon within an expanded bounding box. The resulting polygon from s2 is not very smooth though, so I wouldn't use it anyway, but with this method one would at least get something:

# Special nonconvex_hull when x is in longlat and s2 is active
special <- function(x, convex, concave) {
  x_expanded <- st_buffer(x, convex + concave)
  # Create covering
  x_box <- (st_bbox(x_expanded) + c(-0.1, -0.1, 0.1, 0.1)) %>% st_as_sfc()
  x_inverse <- st_difference(x_box, x_expanded)
  x_inverse_expanded <- st_buffer(x_inverse, concave)
  x_result <- st_difference(x_expanded, x_inverse_expanded)
  x_result
}

# longlat coordinates and S2, so units are in metres:
plot(special(df_sf, 50000, 50000), col = 4)
plot(df_sf, col = 2, add = TRUE)

# no crs, ordinary units:
df2 <- st_as_sf(df, coords = c("x", "y")) %>% st_bbox() %>% 
  st_as_sfc()
plot(special(df2, 0.5, 0.5), col = 4)
plot(df2, col = 2, add = TRUE)

# non-longlat coordinates, units are whatever the crs says (here km)
df_tr <- st_transform(df_sf, crs = fm_crs("mollweide_globe"))
plot(special(df_tr, 50, 50), col = 4, axes = TRUE)
plot(df_tr, col = 2, add = TRUE)

[finnlindgren]

With the new workaround, this works:

plot(fm_nonconvex_hull(df_sf, 50000), col = 2)