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

Posted 2023-03-23

【中文标题】如何在 R 中使用带有 sf/rnaturalearth/ggplot 的 Mollweide 投影旋转世界地图？

【英文标题】：How to rotate world map using Mollweide projection with sf/rnaturalearth/ggplot in R?

【发布时间】：2021-09-17 13:43:52

【问题描述】：

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

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

从新的 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    

另一种选择

然后，在一个新的会话中，我还尝试了使用 Robinson 投影和 GDAL 推荐的方法here，通过这种方法可以裁剪地球，重新定位经度，然后将它们缝合在一起。但是，由于某种原因，连接线仍然存在。顺便说一句，在第二个示例中，我从 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 获得以太平洋为中心的世界地图将是最简单的解决方案，但这似乎不存在？

【问题讨论】：

您是否可以选择使用tmap？

是的，很高兴接受任何建议。这可以使用tmap 解决吗？不过我从来没用过。

嗯，不，tmap 不能解决跨越日期变更线的问题。如果您搜索类似的问题，您会发现很多关于该主题的问题，但没有令人满意的解决方案。

确实如此，这就是我决定在这里发帖的原因，特别是考虑到第二种方法（显然在 2015 年曾经有效）不再有效。

【参考方案1】：

查看可行的潜在解决方案，基于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

^{由reprex package (v2.0.0) 于 2021-07-09 创建}

【讨论】：

非常感谢@dieghernan！它对我来说也很完美。我还可以确认向这个新地图添加点非常简单。这对我有用：`` coords % 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) ```