fze*_*oni 0 r geospatial ggplot2 r-sf
我试图在中央子午线不同于0的Robinson投影中投影世界地图。根据此StackOverFlow线程,这应该很容易(尽管示例使用sp)。
这是我的可复制代码:
library(sf)
library(ggplot2)
library(rnaturalearth)
world <- ne_countries(scale = 'small', returnclass = 'sf')
# Notice +lon_0=180 instead of 0
world_robinson <- st_transform(world, crs = '+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
ggplot() +
geom_sf(data = world_robinson)
Run Code Online (Sandbox Code Playgroud)
这就是结果。多边形从投影的一侧向另一侧闭合。

尝试具有sp相同的效果。我还尝试了一个shapefile,该文件仅包含来自http://www.naturalearthdata.com/的海岸线多边形(无政治边界),效果相似。

我试图在Mac OS X和Ubuntu 18.04上的两个独立的R安装上运行我的代码段。
变换后,跨越子午线的多边形在整个地图上都被拉伸。解决此问题的一种方法是将这些多边形从中间拆分,以便所有多边形都完全位于线的西侧或东侧。
# define a long & slim polygon that overlaps the meridian line & set its CRS to match
# that of world
polygon <- st_polygon(x = list(rbind(c(-0.0001, 90),
c(0, 90),
c(0, -90),
c(-0.0001, -90),
c(-0.0001, 90)))) %>%
st_sfc() %>%
st_set_crs(4326)
# modify world dataset to remove overlapping portions with world's polygons
world2 <- world %>% st_difference(polygon)
# perform transformation on modified version of world dataset
world_robinson <- st_transform(world2,
crs = '+proj=robin +lon_0=180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
# plot
ggplot() +
geom_sf(data = world_robinson)
Run Code Online (Sandbox Code Playgroud)
这是Z.lin答案的扩展(即,首先使用该答案来计算world_robinson)。但是,可以添加另一个有用的步骤。投影后,由多个多边形组成的区域,因为它们在原始投影中从地图的一侧穿过到另一侧(请参见南极洲,斐济和俄罗斯),因此在重新投影后仍会发生这种分裂。例如,这是南极洲的近景,我们可以看到它在本初子午线上有一个边界,不应有任何边界:
为了将这些区域重新缝合在一起,我们首先可以通过找到与本初子午线交叉的多边形来找出哪些多边形是问题:
bbox = st_bbox(world_robinson)
bbox[c(1,3)] = c(-1e-5,1e-5)
polygon2 <- st_as_sfc(bbox)
crosses = world_robinson %>%
st_intersects(polygon2) %>%
sapply(length) %>%
as.logical %>%
which
Run Code Online (Sandbox Code Playgroud)
现在我们可以选择这些多边形并将其缓冲区大小设置为零:
library(magrittr)
world_robinson[crosses,] %<>%
st_buffer(0)
ggplot(world_robinson) + geom_sf()
Run Code Online (Sandbox Code Playgroud)
如我们所见,地图不再分割本初子午线:
| 归档时间: |
|
| 查看次数: |
99 次 |
| 最近记录: |