如何在 R 中使用 Mollweide 投影和 sf/rnaturalearth/ggplot 旋转世界地图？
How to rotate world map using Mollweide projection with sf/rnaturalearth/ggplot in R?
我想使用以太平洋地区（特别是澳大利亚）为中心的 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
另一种选择
然后，从一个新的会话开始，我也尝试了推荐的方法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"))
生成这个数字
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 获得以太平洋为中心的世界地图将是最简单的解决方案，但这似乎不存在？
查看可行的潜在解决方案，基于 :
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 创建
我想使用以太平洋地区（特别是澳大利亚）为中心的 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
另一种选择
然后，从一个新的会话开始，我也尝试了推荐的方法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"))
生成这个数字
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 获得以太平洋为中心的世界地图将是最简单的解决方案，但这似乎不存在？
查看可行的潜在解决方案，基于
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 创建