An example of the robinson maps often used to represent global modelling data

robinson_map.r

library(raster)
library(maps)
library(maptools)
library(sf)
library(scales)

# set coordinate systems
robinson = CRS(" +proj=robin +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs")
latlon = CRS("+init=epsg:4326")

# set the colours for the map
# red and blues in this case for temperatures
cols = dichromat_pal("DarkRedtoBlue.12")(12)
cols = gradient_n_pal(cols)(seq(0, 1, length = 15))

# load data
data("wrld_simpl")

# download a file here: https://data.giss.nasa.gov/gistemp/maps/
temp = raster("amaps.nc") # I downloaded the map for April 2017

# create grid lines
pts=SpatialPoints(rbind(c(-180,-90),
                        c(-180,90),
                        c(180,90),
                        c(180,-90)),
                  CRS("+init=epsg:4326"))
gl = gridlines(pts, easts = seq(-180,180,180), norths = seq(-90,90,90), ndiscr = 100)
bb = gridlines(pts, easts = seq(-180,180,360), norths = seq(-90,90,180), ndiscr = 100)

# reproject grid lines and bounding box
gl_robinson = spTransform(gl, robinson)
bb_robinson = spTransform(bb, robinson)

# transform the wold map and data
wrld_robinson = spTransform(wrld_simpl, robinson)

# reproject the gridded data
temp_robinson = trim(projectRaster(temp, crs = robinson, over = TRUE))
temp_robinson[temp_robinson>7] = 7 # zlim otherwise makes these missing values

par(mar=c(5,0,3,0))
plot(temp_robinson,
     xaxt = "n",
     yaxt = "n",
     bty = "n",
     zlim = c(-7,7),
     col = cols,
     box = FALSE,
     horiz = TRUE,
     legend.args=list(text=expression('Temperature Anomaly (' * degree * 'C)'),
                      side = 3,
                      font = 2,
                      line = 2,
                      cex = 1.5))

lines(wrld_robinson, col = "grey25")
lines(gl_robinson, col = "grey25")
lines(bb_robinson, lwd = 2)

mtext(text = "April 2017 - Anomaly vs 1951-1980",
      line = 1,
      side = 3,
      cex = 1.5)