rggplot2r-sfrgdalproj

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


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

enter image description here

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

enter image description here

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?


Solution

  • 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)