从 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

我觉得我取得了重大进步，所以我发帖了，但这不是一个完整的答案。

# This is the portion containing the international dateline
df <- world[184, ]
# Split MULTIPOLYGON into individuals
df2 <- st_cast(df, "POLYGON")
# The little blob at the top is in df2[36, ] and df[38, ]
# Simplify it with the right tolerance and the line is gone
ggplot()+
  geom_sf(data=st_simplify(st_union(df2[36, ], df2[38, ]), dTolerance = 2), color="black", size=1)+
  coord_sf(crs=3338)+
  xlim(c(xmin,xmax))+ylim(c(ymin,ymax))+
  theme_bw()

结果：

Answer 2

简答

EPSG:3338 是问题所在 - 请改用 UTM（326XX 或 327XX）代码。

长答案

我的直觉是这与将地理（经纬度）数据投影到平面上的挑战有关——投影的 CRS，或者更简单地说是 RStudio 中绘图查看器窗格的平面。

我们知道，在地球的椭圆体模型中，经度 -179 和 +179 之间的（最小）地面距离与 -1 和 +1 之间的距离相同，相差 2 度。不过从数值上看，这两条经线之间的距离是358度。

想象你是一个外星人（或地平论者），看着 world 的投影，而你不知道地球是椭圆体形状（或者你不知道这是一个投影）。如果您认为要从俄罗斯的一部分（红色）到达另一部分，就必须弄湿，这是情有可原的。我想默认情况下，ggplot 是一个扁平地球。

假设上图中的每个多边形都是一块拼图。在您的情节中，我猜您将原点设置为 EPSG:3338 (coord_sf(crs = 3338)) 的中心，我认为它位于 Alaska/Canada 中的某处？（我在这里猜测是因为我不使用这种表示法，而是我更喜欢在发送到 ggplot 之前转换数据）。无论如何，ggplot 知道它应该重新排列它的 'puzzle pieces'，所以经度 -179 和 +179 彼此相邻 - 但这纯粹是视觉上的，如您的情节：

所以，我的猜测是，当您尝试使用 st_union() 或 st_simplify() 时，多边形在 space 中实际上并没有彼此相邻，因此没有连接。这是投影 CRS 应该解决问题的地方，将坐标转换为相对于 (long 0, lat 0) 以外的原点的值。

我认为这对您来说是一个麻烦源 - EPSG:3338 的快速 google 说这对阿拉斯加有好处，但没有提到俄罗斯。当我 googled '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 个投影的经度坐标中）。我对 EPSG:3338 的了解还不足以解释为什么它也这样做，但怀疑是因为它以阿拉斯加为中心，所以他们没有考虑穿越第 180 条子午线。

如果我们绘制 russia_32635 我们可以看到这些片段彼此相邻，但请记住我们还不信任 ggplot。当我们使用 st_simplify() 时，这条日期线（红色）消失了，证明 2 个多边形彼此相邻并且可以是 simplified/unioned.

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) 它覆盖了如此广阔的区域，这意味着大多数预测将从该国的一端倾斜到另一端。

但是对我来说，确定情节的方向是有意义的 'north-up'。一种方法是制作您自己的 'Mollweide' 投影并将原点分配给俄罗斯的大致中心（经度 99，纬度 65）。在没有 st_buffer(0) 的情况下，出于某种原因，此图带有日期线（请参阅 and for examples, and 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()

奖金

我尝试用 tmap 和 leaflet 绘制 russia_32635 %>% st_simplify，但没有得到想要的结果。我认为这是因为这些包更喜欢地理（lon-lat）坐标；据我所知，leaflet 只接受 longlat 格式，虽然 tmap 肯定可以处理投影数据，但我的猜测是在引擎盖下它会将它（或类似的）转换为它的首选投影。如果您真的想要这种可视化效果 (, and 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 的主要特征：面积、方向或距离。当投射像俄罗斯这样大而极的东西时，这一点变得更加明显。希望这些选项之一适合您的问题。

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

Remove line from polygon crossing the international dateline in R (e.g. Russia in rnaturalearth)

r

mapping

sf

简答

长答案

奖金