从 R 中与国际日期变更线交叉的多边形中删除线（例如 rnaturalearth 中的俄罗斯）

Question

问题：跨越国际日期变更线的多边形经常有南北线穿过。rnaturalearth 包中的俄罗斯东部就是一个很好的例子，但我在其他空间数据中也遇到过这种情况。我希望能够删除这条线以进行绘图。

尝试： 我主要使用R中的sf包进行映射。我尝试了各种解决方案，包括 st_union、st_combine、st_wrap_dateline、st_remove_holes，以及使用其他包中的函数（例如aggregate、merge 和 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 中绘图查看器窗格的平面。

我们知道，在地球的椭球模型上，经度 -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 阿拉斯加）和 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 个投影的经度坐标）。我对 EPSG:3338 的了解不够，无法解释为什么它也会这样做，但怀疑这是因为它以阿拉斯加为中心，所以他们没有考虑穿越 180 度子午线。

如果我们绘图，russia_32635我们可以看到这些片段彼此相邻，但请记住我们还不信任ggplot。当我们使用st_simplify()这条日期变更线（红色）时，它就会消失，证明两个多边形彼此相邻并且可以简化/合并。

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：

如果by_feature为 TRUE，则每个要素几何图形都会合并。例如，这可以用于在使用 组合多边形后解析内部边界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”投影，并将原点指定为俄罗斯的大致中心（经度 99，纬度 65）。如果没有st_buffer(0)，则由于某种原因，这会与日期线一起绘制（请参阅此处和此处的示例，以及此处第 6.5 节的解释）。

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仅接受格式，虽然肯定可以处理投影数据，但我的猜测是，它会在引擎盖下将其（或类似的）转换为首选投影。如果您确实想要此可视化（此处、此处和此处），则可以在与上述相同的链接中找到解决方法。longlattmap

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 主要特征：面积、方向或距离。当预测像俄罗斯这样的大国和极地国家时，这一点变得更加明显。希望这些选项之一适合您的问题。