recenter projections over Pacific with sf package

Question

I'm trying to create projected maps (Robinson in the example below) centered on the Pacific, but I'm getting strange unwanted horizontal lines. The st_bbox() solution doesn't seem to help

rm(list=ls())
library(giscoR)
library(rworldmap)
library(tidyverse)
library(sf)

crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs"
crsrobin <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

# Using giscoR
world_sf <- gisco_get_countries(resolution = 10)

ggplot() +
  geom_sf(data = world_sf)

world_robinson <- st_transform(world_sf, 
                               crs = crsrobin) ## looks great
ggplot() +
  geom_sf(data = world_robinson) ## looks great

## Pacific centered Robinson
crsrobin <- "+proj=robin +lon_0=-180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
world_robinson <- st_transform(world_sf, 
                               crs = crsrobin) ## horizontal artifacts

ggplot() +
  geom_sf(data = world_robinson) ## horizontal artifacts

I get the same result with different map data

world_sp <- fortify(getMap())
world_sf <- world_sp %>% st_as_sf(coords = c("long", "lat"), crs = crsLONGLAT, row.names="group") %>%
  group_by(group) %>% summarise(geometry = st_combine(geometry)) %>% st_cast("POLYGON")

This bbox hack looks fine for coordinates, but doesn't work with the projection

world_box <- st_cast(st_combine(st_as_sf(data.frame("x"=c(-180, -180, 180, 180),"y"=c(90, -90, -90, 90)),coords = c("x","y"), crs=crsLONGLAT)),"POLYGON")
world_box_robin <- st_transform(world_box, crs=crsrobin)

ggplot() + geom_sf(data=world_box) + geom_sf(data=world_sf)

b <- st_bbox(world_box)
ggplot() + geom_sf(data=world_sf) +
  coord_sf(crs = crsLONGLAT, xlim = c(b["xmin"], b["xmax"]), ylim = c(b["ymin"], b["ymax"]))


ggplot()+ geom_sf(data=world_box_robin)  + geom_sf(data=world_robinson) 


b <- st_bbox(world_box_robin)
ggplot() + geom_sf(data=world_robinson) + 
  coord_sf(crs = crsrobin, xlim = c(b["xmin"], b["xmax"]), ylim = c(b["ymin"], b["ymax"]))

What am I missing?

Answer 1

You might be looking for st_break_antimeridian:

crsrobin <- "+proj=robin +lon_0=-180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

world_robinson2 <- world_sf %>%
  st_break_antimeridian(lon_0 = 180) %>% # insert this before transformation
                                         # (& ignore the warning)
  st_transform(crs = crsrobin)

ggplot() +
  geom_sf(data = world_robinson2)