How to rotate world map using Mollweide projection with sf/rnaturalearth/ggplot in R?

Question

I would like to plot a world map in R using the Mollweide projection centred on the Pacific region (specifically, Australia), using a rnaturalearth --> sf --> ggplot pipeline.

I have been running into the annoying issue of having connected lines across the globe.

From a fresh R session, I run

library(tidyverse)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)

target_crs <- st_crs("+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=133")
worldrn <- ne_countries(scale = "medium", returnclass = "sf") %>%
  sf::st_transform(crs = target_crs)
ggplot(data = worldrn, aes(group = admin)) +
  geom_sf()

which generates this plot

This is my sessionInfo():

R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] sf_0.9-8                rnaturalearthdata_0.1.0 rnaturalearth_0.1.0     forcats_0.5.1           stringr_1.4.0           dplyr_1.0.7             purrr_0.3.4             readr_1.4.0             tidyr_1.1.3            
[10] tibble_3.1.2            ggplot2_3.3.5           tidyverse_1.3.1        

loaded via a namespace (and not attached):
 [1] tidyselect_1.1.1   lattice_0.20-44    haven_2.4.1        colorspace_2.0-2   vctrs_0.3.8        generics_0.1.0     utf8_1.2.1         rlang_0.4.11       e1071_1.7-7        pillar_1.6.1       glue_1.4.2        
[12] withr_2.4.2        DBI_1.1.1          sp_1.4-5           dbplyr_2.1.1       modelr_0.1.8       readxl_1.3.1       lifecycle_1.0.0    rgeos_0.5-5        munsell_0.5.0      gtable_0.3.0       cellranger_1.1.0  
[23] rvest_1.0.0        class_7.3-19       fansi_0.5.0        broom_0.7.8        Rcpp_1.0.6         KernSmooth_2.23-20 scales_1.1.1       backports_1.2.1    classInt_0.4-3     jsonlite_1.7.2     farver_2.1.0      
[34] fs_1.5.0           hms_1.1.0          stringi_1.6.2      grid_4.1.0         cli_3.0.0          tools_4.1.0        magrittr_2.0.1     proxy_0.4-25       crayon_1.4.1       pkgconfig_2.0.3    ellipsis_0.3.2    
[45] xml2_1.3.2         reprex_2.0.0       lubridate_1.7.10   assertthat_0.2.1   httr_1.4.2         rstudioapi_0.13    R6_2.5.0           units_0.7-1        compiler_4.1.0

An alternative

Then, from a fresh session, I also tried the approach recommended here using the Robinson projection and GDAL, whereby one crops the globe, recentres the longitudes, and stitches it all back together. However, for some reason the connecting lines remain. By the way, in this second example I downloaded the country layers from Natural Earth

library(ggplot2)
library(dplyr)
library(ggthemes)
library(sp)
library(rgdal)

# assumes you are in the ne_110m... directory
# split the world and stitch it back together again

system("ogr2ogr world_part1.shp ne_110m_admin_0_countries.shp -clipsrc -180 -90 0 90")
system("ogr2ogr world_part2.shp ne_110m_admin_0_countries.shp -clipsrc 0 -90 180 90")
system('ogr2ogr world_part1_shifted.shp world_part1.shp -dialect sqlite -sql "SELECT ShiftCoords(geometry,360,0), admin FROM world_part1"')
system("ogr2ogr world_0_360_raw.shp world_part2.shp")
system("ogr2ogr -update -append world_0_360_raw.shp world_part1_shifted.shp -nln world_0_360_raw")

proj_str <- "+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
my_prj <- sp::CRS(proj_str)
world <- rgdal::readOGR("world_0_360_raw.shp", "world_0_360_raw") %>%
  sp::spTransform(my_prj) %>%
  ggplot2::fortify()

ggplot(data = world, mapping = aes(x = long, y = lat)) +
  geom_map(map = world, mapping = aes(map_id = id),
           fill = "khaki", colour = "black", size = 0.25) +
  coord_equal() +
  theme_map() +
  theme(plot.background = element_rect(fill = "azure2"))

Which generates this figure

sessionInfo() for this second example:

R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Catalina 10.15.7

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_AU.UTF-8/en_AU.UTF-8/en_AU.UTF-8/C/en_AU.UTF-8/en_AU.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] rgdal_1.5-23   sp_1.4-5       ggthemes_4.2.4 dplyr_1.0.7    ggplot2_3.3.5 

loaded via a namespace (and not attached):
 [1] magrittr_2.0.1   tidyselect_1.1.1 munsell_0.5.0    lattice_0.20-44  colorspace_2.0-2 R6_2.5.0         rlang_0.4.11     fansi_0.5.0      stringr_1.4.0    tools_4.1.0      grid_4.1.0       gtable_0.3.0    
[13] utf8_1.2.1       DBI_1.1.1        withr_2.4.2      ellipsis_0.3.2   assertthat_0.2.1 tibble_3.1.2     lifecycle_1.0.0  crayon_1.4.1     purrr_0.3.4      vctrs_0.3.8      glue_1.4.2       stringi_1.6.2   
[25] compiler_4.1.0   pillar_1.6.1     generics_0.1.0   scales_1.1.1     pkgconfig_2.0.3

