从 R 中跨越国际日期变更线的多边形中删除线(例如 rnaturalearth 中的俄罗斯)

Question

问题:跨越国际日期变更线的多边形通常有一条南北线穿过它们。 rnaturalearth 包中的俄罗斯东部就是一个很好的例子，但我也遇到过其他空间数据。我希望能够删除这条线进行绘图。
尝试:
我主要使用 R 中的 sf 包进行映射。我尝试了各种解决方案，包括 st_union、st_combine、st_wrap_dateline、st_remove_holes，以及使用其他包中的函数，如聚合、合并和 gUnaryUnion，但到目前为止我的努力都没有结果。
示例:以下代码使用流行的 rnaturalearth 包演示了俄罗斯沿国际日期变更线的问题线。

library(tidyverse)
library(rnaturalearth)
library(sf)

#Import data
world <- ne_countries(scale = "medium",
                       returnclass = "sf") 

#I use the Alaska albers projection for this map,
#limit extent (https://spatialreference.org/ref/epsg/nad83-alaska-albers/)
xmin <- -2255938
xmax <- 1646517
ymin <- 449981
ymax <- 2676986

#plot
ggplot()+
  geom_sf(data=world, color="black", size=1)+
  coord_sf(crs=3338)+
  xlim(c(xmin,xmax))+ylim(c(ymin,ymax))+
  theme_bw()

谢谢!

Answer 1

简答
EPSG:3338 是问题所在 - 请改用 UTM(326XX 或 327XX)代码。
长答案
我的直觉是，这与将地理(长纬度)数据投影到平面上的挑战有关 - 无论是投影的 CRS，还是更简单的 RStudio 中绘图查看器 Pane 的平面。
我们知道，在地球的椭球模型上，-179 和 +179 经度之间的(最小)地面距离与 -1 和 +1 之间的距离相同，为 2 度的距离。然而，从数字的角度来看，这两条经线之间的距离是 358 度。
想象你是一个外星人(或一个扁平地球)，看着下面的 投影的 world ，而且你不知道地球是椭圆形的(或者你不知道这是一个投影)。如果您认为要从俄罗斯的一个地区(红色)到达另一个地区，您将不得不弄湿，那是情有可原的。我猜默认情况下，ggplot是平地。

想象上图中的每个多边形都是一块拼图。在您的情节中，我猜您将原点设置为 EPSG:3338 ( coord_sf(crs = 3338) ) 的中心，我认为它位于阿拉斯加/加拿大的某个地方？ (我在这里猜测是因为我不使用这种表示法，而是我更喜欢在发送到 ggplot 之前转换数据)。无论如何，ggplot知道它应该重新排列它的“拼图”，所以经度 -179 和 +179 彼此相邻 - 但这纯粹是视觉上的，就像你的情节一样:

所以，我的猜测是，当您尝试使用 st_union() 时或 st_simplify() ，多边形在空间中实际上并不相邻，因此没有连接。这是投影 CRS 应该解决问题的地方，将坐标转换为相对于 (long 0, lat 0) 以外的原点的值。
我认为这对您来说是麻烦的一个来源 - EPSG:3338 的快速谷歌说它对阿拉斯加有好处，但没有提到俄罗斯。当我在谷歌上搜索“utm russia”时出现的第一件事是 EPSG:32635。因此，让我们看一下 EPSG 代码 4326(WGS84 longlat)、3338(NAD83 Alaska)和 32635 的经度值。

# pull out russia
world %>% 
  filter(
    str_detect(name_long, 'Russia')
  ) %>% 
  select(name_long, geometry) %>% 
  {. ->> russia}

# extract coords of each projection
russia %>% 
  st_transform(3338) %>% 
  {. ->> russia_3338} %>% 
  st_coordinates %>% 
  as_tibble %>% 
  select(X) %>% 
  mutate(
    crs = 'utm_3338'
  ) %>% 
  {. ->> russia_coords_3338}

russia %>% 
  st_transform(4326) %>% 
  {. ->> russia_4326} %>% 
  st_coordinates %>% 
  as_tibble %>% 
  select(X) %>% 
  mutate(
    crs = 'utm_4326'
  ) %>% 
  {. ->> russia_coords_4326}

russia %>% 
  st_transform(32635) %>% 
  {. ->> russia_32635} %>% 
  st_coordinates %>% 
  as_tibble %>% 
  select(X) %>% 
  mutate(
    crs = 'utm_32635'
  ) %>% 
  {. ->> russia_coords_32635}

让我们将它们组合起来，看看经度值的直方图

