A quick demo on how to use rayshader to create 3D thematic views of the Earth, e.g. here a visualization of gross primary production at 0.5 degree resolution

[related to my twitter post https://twitter.com/Reichstein_BGC/status/1330583433035337730]

Loading libraries and define functions

options(rgl.useNULL=TRUE) ## needed with knitR cf. https://community.rstudio.com/t/rgl-scene-not-showing-in-knitted-html-document-under-rstudio-1-2-5001/43582/4
library(rayshader)
library(raster)
library(rgl)
library(scales)
library(tidyverse)

values2rgb  <- function(values, pal=pals::viridis(64)) {
  
  
  ncol <- length(pal)
  cols <- pal[scales::rescale(values, c(1,ncol))] %>% col2rgb(alpha=TRUE)
  cols  <- cols %>% array(dim=c(4, values %>% dim)) %>% aperm(c(3,2,1))
  cols/255.
}

## Function to rayshade an elevation matrix with specified coloring values, which will be mapped onto a color scale with values2rgb
rs_surface  <- function(elmat, coloring, water = FALSE, alpha_over=0.75, 
                        texture="bw", zscale=250, shadowintens=0.8, shadowReachFactor=1, windowsize=c(1920, 1080), zoom=0.5, ...)
{
  elmat %>% 
    
    sphere_shade(texture=texture,  colorintensity = 1, zscale = zscale) %>% 
    #ambient_shade() %>%
    add_shadow(ray_shade(elmat,zscale=zscale/shadowReachFactor),max_darken = 1-shadowintens) %>% 
    add_overlay( coloring %>% values2rgb, alphalayer = alpha_over) %T>%
    plot_map() %>%
    plot_3d(zscale = zscale, heightmap = elmat, wateralpha = 0.7, water=water, windowsize = windowsize, zoom = zoom, ...)
  
  
  }

Preparing data

## Elevation matrix from NOAA http://www.ngdc.noaa.gov/mgg/global/global.html (downsampled by U. Weber)
  elmat <- raster("M:/data/DataStructureMDI/DATA/grid/Global/0d50_static/ETOPO/ETOPO1/Data/ETOPO1.halfdegree.nc") %>%  raster_to_matrix() 
## Thematic values: here, Gross Primary Production from Beer et al. 2010, Science
  gpp  <-  raster("D:/markusr/_publi/onTheWay/_FLUXNET_LaThuile/_Beer_GPP/data/GPPdata_Beer_etal_2010_Science.nc") %>%  raster_to_matrix()

Plotting

 ### Now plot map of gpp where gpp also gives the elevation...
 ### Some tweaking to make plot nicer (not these steep grey cliffs)
 ### !!! for the filter needs to install imagine package
  gppFilt <- imagine:: quantileFilter(gpp, 3,0.001)
  gpp2Plot <- gpp
  gpp2Plot[is.na(gpp)]  <- gppFilt[is.na(gpp)] * 0.0
  
  rgl.clear(); rs_surface(gpp2Plot, zscale = 50, gpp, water=FALSE, shadowReachFactor = 2, shadowintens = 0.8, 
                          zoom=0.7, solid=T, soliddepth=-10, theta=0, windowsize = c(1920, 1080))

Next is the WebGL generated with rglwidget().

  rglwidget()
  ### Now plot map of gpp where gpp als gives the elevation but inlcude ocean...
  ## Missing values kind of a problem still
  elmat4plot <- ifelse(elmat < 0, elmat/5, gpp2Plot )
  rgl.clear(); rs_surface(elmat4plot, zscale = 100, gpp, water=TRUE, shadowReachFactor = 2, shadowintens = 0.8, theta=0)

  #rglwidget()
  
  #use rayshader::render_movie() to render a movie, I used frames=720

Next would be the WebGL generated with rglwidget() (too big to upload)