如何在 R 中使用 Mollweide 投影和 sf/rnaturalearth/ggplot 来旋转世界地图？

Question

我想使用以太平洋地区（特别是澳大利亚）为中心的 Mollweide 投影在 R 中绘制世界地图，使用

rnaturalearth

-->

sf

-->

ggplot

管道。

我遇到了一个烦人的问题，即在全球范围内连接线路。

从新的 R 会话中，我运行

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

生成此图

这是我的

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    

替代方案

然后，在新的会议中，我还尝试了使用罗宾逊投影和 GDAL 推荐的方法，即裁剪地球，重新调整经度，然后将其全部缝合在一起。然而，由于某种原因，连接线仍然存在。顺便说一句，在第二个示例中，我从 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"))

生成这个数字

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] 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 

如果您能帮助我摆脱这些线条，我将不胜感激！理想情况下，通过 

rnaturalearth

拥有以太平洋为中心的世界地图将是最简单的解决方案，但这似乎不存在？

Answer 1

此答案

查看可行的潜在解决方案：

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

由

reprex 包于 2021-07-09 创建^（v2.0.0）

Answer 2

你能帮我解决一下吗？

library("rnaturalearth") # World shapefile
#> Warning: il pacchetto 'rnaturalearth' è stato creato con R versione 4.3.2
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, will retire in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> The sp package is now running under evolution status 2
#>      (status 2 uses the sf package in place of rgdal)
#> Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library("sf")
#> Warning: il pacchetto 'sf' è stato creato con R versione 4.3.2
#> Linking to GEOS 3.11.2, GDAL 3.7.2, PROJ 9.3.0; sf_use_s2() is TRUE

world <- ne_countries(scale = 110, returnclass = "sf") |> 
  st_make_valid()

offset <- 180 


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) |> 
  st_make_valid()
#> Warning in st_is_longlat(x): bounding box has potentially an invalid value
#> range for longlat data


# modify world dataset to remove overlapping portions with world's polygons
world <- world |>  st_difference(polygon)
#> Error in wk_handle.wk_wkb(wkb, s2_geography_writer(oriented = oriented, : Loop 0 is not valid: Edge 0 crosses edge 78

创建于 2023-12-04，使用

reprex v2.0.2