rggplot2geospatialr-sf

Visual bug when changing Robinson projection's central meridian with ggplot2

I am attempting to project a world map in a Robinson projection where the central meridian is different from 0. According to this StackOverFlow thread, it should be an easy thing (albeit the example uses sp).

Here is my reproducible code:

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

world <- ne_countries(scale = 'small', returnclass = 'sf')

# Notice +lon_0=180 instead of 0
world_robinson <- st_transform(world, crs = '+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

ggplot() +
geom_sf(data = world_robinson)

This is the result. Polygons are closing themselves from one side to the other of the projection.

image

Trying with sp gives the same effect. I also tried with a shapefile including only polygons from coastlines (no political borders) from http://www.naturalearthdata.com/ and the effect is similar.

image

I tried to run my snippet on two independent R installations on Mac OS X and Ubuntu 18.04.

Solution

  • Polygons that straddle the meridian line get stretched all the way across the map, after the transformation. One way to get around this is to split these polygons down the middle, so that all polygons are either completely to the west or east of the line.

    # define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world
polygon <- st_polygon(x = list(rbind(c(-0.0001, 90),
                                     c(0, 90),
                                     c(0, -90),
                                     c(-0.0001, -90),
                                     c(-0.0001, 90)))) %>%
  st_sfc() %>%
  st_set_crs(4326)

# modify world dataset to remove overlapping portions with world's polygons
world2 <- world %>% st_difference(polygon)

# perform transformation on modified version of world dataset
world_robinson <- st_transform(world2, 
                               crs = '+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')

# plot
ggplot() +
  geom_sf(data = world_robinson)

    plot