Why do some xlims and ylims produce this error in ggplot and sf?

Question

I am learning to use ggplot with spatial data using sf.

When I try to produce the following plot I get an error:

library(sf)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)

world <- ne_countries(scale = "medium", returnclass = "sf")

ggplot(data = world) +
   geom_sf() +
   coord_sf(xlim = c(-102.15, -74.12), ylim = c(0, 33.97), expand = FALSE)

# Error in st_cast.POINT(x[[1]], to, ...) : 
#  cannot create MULTILINESTRING from POINT

However, if I just very slightly adjust the ylims in either direction, it works!:

ggplot(data = world) +
   geom_sf() +
   coord_sf(xlim = c(-102.15, -74.12), ylim = c(0.01, 33.97), expand = FALSE)

or

ggplot(data = world) +
   geom_sf() +
   coord_sf(xlim = c(-102.15, -74.12), ylim = c(-0.01, 33.97), expand = FALSE)

So that's a fine if rather hacky solution, but I am wondering what the issue is? On the final plot, I notice a little islands creeping into the map (the galapagos?) Is the tip of this polygon somehow interfering with the first piece of code by looking like a point rather than a polygon? Or is something else going on, and how can I address is a bit more elegantly?

Answer 1

I don't believe that the problem you describe is caused by the Galapagos - as it is present at scale = "small" which omits these islands.

I have tried another source of the world country dataset (the one from natural earth is not spherically valid, so I wanted to rule this issue out). The error remains, and now I suspect it has something to do with the code in {ggplot2} - specifically in how it handles the equator (i.e. zero) as a limit when the default expand = TRUE is overriden by user (it might be worth filing an issue at their github).

In the meantime I suggest two possible approaches:

• crop your area of interest outside of your ggplot call (this may be preferable for a variety of reasons - see ggplot2 and sf: geom_sf_text within limits set by coord_sf)
• avoid the expand = FALSE in your coord_sf() call

For cropping of your world object at data level (and not at presentation level in your ggplot call) consider this code:

library(sf)
library(ggplot2)
library(rnaturalearth)
library(rnaturalearthdata)

world <- ne_countries(scale = "medium", returnclass = "sf") %>% 
  st_make_valid()

bounds <- matrix(c(-102.15, 0, 
                   -74.12, 0, 
                   -74.12, 33.97, 
                   -102.15, 33.97,
                   -102.15, 0), 
                 byrow = TRUE, ncol = 2) %>%
  list() %>% 
  st_polygon() %>% 
  st_sfc(crs = 4326) 

small_world <- st_intersection(world, bounds) 

ggplot(data = small_world) +
  geom_sf()