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

Question

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.

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.

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

This thread is related: https://gis.stackexchange.com/questions/78346/ortho-projection-produces-artifacts. Does anyone have suggestions on how to do that in R? — fzenoni, May 15 '19 at 14:02

score 9 · Accepted Answer · answered May 15 '19 at 15:24

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)

Thank you to both @Z.Lin and @dww. I assume that the `0` +/- epsilon value in the polygon definition can be dynamically defined as the `lon_0` value minus 180. Moreover, I still have an issue in the map display. Copy-pasting and executing your code I get the following: https://imgur.com/a/lhtWOpt. The `polygon` leftovers seem to be visible on the left-hand side of the map. Why do I see that and is there a way to remove it? — fzenoni, May 16 '19 at 07:20
This is really useful, thanks - I'm using the medium scale NE world map, and at the st_difference step I get an error: ```"Error in CPL_geos_op2(op, x, y) : Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 78.719726562500085 31.887646484374983 at 78.719726562500085 31.887646484374983."``` I think because the CRS is not 4326 - any ideas how to find the right CRS for the medium scale map? — mlcyo, Mar 09 '21 at 04:52

dww · Answer 2 · 2019-05-16T00:28:47.247

This is an extension to Z.lin's answer (i.e. use that answer first to calculate world_robinson). However, there is another useful step that can be added. After projecting, regions that were comprised of more than one polygon because they cross from one side of the map to the other in the original projection (see Antarctica, Fiji and Russia) still have this split after reprojection. For example, here is a close up of Antarctica where we can see that it has a boundary on the prime meridian where none should be:

To stitch these regions back togther, we can first find out which polygons are the problems by finding those that cross a the prime meridian:

bbox = st_bbox(world_robinson)
bbox[c(1,3)] = c(-1e-5,1e-5)
polygon2 <- st_as_sfc(bbox)

crosses = world_robinson %>%
  st_intersects(polygon2) %>%
  sapply(length) %>%
  as.logical %>%
  which

Now we can select those polygons and set their buffer size to zero:

library(magrittr)
world_robinson[crosses,] %<>%
  st_buffer(0) 

ggplot(world_robinson) + geom_sf()

As we can see, the map no longer has splits down the prime meridian:

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

2 Answers2

Linked