C’est hyper fun mais…

C’est vrai que faire de la 3d dans R est assez grisant et inhabituel. rayshader offre un nombre de fonctions pour arriver à développer son modèle 3d texturé mais les explications fournies sur la page du package sont assez ésotériques de prime abord. J’ai du pas mal chipoter pour arriver à réaliser une image intéressante. Le plus fastidieux à été de trouver une source pour les images d’élévation des terrains.

Le package rayshader

La première étape est de prendre connaissance avec le package rayshader. La version au moment de l’écriture de ce post est la version 0.10.1. Puis de charger le package dans R. Vous aurez aussi besoin du package raster pour transformer votre image en une matrice.

library(rayshader)
library(raster)

## Loading required package: sp

## 
## Attaching package: 'raster'

## The following object is masked from 'package:dplyr':
## 
##     select

## The following object is masked from 'package:tidyr':
## 
##     extract

library(png)

1.1. Appliquer la fonction raster()

Nous allons “rasterizer” cette image.

elevation <- "data/rayshader_tuto/assets/tuto_base/square_elevation.tif"
elevation.raster <- raster::raster(elevation)

La fonction raster crée un objet raster contenant les informations nécessaires pour la transformation en matrix.

1.2. Transformer en matrix.

Cette étape transforme l’image d’élévation en matrix. Vous pouvez observer l’inversion au niveau nrow() et ncol(). La fonction extract extrait les valeurs de l’objet, quant à la fonction extent elle définit la taille de la matrice.

elevation.matrix <- matrix(extract(elevation.raster, extent(elevation.raster)), nrow = ncol(elevation.raster), ncol = nrow(elevation.raster))

Voyons élément par élément l’output.

extent(elevation.raster)

## class      : Extent 
## xmin       : 0 
## xmax       : 600 
## ymin       : 0 
## ymax       : 600

ncol(elevation.raster)

## [1] 600

nrow(elevation.raster)

## [1] 600

Affichons la dimension de notre variable elevation.matrix.

dim(elevation.matrix)

## [1] 600 600

Nous avons une matrix avec 600 x 600 valeurs soit 360,000 valeurs.

1.3. Calculer la “surface color map”.

À partir de notre matrix, nous allons calculer une couleur pour chaque point. La couleur est calculée selon la texture par défault incluse dans le package et l’orientation du point. Il y 6 textures : imhof1,imhof2,imhof3,imhof4,desert, bw et unicorn. Vous pouvez même faire votre propre texture avec la fonction create_texture().

À quoi cela sert-il? Et bien à simuler la lumière sur l’objet 3d qui sera créé après. La lumière dans ce cas-ci est paramétrée à 35°. La texture choisie est desert et la relation x, y avec z est de 0.01.

elevation.matrix %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = .05) %>%
  plot_map()

1.4. Calculer l’“ambiant occulsion mapping”.

L’ambiant occulusion va assombrir les zones difficilement accessibles à la lumière. Ce qui va renforcer le réalisme. Dans ce cas-ci ce n’est pas très parlant, mais sur une élévation sillonnée c’est plus intéressant.

elevation.matrix %>% 
  ambient_shade() %>%
  plot_map()

1.5. Calculer les ombres portées en raytracing.

La 4e étape consiste à projeter les ombres par raytracing. Le paramètre anglebreaks va positionner la source de lumière en degré à partir de l’horizon. Essayons quelques valeurs.

sunangle = 0

elevation.matrix %>% 
  ray_shade(anglebreaks = 0, sunangle = 35, zscale = 0.05) %>%
  plot_map()

sunangle = 45

elevation.matrix %>% 
  ray_shade(anglebreaks = 45, sunangle = 35, zscale = 0.05) %>%
  plot_map()

sunangle = 90

elevation.matrix %>% 
  ray_shade(anglebreaks = 90, sunangle = 35, zscale = 0.05) %>%
  plot_map()

Pour obtenir une ombre moins dure, le paramètre autorise une séquence de valeur. Partons de 0 à 45° par incrément de 5.

