Combine sf geometry with transformed geometry

I am replicating xkcd's "Bad Map Projection: ABS(Longitude)" with R. I have a decent start but I'd like to do it the "right" ways.

Here, I take a world map and make a copy with a horizontal reflection. Then I plot both and crop the viewing window to positive longitude.

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

world <- ne_countries(scale = "small", returnclass = "sf") |>
  sf::st_as_sf() %>% 
  sf::st_set_crs(value = 4326) 

# matrix multiplication to reflect longitude
world_rev <- st_geometry(world)  * matrix(c(-1,0,0,1), 2, 2) 
world2 <- world |>
  st_set_geometry(world_rev) |>
  sf::st_as_sf() |>
  sf::st_set_crs(value = 4326) 

ggplot() +
  geom_sf(data = world, fill = "#336666", alpha = 0.5) +
  geom_sf(data = world2, fill = "#663366", alpha = 0.5) +
  coord_sf(xlim = c(0, 180), expand = 0)

A shortfall of this is that I can't seem to apply different projections. For instance, if I end with

coord_sf(crs= "+proj=robin", xlim = c(0, 180), expand = 0)

I get a blank plot. Is there a more robust way to create the reversed object?

Somewhat closer, though still not ideal:

• move longitude of Robinson projection center from default 0 to 90
• to handle "over the edge" geometries, crop those to 0 ... 180 longitude
• for this to work, disable s2

library(ggplot2)
library(sf)
#> Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(rnaturalearth)
library(rnaturalearthdata)

world <- ne_countries(scale = "small", returnclass = "sf") |>
  sf::st_as_sf() %>% 
  sf::st_set_crs(value = 4326) 

# matrix multiplication to reflect longitude
world_rev <- st_geometry(world)  * matrix(c(-1,0,0,1), 2, 2) 

world2 <- world |>
  st_set_geometry(world_rev) |>
  sf::st_as_sf() |>
  sf::st_set_crs(value = 4326) 

sf_use_s2(use_s2 = FALSE)
#> Spherical geometry (s2) switched off
bbox_0_180 <- st_bbox(c(xmin = 0, ymin = -90, xmax = 180, ymax = 90), crs = 4326)
world_crp  <- st_crop(world,  bbox_0_180)
world2_crp <- st_crop(world2, bbox_0_180)

ggplot() +
  geom_sf(data = world_crp,  fill = "#336666", alpha = 0.5) +
  geom_sf(data = world2_crp, fill = "#663366", alpha = 0.5) +
  coord_sf(crs= "+proj=robin +lon_0=90", expand = FALSE) +
  scale_x_continuous(breaks = seq(0,180,30))

<!-- -->

^{Created on 2023-07-27 with reprex v2.0.2}

You got blank plot with

coord_sf(crs= "+proj=robin", xlim = c(0, 180), expand = 0)

because Robinson is projected and units are meters, applies for xlim too, you could try this instead:

coord_sf(crs= "+proj=robin", xlim = c(-18e6, 18e6), expand = 0)

, with world and world2 would produce something like: