ggplot2（和sf）中世界地图的整个地球多边形

Posted 2023-02-16

【中文标题】ggplot2（和sf）中世界地图的整个地球多边形

【英文标题】：whole earth polygon for world map in ggplot2 (and sf)

【发布时间】：2017-08-29 16:42:18

【问题描述】：

在绘制世界地图时，ggplot2 存在问题：它用相同的颜色为整个背景着色，包括实际上不属于地球的图的角落，请参阅下面由以下生成的快照代码（它使用前沿 sf abd ggplot2 版本，但问题是一般性的，请参阅下面提到的博客文章）：

    #install.packages("devtools")
    #devtools::install_github("tidyverse/ggplot2")
    #devtools::install_github("edzer/sfr")

    library(ggplot2)
    library(sf)
    library(rnaturalearth)
    library(dplyr)

    theme_map <- function(...) 
      theme_minimal() +
      theme(
        text = element_text(family = "Ubuntu Regular", color = "#22211d"),
        axis.line = element_blank(),
        axis.text.x = element_blank(),
        axis.text.y = element_blank(),
        axis.ticks = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        panel.grid.minor = element_line(color = "#ebebe5", size = 0.2),
        panel.grid.major = element_line(color = "#ebebe5", size = 0.2),
        plot.background = element_rect(fill = "#f5f5f2", color = NA),
        panel.background = element_rect(fill = "#f5f5f2", color = NA),
        legend.background = element_rect(fill = "#f5f5f2", color = NA),
        panel.border = element_blank(),
        ...
      )
    

    crs <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"
    ctrys50m <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
       select(iso_a3, iso_n3, admin)

    ggplot() + 
      geom_sf(data = ctrys50m, alpha = 0.15, fill="grey") +
      coord_map() +
      coord_sf(crs = crs) +
      theme_map()

为了能够很好地绘制地球轮廓，在D3.js 中添加了一个特殊的GeoJSON type、type: "Sphere"，参见this thread 可以在操作中看到here：它是外部以下快照中的整个地球黑色边框：

我在R/ggplot2 中发现的唯一技巧是 Matt Strimas-Mackey 在他的博客条目 Mapping the Longest Commericial Flights in R 中发布的技巧，请参阅 Bounding box and graticules 部分和 @ 987654341@ 和 project_recenter 函数。

这些是相当多的代码，我想知道是否有一些 sf 或 geom_sf 代码会使代码更简洁/更简单，所以我尝试了：

    # whole world WSG84 bounding box
    sphere <- ne_download(category = "physical", type = "wgs84_bounding_box", returnclass = "sf")
    sphere_laea <- st_transform(sphere, crs)
    ggplot() + 
      geom_sf(data = sphere, fill = "#D8F4FF") +
      coord_sf(crs = crs) +
      geom_sf(data = ctrys50m, alpha = 0.15, fill="grey") +
      coord_map() +
      coord_sf(crs = crs) +
      theme_map()

我得到的只是一个额外的“反子午线”（注意从北极开始的线...），没有充满#D8F4FF...的海洋 而且底部的多边形非常不规则（D3.js 大师做了一些聪明的adaptive resampling 以提高投影线的准确性......）

关于我为 ggplot2 世界地图获取整个世界多边形的尝试有什么问题吗？ （感谢您阅读本文！）

【参考方案1】：

按照 Dave 的建议，我创建了一个足够小的新标线，以便凸包操作可以很好地逼近 地球的边界。 以下代码给出了非常好的结果（见后图）：

library(ggplot2)
library(sf)
library(rnaturalearth)
library(dplyr)

crs <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +datum=WGS84 +units=m +no_defs"

ctrys50m <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
  select(iso_a3, iso_n3, admin)

sphere <- st_graticule(ndiscr = 10000, margin = 10e-6) %>%
  st_transform(crs = crs) %>%
  st_convex_hull() %>%
  summarise(geometry = st_union(geometry))

ggplot()  +
  geom_sf(data = sphere, fill = "#D8F4FF", alpha = 0.7) +
  geom_sf(data = ctrys50m, fill="grey") +
  theme_bw()

【参考方案2】：

我还在探索sf 的神奇之处，所以可能会有更直接的解决方案来解决您的问题。不过，您可能会发现这种方法很有用：

为了提取地球球体的多边形，只需使用所需投影的经纬网，创建一个凸包并合并生成的多边形。

library(ggplot2)
library(sf)
library(rnaturalearth)
library(dplyr)

crs <- "+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000
        +datum=WGS84 +units=m +no_defs"

ctrys50m <- ne_countries(scale = 50, type = "countries", returnclass = "sf") %>%
  select(iso_a3, iso_n3, admin)

sphere <- st_graticule(st_transform(ctrys50m, crs = crs)) %>%
  st_convex_hull() %>%
  summarise(geometry = st_union(geometry))

ggplot()  +
  geom_sf(data = sphere, fill = "#D8F4FF", alpha = 0.7) +
  geom_sf(data = ctrys50m, fill="grey") +
  #coord_sf(crs = crs) + # not necessary because the crs from the first layer is being used.
  theme_bw()

这将为您提供以下情节，该情节仍然有一些可见的角落，但根据您的目的可能就足够了。

【讨论】：

这是一个绝妙的建议。我将使用ne_countries(scale = 10, ...) 来获得更好的多边形，否则它是一个可行的解决方案。谢谢！

我得到了非常好的结果，即没有可见的角落，sphere的定义如下：sphere <- st_graticule(ndiscr = 10000, margin = 10e-6) %>% st_transform(crs = crs) %>% st_convex_hull() %>% summarise(geometry = st_union(geometry))