Unexpected behaviour of st_covers() with increasing size of x

I'm working with spatial data in the sf package in R, where I need to check if a set of randomly generated points falls within a bounding box. I'm using the st_covers() function to perform the coverage test. It works fine with the original bounding box and even with a slightly expanded bounding box, but behaves unexpectedly when the bounding box is expanded beyond a certain threshold.

library(sf)

set.seed(123)

# Define expansion function
expand_bbox <- function(bbox, expand_factor = 0.1) {
  # Get current bbox dimensions
  x_range <- bbox["xmax"] - bbox["xmin"]
  y_range <- bbox["ymax"] - bbox["ymin"]
  
  # Expand the limits, adjusting both directions correctly
  bbox["xmin"] <- max(bbox["xmin"] - (expand_factor * x_range), -180)
  bbox["xmax"] <- min(bbox["xmax"] + (expand_factor * x_range), 180)
  bbox["ymin"] <- max(bbox["ymin"] - (expand_factor * y_range), -90)
  bbox["ymax"] <- min(bbox["ymax"] + (expand_factor * y_range), 90)
  
  return(bbox)
}

# Create random longitude and latitude points within bbox of South America
points_sf <- data.frame(
  lon = runif(100, min = -81.41094, max = -34.72999), 
  lat = runif(100, min = -55.61183, max = 12.4373)
) %>%  
  st_as_sf(coords = c("lon", "lat"), crs = sf::st_crs(4326))

# Get original bounding box
bbox <- st_bbox(points_sf) %>% 
  st_as_sfc() %>% 
  st_set_crs(st_crs(points_sf)) 

st_covers(bbox, points_sf) # Looks fine

# Expand the bounding box by factor 1
bbox2 <- st_bbox(points_sf) %>% 
  expand_bbox(expand_factor = 1) %>% 
  st_as_sfc() %>% 
  st_set_crs(st_crs(points_sf)) 

st_covers(bbox2, points_sf) # Still fine

# Expand the bounding box by factor 2
bbox3 <- st_bbox(points_sf) %>% 
  expand_bbox(expand_factor = 2) %>% 
  st_as_sfc() %>% 
  st_set_crs(st_crs(points_sf)) 

st_covers(bbox3, points_sf) # Unexpected behavior

# Plot the results
plot(bbox3, col = "red")
plot(bbox2, col = "orange", add = TRUE)
plot(bbox, col = "yellow", add = TRUE)
plot(points_sf, add = TRUE)

With an expansion factor of 1 (orange), everything works as expected. However, when I increase the expansion factor to 2 (red), st_covers() returns empty results, and the coverage test seems to break.

Other predicate functions like st_contains() or st_within() behave similarly. Any advice is greatly appreciated.

EDIT: I figured out that the break happens between an expand_factor of 1.44 and 1.45. In the original dataset, it was somewhere between 0.2 and 0.25, so it's not a fixed threshold.

Regarding your last edit, this is where bbox longitude range exceeds 180 degrees and this can introduce somewhat unexpected effects when using sf with s2 (geometries on sphere, enabled by default). Try switching it off with sf::sf_use_s2(use_s2 = FALSE). And this is what you are probably after when using binary predicates with bounding boxes and geographical coordinates (somewhat related)

180 vs 181 bbox longitude range:

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

p0   <- st_as_sfc("POINT(0 0)",   crs = "WGS84")
p100 <- st_as_sfc("POINT(100 0)", crs = "WGS84")

# x range -90 .. 90, covers (0 0) but not (100 0) -- expected
bbox_le180 <- 
  st_bbox(c(xmin = -90, ymin = -10, xmax = 90, ymax = 10), crs = "WGS84") |> 
  st_as_sfc() 
st_covers(bbox_le180, c(p0, p100), sparse = FALSE)
#>      [,1]  [,2]
#> [1,] TRUE FALSE


# x range -90 .. 91, does not cover (0 0) but covers (100 0) -- unexpected
bbox_gt180 <- 
  st_bbox(c(xmin = -90, ymin = -10, xmax = 91, ymax = 10), crs = "WGS84") |> 
  st_as_sfc() 
st_covers(bbox_gt180, c(p0, p100), sparse = FALSE)
#>       [,1] [,2]
#> [1,] FALSE TRUE

sf_use_s2(use_s2 = FALSE)
#> Spherical geometry (s2) switched off
# x range -90 .. 91, covers (0 0) but not (100 0) -- expected
st_covers(bbox_gt180, c(p0, p100), sparse = FALSE)
#> although coordinates are longitude/latitude, st_covers assumes that they are
#> planar
#>      [,1]  [,2]
#> [1,] TRUE FALSE

When checking WKT of s2 geometries, vertex order is reversed for bbox_gt180, this [ring direction] defines interior / exterior of the polygon. By default s2 checks polygon ring direction and "corrects" it when polygon area becomes bigger than half the globe, so bbox becomes "hole".

For bbox_le180, GEOS & s2 WKTs are the same:

st_as_text(bbox_le180)
#> [1] "POLYGON ((-90 -10, 90 -10, 90 10, -90 10, -90 -10))"
st_as_s2(bbox_le180) |> s2::s2_as_text()
#> [1] "POLYGON ((-90 -10, 90 -10, 90 10, -90 10, -90 -10))"

For bbox_gt180, GEOS & s2 WKTs are different (vertex order reversed):

st_as_text(bbox_gt180)
#> [1] "POLYGON ((-90 -10, 91 -10, 91 10, -90 10, -90 -10))"
st_as_s2(bbox_gt180) |> s2::s2_as_text()
#> [1] "POLYGON ((-90 10, 91 10, 91 -10, -90 -10, -90 10))"

There's also s2_oriented option to disable ring direction check:

sf_use_s2(use_s2 = TRUE)
#> Spherical geometry (s2) switched on
withr::with_options(
  list(s2_oriented = TRUE),
  {
    st_as_s2(bbox_gt180) |> s2::s2_as_text() |> print()
    st_covers(bbox_gt180, c(p0, p100), sparse = FALSE)
    
  }
)
#> [1] "POLYGON ((-90 -10, 91 -10, 91 10, -90 10, -90 -10))"
#>      [,1]  [,2]
#> [1,] TRUE FALSE

There are more details and examples in sf s2geometry Vignette - https://r-spatial.github.io/sf/articles/sf7.html#polygons-on-s2-divide-the-sphere-in-two-parts

^{Created on 2024-10-14 with reprex v2.1.1}