在R中从多边形中移除跨越国际日期线的线(例如俄罗斯在稀土中)

Question

问题:跨越国际日期线的多边形经常有一条南北线穿过它们.俄罗斯东部的稀土包裹就是一个很好的例子，但我也遇到过其他空间数据。我希望能够删除这条线来作图。

尝试：I主要使用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 (326 is或327 is)代码代替。

长答案

我的直觉是，这与将地理(长lat)数据投影到平面的挑战有关--或者是投影CRS，或者更简单地说是RStudio中的绘图查看器面板的平面。

我们知道，在地球的椭球模型中，-179和+179经度之间的(最小)地面距离与-1到+1之间的距离相同，距离为2度。然而，从数字的角度来看，这两条经线之间的距离是358度。

假设你是一个外星人(或一个扁平地球)，看着下面的投影 of world，而你不知道地球是椭球形的(或者你不知道这是一个投影)。如果你认为要从俄罗斯的一个地区(红色)进入另一个地区，你就得淋湿，这是可以原谅的。我想在默认情况下，ggplot是一个平地者.

​

​

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

​

​

所以，我的猜测是，当你尝试使用st_union()或st_simplify()时，在空间中，多边形实际上并不是相邻的，所以没有连接。这是投影CRS应该解决问题的地方，将和弦转换为相对于原点的值(long 0，lat 0)。

我认为这是给你带来麻烦的原因之一--EPSG的一个快速谷歌( google )：3338表示，这对阿拉斯加是件好事，但没有提到俄罗斯。当我在谷歌上搜索“utm俄罗斯”时，第一件事是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在地球两端有两组不同的同弦，在两者之间有一个很大的断裂(跨越x = 0)。投影32635，但只有一组和弦，这意味着俄罗斯的两个部分，根据这一投影，是数字位置相邻的对方。投影32635之所以有效，是因为它将和弦转换为‘(最小值？)距离原点的距离；它的起源(不像长弦和弦)不是在世界的另一边，它不需要在世界各地走两个不同的方向来确定到国家两端的最小距离(这就是导致另外两个投影的经度和弦断裂的原因)。我对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()已经解散了日期线上的两个边界，将单个多边形的数量从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是真的，那么每个特征几何都是统一的。例如，这可以用于在使用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，lat 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()

​

​

奖金

我试着用tmap和leaflet绘制tmap，但是没有得到想要的结果。我认为这是因为这些包更喜欢地理(lon-lat)和弦；据我所知，leaflet只接受longlat格式，尽管tmap当然可以处理投影数据，但我猜想在tmap下它会将它(或类似的)转换成它的首选投影。如果您真的想要这个可视化(这里、这里和这里)，那么可以在上面的链接中找到解决方案。

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的主要特征:区域、方向或距离。在投射像俄罗斯这样大而极的东西时，这一点就更加明显了。希望这些选项之一适合您的问题。

Answer 2

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

# 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()

结果：

​