Search code examples
rrasterrgdaltmapr-stars

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

I've been plotting some global rasters lately using mainly raster and tmap. I'd like to plot the maps in Robinson projection instead of lat-lon. Simple projection to Robinson however duplicates some areas on the edges of the map as you can see from the figures below (Alaska, Siberia, NZ).

Previously, I found a workaround with PROJ.4 code parameter "+over" as outlined in here and here.

With the latest changes to rgdal using GDAL > 3 and PROJ >= 6, this workaround seems to be obsolete. Has anyone found a new way on how to plot global rasters in Robinson/Eckert IV/Mollweide without duplicated areas?

I'm running 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 and PROJ 6.3.1 on a 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

enter image description here

I tried to do the same with stars instead of raster but no resolution was found, supposedly since tmap uses stars since version 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")

Thanks for any insights - hoping someone else is thinking about this issue!

Solution

  • With raster you can do

    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)

    Create a mask

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

    And use it

    mrob <- mask(rrob, re)

    The new package terra has a mask argument for that (you need version >= 0.8.3 for this, 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)