elevation.matrix %>% 
  ray_shade(anglebreaks = seq(0,45, 5), sunangle = 35, zscale = 0.05) %>%
  plot_map()

Pour obtenir une ombre moins dure, le paramètre autorise une séquence de valeur. Partons de 0 à 20° par incrément de 1.

elevation.matrix %>% 
  ray_shade(anglebreaks = seq(0,20, 1), sunangle = 35, zscale = 0.05) %>%
  plot_map()

1.5. Le rendu 3d.

Nous voici à la dernière étape. Assemblons l’ensemble pour constituer notre modèle 3d. Dans l’ordre :

• sphere_shade
  
• ambiant_shade
  
• ray_shade

Pour ajouter les ombres (comme pour une texture, nous verrons cela après), il faut le faire avec la fonction add_shadow. Son paramètre max_darken permet de régler l’opacité du calque.

Enfin, nous appliquons la fonction plot_3d qui imprime le modèle 3d à l’écran. theta ajuste la rotation de la camera et phi l’angle à partir de l’azimut.

render_snapshot() affiche un snapshot du rendu.

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = 0.01) %>%
  add_shadow(ambient_shade(elevation.matrix)) %>% 
  add_shadow(ray_shade(elevation.matrix, anglebreaks = seq(0,20, 1), sunangle = 35, zscale = 0.01), max_darken = 0) %>%
  plot_3d(heightmap = elevation.matrix , 
          zscale = 0.01, 
          fov = 90,
          lineantialias = TRUE,
          theta = 45,
          phi = 15,
          zoom = 0.7)

render_snapshot()

rgl::rgl.clear()

1.6. Ajout d’une texture.

Le package permet aussi d’ajouter une texture, comprenez par là une image sur le modèle 3d. La procédure est assez similaire que pour transformer notre élévation en matrix. Voici l’image que nous allons utiliser.

Nous utilisons la fonction readPNG du package png.

elevation.texture.map <- readPNG("data/rayshader_tuto/assets/tuto_base/square_texture.png")

Assemblons à nouveau le tout avec la fonction add_overlay pour ajouter la texture. Il est possible de rendre transparente une couleur particulière. Dans ce cas, nous n’en avons pas besoin. Réglons l’opacité à .99 avec alphalayer.

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = 0.01) %>%
  add_overlay(elevation.texture.map, alphacolor = NULL, alphalayer = 0.99) %>% 
  add_shadow(ambient_shade(elevation.matrix)) %>% 
  add_shadow(ray_shade(elevation.matrix, anglebreaks = seq(0,20, 1), sunangle = 35, zscale = 0.01), max_darken = 0.1) %>%
  plot_3d(heightmap = elevation.matrix , 
          zscale = 0.01, 
          fov = 90,
          lineantialias = TRUE,
          theta = 45,
          phi = 15,
          zoom = 0.7)

render_snapshot()

rgl::rgl.clear()

BOUM ! Voici notre rendu complet.

Pour régler la hauteur de l’élévation, nous pouvons enregistrer la valeur de l’échelle z dans un vecteur.

my.z <- 0.1

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = my.z) %>%
  add_overlay(elevation.texture.map, alphacolor = NULL, alphalayer = 0.99) %>% 
  add_shadow(ambient_shade(elevation.matrix)) %>% 
  add_shadow(ray_shade(elevation.matrix, anglebreaks = seq(0,20, 1), sunangle = 35, zscale = my.z), max_darken = 0.1) %>%
  plot_3d(heightmap = elevation.matrix , 
          zscale = my.z, 
          fov = 90,
          lineantialias = TRUE,
          theta = 45,
          phi = 15,
          zoom = 0.7)

render_snapshot()

rgl::rgl.clear()

