时间戳

首页 > 解决方案 > 从多边形中删除与 R 中国际日期变更线相交的线（例如俄罗斯在 rnaturalearth）

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

问题描述

问题：跨越国际日期变更线的多边形经常有一条南北线穿过它们。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()

谢谢！

标签： rmappingsf

解决方案

简短的回答

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()orst_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 在地球两端有两组不同的坐标，有一个很大的中断（跨越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：

如果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”投影并将原点分配到俄罗斯的大致中心（lon 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()

在此处输入图像描述

奖金

russia_32635 %>% st_simplify我尝试使用tmapand进行绘图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 的主要特征：面积、方向或距离。当投影像俄罗斯这样大而极地的东西时，这一点变得更加明显。希望这些选项之一适合您的问题。

推荐阅读