如何在没有重复区域的情况下在 Robinson 投影中使用 tmap 绘制全局栅格?

How to plot global rasters with tmap in Robinson projection without duplicated areas?

我最近一直在绘制一些全球栅格,主要使用栅格和 tmap。我想用罗宾逊投影而不是经纬度绘制地图。然而,简单投影到罗宾逊会复制地图边缘的某些区域,如下图所示(阿拉斯加、西伯利亚、新西兰)。

之前,我找到了一个使用 PROJ.4 代码参数“+over”的解决方法,如 here and here 中所述。

随着使用 GDAL > 3 和 PROJ >= 6 对 rgdal 的最新更改,此解决方法似乎已过时。有没有人找到如何在 Robinson/Eckert IV/Mollweide 中绘制没有重复区域的全局栅格的新方法?

我是 运行 R 4.0.1,tmap 3.1,stars 0.4-3,raster 3.3-7,rgdal 1.5-12,sp 1.4-2,GDAL 3.1.1 和 PROJ 6.3.1在 macOS Catalina 10.15.4

require(stars)
require(raster)
require(tmap)
require(dplyr)

# data
worldclim_prec = getData(name = "worldclim", var = "prec", res = 10)
jan_prec <- worldclim_prec$prec1

# to Robinson and plot - projection outputs a warning
jp_rob <- jan_prec %>%
  projectRaster(crs = "+proj=robin +over")
tm_shape(jp_rob) + tm_raster(style = "fisher")
Warning messages:
1: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded ellps WGS 84 in CRS definition: +proj=robin +over
2: In showSRID(uprojargs, format = "PROJ", multiline = "NO") :
  Discarded datum WGS_1984 in CRS definition

我尝试用星星而不是光栅做同样的事情,但没有找到分辨率,据说是因为 tmap 从 3.0 版开始使用星星。

# new grid for warping stars objects
newgrid <- st_as_stars(jan_prec) %>%
  st_transform("+proj=robin +over") %>%
  st_bbox() %>%
  st_as_stars()

# to stars object - projection outputs no warning
jp_rob_stars <- st_as_stars(jan_prec) %>%
  st_warp(newgrid)

tm_shape(jp_rob_stars) + tm_raster(style = "fisher")

感谢您提供任何见解 - 希望其他人正在考虑这个问题!

有了raster你可以做到

library(raster)
prec <- getData(name = "worldclim", var = "prec", res = 10)[[1]]
crs <- "+proj=robin +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m"
rrob <- projectRaster(prec, crs=crs)

创建遮罩

library(geosphere)
e <- as(extent(prec), "SpatialPolygons")
crs(e) <- crs(prec)
e <- makePoly(e)  # add additional vertices
re <- spTransform(e, crs)

并使用它

mrob <- mask(rrob, re)

新包 terra 有一个 mask 参数(你需要版本 >= 0.8.3,available from github

prec <- getData(name = "worldclim", var = "prec", res = 10)[[1]]
jp <- rast(prec$prec1) 
jp <- jp * 1 # to deal with NAs in this datasaet
rob <- project(jp, crs, mask=TRUE)