R Make choropleth map with Robinson projection

I'm trying to make a choropleth map with a Robinson projection, mapping means onto each country. But something in the projection goes very wrong.

Here's a reproducible example and excerpts of full outputs:

mean <- structure(list(NAME = "Argentina", region = 0.000294752622638517), row.names = c(NA, 
-1L), class = c("data.table", "data.frame"), .internal.selfref = <pointer: 0x0000019a5c98d390>)

> mean
              NAME         mean
  1:     Argentina 0.0002947526
  2:    Kazakhstan 0.0083771737
  3:      Pakistan 0.0023414069
  4: United States 0.0002993637
  5:        Brazil 0.0003627592
 ---                           
176:      Kiribati 0.0621971274
177:    Uzbekistan 0.0224870460
178:  Turkmenistan 0.0639485852
179:    Kyrgyzstan 0.0174771136
180:    Tajikistan 0.0387022176

And shapefile as sf object:

shp <- structure(list(FIPS = "AC", ISO2 = "AG", ISO3 = "ATG", UN = 28L, 
    NAME = "Antigua and Barbuda", AREA = 44L, POP2005 = 83039, 
    REGION = 19L, SUBREGION = 29L, LON = -61.783, LAT = 17.078, 
    geometry = structure(list(structure(list(list(structure(c(-123.373336, 
    -123.47612, -123.65834, -123.752228, -123.761124, -123.767228, 
    -123.771668, -123.774444, -123.782226, -123.774444, -123.768342, 
    -123.665558, -123.652786, -123.588898, -123.568344, -123.488342, 
    -123.348342, -123.34056, -123.337784, -123.332778, -123.335006, 
    -123.365006, -123.373336, 34.0488820000003, 33.979438, 33.9938880000001, 
    34.0338820000002, 34.0394439999999, 34.0472180000002, 34.0561060000001, 
    34.0661080000002, 34.1883319999999, 34.2105479999999, 34.2194439999998, 
    34.327774, 34.334442, 34.3266600000002, 34.3166659999999, 
    34.2744359999999, 34.1872180000001, 34.1805500000002, 34.169998, 
    34.0916600000002, 34.081108, 34.054992, 34.0488820000003), dim = c(23L, 
    2L))), list(structure(c(-123.458344, -123.462234, -123.46556, 
    -123.477784, -123.50389, -123.631118, -123.669448, -123.678894, 
    -123.685562, -123.695008, -123.706116, -123.713348, -123.747788, 
    -123.750564, -123.746124, -123.700562, -123.691116, -123.678344, 
    -123.574448, -123.56723, -123.48668, -123.48056, -123.47612, 
    -123.463348, -123.458344, 35.217216, 35.0944440000002, 35.082222, 
    35.0811080000002, 35.0988840000001, 35.1677700000001, 35.1772160000002, 
    35.173332, 35.1655500000002, 35.1616600000002, 35.1661080000001, 
    35.184998, 35.377778, 35.397216, 35.4077760000001, 35.4455500000003, 
    35.4499960000001, 35.4494400000002, 35.4011080000001, 35.3944400000001, 
    35.306106, 35.2983320000002, 35.2894440000002, 35.2499920000001, 
    35.217216), dim = c(25L, 2L)))), class = c("XY", "MULTIPOLYGON", 
    "sfg"))), class = c("sfc_MULTIPOLYGON", "sfc"), precision = 0, bbox = structure(c(xmin = -123.782226, 
    ymin = 33.979438, xmax = -123.332778, ymax = 35.4499960000001
    ), class = "bbox"), crs = structure(list(input = "WGS 84", 
        wkt = "GEOGCRS[\"WGS 84\",\n    DATUM[\"World Geodetic System 1984\",\n        ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n            LENGTHUNIT[\"metre\",1]]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"latitude\",north,\n            ORDER[1],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        AXIS[\"longitude\",east,\n            ORDER[2],\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n    ID[\"EPSG\",4326]]"), class = "crs"), n_empty = 0L)), row.names = c(NA, 
-1L), sf_column = "geometry", agr = structure(c(FIPS = NA_integer_, 
ISO2 = NA_integer_, ISO3 = NA_integer_, UN = NA_integer_, NAME = NA_integer_, 
AREA = NA_integer_, POP2005 = NA_integer_, REGION = NA_integer_, 
SUBREGION = NA_integer_, LON = NA_integer_, LAT = NA_integer_
), levels = c("constant", "aggregate", "identity"), class = "factor"), class = c("sf", 
"tbl_df", "tbl", "data.frame"))

> shp

> shp
Simple feature collection with 246 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -360 ymin: -180 xmax: 360 ymax: 167.2472
Geodetic CRS:  WGS 84
# A tibble: 246 × 12
   FIPS  ISO2  ISO3     UN NAME    AREA POP2005 REGION SUBREGION     LON   LAT
   <chr> <chr> <chr> <int> <chr>  <int>   <dbl>  <int>     <int>   <dbl> <dbl>
 1 AC    AG    ATG      28 Anti…     44  8.30e4     19        29  -61.8   17.1
 2 AG    DZ    DZA      12 Alge… 238174  3.29e7      2        15    2.63  28.2
 3 AJ    AZ    AZE      31 Azer…   8260  8.35e6    142       145   47.4   40.4
 4 AL    AL    ALB       8 Alba…   2740  3.15e6    150        39   20.1   41.1
 5 AM    AM    ARM      51 Arme…   2820  3.02e6    142       145   44.6   40.5
 6 AO    AO    AGO      24 Ango… 124670  1.61e7      2        17   17.5  -12.3
 7 AQ    AS    ASM      16 Amer…     20  6.41e4      9        61 -171.   -14.3
 8 AR    AR    ARG      32 Arge… 273669  3.87e7     19         5  -65.2  -35.4
 9 AS    AU    AUS      36 Aust… 768230  2.03e7      9        53  136.   -25.0
10 BA    BH    BHR      48 Bahr…     71  7.25e5    142       145   50.6   26.0
# ℹ 236 more rows
# ℹ 1 more variable: geometry <MULTIPOLYGON [°]>
# ℹ Use `print(n = ...)` to see more rows

Here's what I've tried:

library("sf")
shp = sf::read_sf("World Borders Shapefiles/WB.shp")  
shp = sf::st_transform(shp, CRS("ESRI:54030"))

sf <- merge(shp, mean, by="NAME")

sf %>% 
  ggplot(aes(fill = mean, color = mean)) +
  geom_sf() +
  coord_sf(crs = "+proj=robin +lon_0=0w", datum = NA)

As you see, I try in two places to specify the projection, but it still looks like this:

I also don't get it to work with ggplot. Any suggestions would be greatly appreciated!

Solution

It is very likely that your shp object is not centered in longitude = 0°; thus, that is the reason why some polygons seem to cross from -180° to 180°. You can use rnaturalearth package to obtain that info and then just simply make the same join and plot your data. Here is a solution with some dummy data.

library(sf)
library(rnaturalearth)
library(dplyr)

# Download countries polygons
shp <- ne_countries(type = "countries", 
                   scale = "medium",
                   returnclass = "sf") |>
  # transform to desired projection
  st_transform("ESRI:54030") |>
  # This is a dummy column used to fill in mean values
  mutate(mean = rnorm(241,0.5,0.1))

# Make plot
shp |>
  ggplot(aes(fill = mean, color = mean)) +
  geom_sf() +
  coord_sf(expand = FALSE)