使用 sf 和 ggplot 显示地球的边缘或海洋背景

Question

我正在尝试使用 ggplot 和 sf 获取地球边缘的一条线（或等效的海洋多边形）。我的代码适用于基本的 Mollweide 投影...

library(tmap)
library(sf)
library(s2)
sf_use_s2(TRUE) 
data("World") ## get built-in data from the tmap package
world.sf <- World

lat <- c(seq(from=-90,to=90,by=0.25)) ## creates a sequence -90 to 90 by 0.25
lon <- c(rep(180,length(lat)),rep(-180,length(lat))) # creates a sequence of 180 and -180 longitudes
oceans.df <- data.frame("lon"=lon, "lat"=c(lat,lat)) ## combines the point in a data frame
oceans.sf <- st_as_sf(oceans.df,coords = c("lon","lat"), crs=4326) %>% 
  st_union() %>% # combines the points into a MULTIPOINT
  st_convex_hull() # gets the outline of the points

ggplot() + 
  geom_sf(data=oceans.sf, fill="blue") ## fine for Mercator


oceans_moll.sf <- st_as_sf(oceans.df,coords = c("lon","lat"), crs=4326) %>% 
  st_transform(crs="+proj=moll") %>% ## 
  st_union() %>% # combines the points into a MULTIPOINT
  st_convex_hull() # gets the outline of the points

world_moll.sf <- st_transform(world.sf, crs="+proj=moll")
ggplot() + 
  geom_sf(data=oceans_moll.sf, fill="blue") + 
  geom_sf(data=world_moll.sf, fill="white")

但对于以太平洋为中心的摩尔韦德来说，它失败了

oceans_moll_pacific.sf <- st_as_sf(oceans.df,coords = c("lon","lat"), crs=4326) %>% 
  st_transform(crs="+proj=moll +lon_0=-150 +x_0=0 +y_0=0") %>% ## 
  st_union() %>% # combines the points into a MULTIPOINT
  st_convex_hull() # gets the outline of the points

ggplot() + 
  geom_sf(data=oceans_moll_pacific.sf, fill="blue")

不太确定出了什么问题。我尝试使用下面的代码来定义地球 sf::st_transform() 返回空几何体 但这也不适用于重投影。

ocean <- st_point(x = c(0,0)) %>%
  st_buffer(dist = 6371000) %>%
  st_sfc(crs = "+proj=ortho +lat_0=20 +lon_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")

不确定这里出了什么问题......

Answer 1

基于https://gis.stackexchange.com/questions/468260/oceans-as-global-background-wont-reproject-properly/468292#468292发布的一个很好的答案，这个通用黑客修复了多边形并创建了一个任何本初子午线的海洋/轮廓。

rm(list=ls())
library(rnaturalearthdata)
library(rnaturalearth)
library(sf)
library(tidyverse)

ante_meridian <- -20 ## change as needed
offset <- ante_meridian+180

jitter = 0.01
lat <- c(seq(from=-90,to=90,by=0.25))
lon <- c(rep(offset-jitter,length(lat)),rep(offset+jitter,length(lat)))


oceans.sf <- cbind(lon, c(lat, rev(lat))) %>% st_linestring() %>% st_cast("POLYGON") %>% 
  st_sfc(crs=4326) %>% st_transform(crs=paste0("+proj=moll +lon_0=",ante_meridian," +x_0=0 +y_0=0"))

world_ne.sf <- st_as_sf(ne_countries(scale="medium"))


world_Mollweide_rotated.sf <- world_ne.sf %>%
  st_break_antimeridian(lon_0 = ante_meridian) %>% # insert this before transformation
  # (& ignore the warning)
  st_transform(crs=paste0("+proj=moll +lon_0=",ante_meridian," +x_0=0 +y_0=0"))


ggplot() + geom_sf(data=oceans.sf, fill="blue") + geom_sf(data=world_Mollweide_rotated.sf, 
                                                          fill="grey80")