ggspatial: geom_sf fails with certain xlim values in coord_sf

What I'm trying to do: plot spatial objects with ggspatial::geom_sf(), using coord_sf() or a combination of ggspatial::layer_spatial() and ggspatial::annotation_spatial() to specify the extent of the of the plot.

This ggplot/ggspatial behavior has been described in several posts, but the "solutions" have been just ad-hoc hacks that do nothing to ensure the issue doesn't re-occur. See: Setting limits for x and y using coord_sf after updating ggplot2 Error when plotting latitude and longitude points on US map in RStudio

This code works just fine:

require(sf); require(ggspatial); require(rnaturalearth)    
country_polygons <- st_as_sf(ne_countries())
ggplot() + geom_sf(data=country_polygons)

As does this:

ggplot() + geom_sf(data=country_polygons) + coord_sf(xlim=c(-100,100), ylim=c(-60,60))

But:

ggplot() + geom_sf(data=country_polygons) + coord_sf(xlim=c(-160,150), ylim=c(-60,60))

Results in the error: Error in st_cast.POINT(X[[i]], ...) : cannot create MULTILINESTRING from POINT

Obviously this should not be happening. I don't understand why geom_sf calls st_cast.POINT() since there should be no points in the input. I've tried with three different world maps (the GADM 3.6 shapefile, rworldmap::countriesLow and the rnaturalearth one), so it doesn't seem to be dataset specific.

Using these x and y limits values, even a simple sf point object fails to plot!

set.seed(80085)
tibble(Lon=runif(1000,-180,180),
       Lat=runif(1000,-90,90)) %>% 
  st_as_sf(coords=1:2, remove=F, crs=4326) -> random_points

random_points %>% 
      ggplot() + geom_sf() + coord_sf(xlim=c(-160,150), ylim=c(-60,60))

The error message is again "cannot create MULTILINESTRING from POINT" which I have no idea why it would try to create a MULTILINESTRING.

UPDATE: The solution suggested here doesn't work: Why do some xlims and ylims produce this error in ggplot and sf?

random_points %>% 
  st_crop(xmin=-160, xmax=150, ymin=-60, ymax=60) %>% 
  ggplot() + geom_sf()

Results, weirdly, in only the points between 150oE and 160oW, i.e. over the Pacific Ocean around the 180th longitude, being preserved. I tried using longitudes from 0 to 360 and swapping the xmin and xmax, to no avail.

Faulty behavior of st_crop aside, passing a correct cropped object gives the same error again:

random_points %>% 
  filter(Lon<150, Lon>-160, Lat>-60, Lat<60) %>% 
  ggplot() + geom_sf()

> Error in st_cast.POINT(x[[1]], to, ...) : 
  cannot create MULTILINESTRING from POINT
In addition: Warning message:
In st_cast.GEOMETRYCOLLECTION(X[[i]], ...) :
  only first part of geometrycollection is retained

END UPDATE

Is there a solution or a workaround?

I believe this has something to do with the s2 geometry engine behaving in an unexpected fashion.

TL;DR: consider using sf_use_s2(F) in your code.

Long version:

library(sf)
library(dplyr)
library(mapview)

set.seed(80085)

random_points <- tibble(Lon = runif(1000,-180,180),
                        Lat = runif(1000,-90,90)) %>% 
  st_as_sf(coords=1:2, remove=F, crs=4326) 


mapview(random_points)  # this is "truth"

area_of_interest <- matrix(c(-160, -60, 
                             150, -60, 
                             150, 60, 
                             -160, 60,
                             -160, -60), 
                           byrow = TRUE, 
                           ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc(crs = 4326) 

mapview(area_of_interest) # this is our area of interest

# an atribute of points - in AOI or not?
random_points$aoi <- st_contains(area_of_interest,
                                 random_points,
                                 sparse = F) %>% 
  t() %>% 
  c()

# a visual overview; this is _not_ expected!
# the topology of AOI was not applied correctly
mapview(random_points, zcol = "aoi")

# let us try turning S2 engine off, and fall back to good old GEOS
sf_use_s2(F)

# exactly the same code as before!!
random_points$aoi <- st_contains(area_of_interest,
                                 random_points,
                                 sparse = F) %>% 
  t() %>% 
  c()

# but the output behaves much better!
mapview(random_points, zcol = "aoi")

# a seletion of random points; cropped to area of interest
library(ggplot2)

random_points %>% 
  filter(aoi) %>% 
  ggplot() + geom_sf(pch = 4, color = "red")