I'd appreciate any help to get rid of these lines! Ideally having a Pacific-centred world map via rnaturalearth would be the easiest solution, but that does not seem to exist?

Yep, happy to take any suggestions. Is this solvable using `tmap`? I never used it though. — dbarneche, Jul 07 '21 at 08:53
Hmm, no `tmap` does not solve the problem which lies in the jump across the date line. If you search for similar issues you will find a lot of questions regarding that topic but no satisfactory solution. — mgrund, Jul 07 '21 at 09:43
Indeed, and that's why I decided to post here, particularly considering that the second approach (which apparently used to work in 2015) does not work anymore. — dbarneche, Jul 07 '21 at 12:44

dieghernan · Accepted Answer · 2021-07-09T08:55:45.253

See a potential solution that works, based on this answer:

library(tidyverse)
library(rnaturalearth)
library(rnaturalearthdata)
library(sf)
#> Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

target_crs <- st_crs("+proj=moll +x_0=0 +y_0=0 +lat_0=0 +lon_0=133")


worldrn <- ne_countries(scale = "medium", returnclass = "sf") %>%
  st_make_valid()


# define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world

# Centered in lon 133

offset <- 180 - 133


polygon <- st_polygon(x = list(rbind(
  c(-0.0001 - offset, 90),
  c(0 - offset, 90),
  c(0 - offset, -90),
  c(-0.0001 - offset, -90),
  c(-0.0001 - offset, 90)
))) %>%
  st_sfc() %>%
  st_set_crs(4326)


# modify world dataset to remove overlapping portions with world's polygons
world2 <- worldrn %>% st_difference(polygon)
#> Warning: attribute variables are assumed to be spatially constant throughout all
#> geometries


# Transform
world3 <- world2 %>% st_transform(crs = target_crs)
ggplot(data = world3, aes(group = admin)) +
  geom_sf(fill = "red")

sessionInfo()
#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19041)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=Spanish_Spain.1252  LC_CTYPE=Spanish_Spain.1252   
#> [3] LC_MONETARY=Spanish_Spain.1252 LC_NUMERIC=C                  
#> [5] LC_TIME=Spanish_Spain.1252    
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] sf_1.0-1                rnaturalearthdata_0.1.0 rnaturalearth_0.1.0    
#>  [4] forcats_0.5.1           stringr_1.4.0           dplyr_1.0.7            
#>  [7] purrr_0.3.4             readr_1.4.0             tidyr_1.1.3            
#> [10] tibble_3.1.2            ggplot2_3.3.5           tidyverse_1.3.1        
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.7         lattice_0.20-44    lubridate_1.7.10   class_7.3-19      
#>  [5] assertthat_0.2.1   digest_0.6.27      utf8_1.2.1         R6_2.5.0          
#>  [9] cellranger_1.1.0   backports_1.2.1    reprex_2.0.0       evaluate_0.14     
#> [13] e1071_1.7-7        httr_1.4.2         highr_0.9          pillar_1.6.1      
#> [17] rlang_0.4.11       readxl_1.3.1       rstudioapi_0.13    rmarkdown_2.9     
#> [21] styler_1.4.1       munsell_0.5.0      proxy_0.4-26       broom_0.7.8       
#> [25] compiler_4.1.0     modelr_0.1.8       xfun_0.24          pkgconfig_2.0.3   
#> [29] rgeos_0.5-5        htmltools_0.5.1.1  tidyselect_1.1.1   fansi_0.5.0       
#> [33] crayon_1.4.1       dbplyr_2.1.1       withr_2.4.2        wk_0.4.1          
#> [37] grid_4.1.0         jsonlite_1.7.2     gtable_0.3.0       lifecycle_1.0.0   
#> [41] DBI_1.1.1          magrittr_2.0.1     units_0.7-2        scales_1.1.1      
#> [45] KernSmooth_2.23-20 cli_3.0.0          stringi_1.6.2      farver_2.1.0      
#> [49] fs_1.5.0           sp_1.4-5           xml2_1.3.2         ellipsis_0.3.2    
#> [53] generics_0.1.0     vctrs_0.3.8        s2_1.0.6           tools_4.1.0       
#> [57] glue_1.4.2         hms_1.1.0          yaml_2.2.1         colorspace_2.0-2  
#> [61] classInt_0.4-3     rvest_1.0.0        knitr_1.33         haven_2.4.1

^{Created on 2021-07-09 by the reprex package (v2.0.0)}

Thanks so much @dieghernan ! It worked perfectly on my end as well. I can also confirm that it is quite simple to add points to this new map. this worked for me: ``` coords <- data.frame(lon = 133, lat = -21) %>% sf::st_as_sf(coords = c("lon", "lat"), crs = 4326) %>% sf::st_transform(crs = target_crs) ggplot() + geom_sf(data = world3, aes(group = admin), fill = "red") + geom_sf(data = coords) ``` — dbarneche, Jul 15 '21 at 03:06

How to rotate world map using Mollweide projection with sf/rnaturalearth/ggplot in R?

An alternative

1 Answers1

Linked