使用sf包制作太平洋中心投影地图时出现异常水平线问题

尝试制作以太平洋为中心的Robinson投影地图，但出现奇怪的多余水平线，st_bbox()解决方案无效。更换不同地图数据后仍得到相同结果，bbox技巧在地理坐标系下可行，但在投影坐标系中无法解决问题。请问忽略了什么？

原始代码 #

rm(list=ls())
library(giscoR)
library(rworldmap)
library(tidyverse)
library(sf)

crsLONGLAT <- "+proj=longlat +datum=WGS84 +no_defs"
crsrobin <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"

# 使用giscoR获取数据
world_sf <- gisco_get_countries(resolution = 10)

# 查看原始地理坐标系地图
ggplot() +
  geom_sf(data = world_sf)

# 常规Robinson投影（中央经线0°），显示正常
world_robinson <- st_transform(world_sf, 
                               crs = crsrobin)
ggplot() +
  geom_sf(data = world_robinson)

# 太平洋中心Robinson投影（中央经线-180°），出现水平伪影
crsrobin <- "+proj=robin +lon_0=-180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
world_robinson <- st_transform(world_sf, 
                               crs = crsrobin)

ggplot() +
  geom_sf(data = world_robinson)

更换地图数据的代码 #

world_sp <- fortify(getMap())
world_sf <- world_sp %>% st_as_sf(coords = c("long", "lat"), crs = crsLONGLAT, row.names="group") %>%
  group_by(group) %>% summarise(geometry = st_combine(geometry)) %>% st_cast("POLYGON")

bbox相关尝试代码 #

world_box <- st_cast(st_combine(st_as_sf(data.frame("x"=c(-180, -180, 180, 180),"y"=c(90, -90, -90, 90)),coords = c("x","y"), crs=crsLONGLAT)),"POLYGON")
world_box_robin <- st_transform(world_box, crs=crsrobin)

# 地理坐标系下的bbox限制
ggplot() + geom_sf(data=world_box) + geom_sf(data=world_sf)

b <- st_bbox(world_box)
ggplot() + geom_sf(data=world_sf) +
  coord_sf(crs = crsLONGLAT, xlim = c(b["xmin"], b["xmax"]), ylim = c(b["ymin"], b["ymax"]))

# 投影坐标系下的bbox限制
ggplot()+ geom_sf(data=world_box_robin)  + geom_sf(data=world_robinson) 

b <- st_bbox(world_box_robin)
ggplot() + geom_sf(data=world_robinson) + 
  coord_sf(crs = crsrobin, xlim = c(b["xmin"], b["xmax"]), ylim = c(b["ymin"], b["ymax"]))

解决方案 #

问题出在跨180°经线的多边形被错误拉伸。原始地图数据中，不少国家（如俄罗斯、新西兰）的几何图形跨越180°经线，当直接将中央经线设为-180°进行投影时，这些跨线多边形会被st_transform处理成从投影区域一端直接连接到另一端，从而产生贯穿地图的多余水平线。st_bbox只能限制显示范围，无法修正多边形本身的错误形态。

解决方法是在投影前先处理跨180°的多边形，使用sf包的st_shift_longitude()函数即可自动完成：

# 太平洋中心Robinson投影修正版
crsrobin <- "+proj=robin +lon_0=-180 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"
# 先处理跨180°经线的多边形
world_sf_shifted <- st_shift_longitude(world_sf)
# 再进行投影转换
world_robinson <- st_transform(world_sf_shifted, crs = crsrobin)
# 绘制地图，无多余水平线
ggplot() + geom_sf(data = world_robinson)

st_shift_longitude()的作用是将跨180°经线的多边形拆分为两个独立部分，同时将经度范围转换为0-360°，适配太平洋中心的投影需求，这样投影后就不会出现拉伸线条。

也可以手动转换经度范围来实现相同效果：

# 手动将经度转换为0-360范围
world_sf_360 <- world_sf %>%
  st_transform(crs = "+proj=longlat +datum=WGS84 +no_defs +over") %>%
  mutate(geometry = st_geometry(.) + c(360, 0)) %>%
  st_set_crs(crsLONGLAT)
# 投影到太平洋中心Robinson
world_robinson <- st_transform(world_sf_360, crs = crsrobin)
ggplot() + geom_sf(data = world_robinson)

内容的提问来源于stack exchange，提问作者Mark R