【问题标题】:How to plot global rasters with tmap in Robinson projection without duplicated areas?如何在 Robinson 投影中使用 tmap 绘制全局栅格而没有重复区域?
【发布时间】:2020-11-04 19:21:54
【问题描述】:

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

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

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

我在 macOS Catalina 10.15 上运行 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。 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")

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

【问题讨论】:

    标签: r raster rgdal tmap r-stars


    【解决方案1】:

    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)
    

    【讨论】:

    • 感谢罗伯特的回答;非常感激!通过掩蔽解决方案效果很好。如果要在某个时候替换 raster,也许我也应该考虑迁移到 terra。