my.z <- 0.005

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = my.z) %>%
  add_overlay(elevation.texture.map, alphacolor = NULL, alphalayer = 0.99) %>% 
  add_shadow(ambient_shade(elevation.matrix)) %>% 
  add_shadow(ray_shade(elevation.matrix, anglebreaks = seq(0,20, 1), sunangle = 35, zscale = my.z), max_darken = 0.1) %>%
  plot_3d(heightmap = elevation.matrix , 
          zscale = my.z, 
          fov = 90,
          lineantialias = TRUE,
          theta = 45,
          phi = 15,
          zoom = 0.7)

render_snapshot()

rgl::rgl.clear()

2.1. Choisir sa zone.

Pour choisir sa zone, allez dans la rubrique data / find data. Une carte du monde s’affiche. Les points rouges et verts présentent des datasets Hi-Res disponibles. Ces datasets contiennent des informations plus denses. Néanmoins vous pouvez sélectionner la zone que vous souhaitez.

2.2. Hawai - Hilo

A l’aide de la souris, sélectionner une zone dans les limites affichées sur la carte (en rouge).

Vérifier le nombre de points inclus dans votre sélection pour ne pas avoir un fichier trop volumineux.

Les options à sélectionner :

• 1. Select Data Output Format : GeoTiff
• 2. Select Data Output Resolotion : Max

Vous pouvez également cochez la génération hillsahde et color relief pour avoir une texture à ajouter. Compléter cette texture avec une image de google map que vous pouvez positionner par dessus.

Inscrivez un titre en bas de page et un nom de projet.

Le temps de calcul peut être un peu long. Vous recevrez un email pour vous avertir quand le fichier est prêt. Prenez les éléments suivants :

• Job Results
• Visualization Products

2.3. Modéliser

Une fois les fichiers souhaités téléchargés, nous pouvons procéder à la modélisation.

Raster :

elevation <- "data/rayshader_tuto/assets/hawai/hilo.tif"
elevation.raster <- raster::raster(elevation)

Matrix :

elevation.matrix <- matrix(extract(elevation.raster, extent(elevation.raster), buffer = 100), nrow = ncol(elevation.raster), ncol = nrow(elevation.raster))

zscale :

my.z <- 30

Texture :

elevation.texture.map <- readPNG("data/rayshader_tuto/assets/hawai/hilo_color_text.png")

Test map

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "imhof1", zscale = my.z) %>%
  plot_map()

Ambient occlusion :

elevation.amb.shade <- ambient_shade(elevation.matrix, zscale = my.z)

Ray Shadow :

elevation.ray.shade <- ray_shade(elevation.matrix,  sunangle = 35, zscale = my.z, )

Ajout de la texture :

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = my.z) %>%
  add_overlay(elevation.texture.map, alphacolor = NULL, alphalayer = 0.9, ,gamma_correction = TRUE) %>% 
  plot_map()

Rendu 3d :

elevation.matrix  %>% 
  sphere_shade(sunangle = 35, texture = "desert", zscale = my.z) %>%
  add_overlay(elevation.texture.map, alphacolor = NULL, alphalayer = 0.9) %>% 
  add_shadow(elevation.amb.shade) %>% 
  add_shadow(elevation.ray.shade, 0.7) %>%
  plot_3d(heightmap = elevation.matrix, 
          zscale = my.z, 
          fov = 90,
          lineantialias = TRUE,
          theta = 45,
          phi = 15,
          zoom = 0.7)

render_water(elevation.matrix, zscale = my.z, wateralpha = 0.3, waterlinecolor = "white", watercolor = "turquoise4")
render_camera(theta = 135, phi = 25, zoom = 0.8, fov = 75)
render_label(heightmap = elevation.matrix, x = 797, y = 963-509, z = 10000, linewidth = 1, zscale = my.z, text = "Hilo")
render_label(heightmap = elevation.matrix, x = 90, y = 963-304, z = 10000, linewidth = 1, zscale = my.z, text = "Mauna Kea - 4205m")
render_snapshot("hawai_render.png")
render_snapshot()

rgl::rgl.clear()

Le final :

Je récupère le fichier dans photoshop pour quelques retouches.

A vous de jouer!
Si vous tweetez vos modèles, n’hésitez pas à me mettre en copie.