# inspect X coords on a histogram
bind_rows(
  russia_coords_3338,
  russia_coords_4326,
  russia_coords_32635,
) %>% 
  ggplot(aes(X))+
  geom_histogram()+
  facet_wrap(~crs, ncol = 1, scales = 'free')

因此，如您所见，投影 4326 和 3338 在地球的两端有 2 组不同的坐标，中间有一个很大的中断(跨越 x = 0)。但是，投影 32635 只有一组坐标，这表明根据此投影，俄罗斯的两个部分在数字上彼此相邻。投影 32635 有效，因为它将坐标转换为“(最小？)距原点的距离”；其起源(与长纬度坐标不同)不在世界的另一端，并且不需要绕地球走 2 个不同的方向来确定到该国任何一端的最小距离(这就是导致中断的原因)在其他 2 个投影的经度坐标中)。我对 EPSG:3338 的了解不够，无法解释为什么它也会这样做，但我怀疑是因为它以阿拉斯加为重点，所以他们没有考虑跨越第 180 条子午线。
如果我们绘制 russia_32635我们可以看到这些碎片彼此相邻，但请记住我们不信任 ggplot刚刚。当我们使用 st_simplify()此日期变更线(红色)消失，证明这 2 个多边形彼此相邻并且可以简化/合并。

ggplot()+
  geom_sf(data = russia_32635, colour = 'red')+
  geom_sf(data = russia_32635 %>% st_simplify, fill = NA)

st_simplify()解散了日期变更线上的 2 个边界，将我们的单个多边形数量从 100 减少到 98。

russia_32635 %>% 
  st_cast('POLYGON')

# Simple feature collection with 100 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 21006.08 ymin: 4772449 xmax: 6273473 ymax: 13233690
# Projected CRS: WGS 84 / UTM zone 35N


russia_32635 %>% 
  st_simplify %>% 
  st_cast('POLYGON')

# Simple feature collection with 98 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 21006.08 ymin: 4772449 xmax: 6273473 ymax: 13233690
# Projected CRS: WGS 84 / UTM zone 35N

或者，它看起来像 st_union(..., by_feature = TRUE)也有效 - 见 ?st_union :

If by_feature is TRUE each feature geometry is unioned. This can for instance be used to resolve internal boundaries after polygons were combined using st_combine.

russia_32635 %>% 
  st_union(by_feature = TRUE) %>% 
  st_cast('POLYGON')

# Simple feature collection with 98 features and 1 field
# Geometry type: POLYGON
# Dimension:     XY
# Bounding box:  xmin: 21006.08 ymin: 4772449 xmax: 6273473 ymax: 13233690
# Projected CRS: WGS 84 / UTM zone 35N

因此，从技术上讲，您的俄罗斯图没有日期变更线。我认为俄罗斯很难绘制，因为 a) 它靠近两极，并且 b) 它覆盖了如此广阔的区域，这意味着大多数预测将从该国的一端倾斜到另一端。
然而对我来说，将情节“北上”定位是有意义的。一种方法是制作您自己的“Mollweide”投影并将原点指定为俄罗斯的大致中心(lon 99，lat 65)。无 st_buffer(0) ，出于某种原因，这与日期线一起绘制(参见 here 和 here 的示例，和 section 6.5 here 的解释)。

my_proj <- '+proj=moll +lon_0=99 +lat_0=65 +units=m'

russia_32635 %>% 
  st_buffer(0) %>% 
  st_transform(crs(my_proj)) %>%
  st_simplify %>% 
  ggplot()+
  geom_sf()

奖金
我尝试绘图 russia_32635 %>% st_simplify与 tmap和 leaflet ，但没有得到想要的结果。我认为这是因为这些包更喜欢地理(经纬度)坐标； leaflet只接受 longlat据我所知，格式，尽管 tmap当然可以处理投影数据，我的猜测是在引擎盖下它将它(或类似)转换为它的首选投影。如果您真的想要这种可视化效果( here 、 here 和 here )，可以在与上述相同的链接中找到解决方法。

library(tmap)

russia_32635 %>% 
  st_simplify %>% 
  tm_shape()+
  tm_polygons()


library(leaflet)

russia_32635 %>% 
  st_simplify %>%
  st_transform(4326) %>% # because leaflet only works with longlat projections
  leaflet %>% 
  addTiles %>% 
  addPolygons()

最终，您在投影数据时只能保留 2/3 的主要特征:面积、方向或距离。当投影像俄罗斯这样大而两极的东西时，这一点更加明显。希望这些选项之一适合您的问题。