时间戳

首页 > 解决方案 > 如何在 Robinson 投影中使用 tmap 绘制全局栅格而没有重复区域？

r - 如何在 Robinson 投影中使用 tmap 绘制全局栅格而没有重复区域？

问题描述

我最近主要使用栅格和 tmap 绘制了一些全球栅格。我想用罗宾逊投影而不是经纬度来绘制地图。然而，从下图（阿拉斯加、西伯利亚、新西兰）可以看出，对 Robinson 的简单投影会复制地图边缘的某些区域。

以前，我发现了一种使用 PROJ.4 代码参数“+over”的解决方法，如here和here中所述。

随着使用 GDAL > 3 和 PROJ >= 6 对 rgdal 的最新更改，此解决方法似乎已过时。有没有人找到一种新方法来绘制 Robinson/Eckert IV/Mollweide 中的全球栅格而没有重复区域？

我在 macOS Catalina 10.15.4 上运行 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

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")

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

标签： rrasterrgdaltmapr-stars

解决方案

有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，可从 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)

推荐阅读