Visualizing Water: Interactive 3D Map of River Basins (draft)

Part 3: how to represent water flux over 3D maps of river basins using WebGL.

This notebook is Part 3 of an exploration to visualize results of hydrologic models. In Part 1 we built custom HTML widgets using D3.js, and in Part 2 we looked at rendering water fluxes using Sankey diagrams. Here we test multiple libraries to generate hillshade (3D) views of river basins and water infrastructure, in particular we want to compare Three.js and WebGL implementations.

Aside from rendering topography and water streams in 3D (and potentially other covariate layers), our objective is to overlay custom labels to illustrate the water cycle.

Another objective is to provide basin and sub-basin statistics, starting with precipitation, ET, soil moisture, and adding other covariates, such as population, land use allocation, and crop allocation.

Some inspiration below:

• Tyler Morgan Tutorial: Adding OSM Data to Rayshader Maps
• Flood depth on 3D terrain
• Procedural-gl examples
• Cesium platform for 3D Geospatial
• NASA Web WorldWind
• WebGL with Google Maps API on how to construct and render 3D scenes with WebGL.

Show code

list.files("./fig", full.names=TRUE)[c(2,7)] %>%
  knitr::include_graphics()

WebGL: Sample Scene Processing with Rayshader¹

Data Acquisition

We’ll experiment with a smaller sample scene of the Selingue Dam in the Niger River basin.

basin <- shapefile("~/Projects/WADashboard/shared/mli/srtm/mli_basin.shp")
zoi <- shapefile("~/Projects/WADashboard/shared/mli/srtm/zoi.shp")
ext <- extent(zoi)
center <- coordinates(zoi)

# Get CGIAR SRTM DEM at 90m
srtm <- getData("SRTM", lon=center[,1], lat=center[,2], path="_data") %>%
  crop(zoi) %>%
  mask(zoi)
srtm

class      : RasterLayer 
dimensions : 838, 973, 815374  (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333  (x, y)
extent     : -8.6725, -7.861667, 11.06917, 11.7675  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : srtm_35_10 
values     : 326, 476  (min, max)

# Satellite basemaps
bmap <- maptiles::get_tiles(terra::ext(zoi), "Esri.WorldImagery", zoom=10) %>%
  stack() %>%
  crop(zoi) %>%
  mask(zoi)

The basin covers a large area, so we would need 8 SRTM tiles, but 1 is enough for a proof of concept. Next we’ll get a satellite basemap.

plot(terra::vect(basin), col=pal[2], border=pal[1], lwd=2,
  main="Niger River Basin (Mali)")
plot(zoi, lty=3, col=alpha(pal["red"], .6), border=pal["red"], lwd=2, add=T)
text(-8, 10.5, "Selingue Dam\n(Mali)", col=pal["red"], cex=.7, font=2)
grid()

plot(terra::ext(zoi), 
  main="Selingue Dam (Mali) - ESRI World Imagery")
plotRGB(bmap, add=T)
grid(col="white")

plot(terra::rast(srtm),
  main="Selingue Dam (Mali) - SRTM 90m")
grid(col="white")

Scene Rendering

Next we convert the 2 rasters to a matrix format that’s compatible with Rayshader hillshading and raytracing algorithms.

Show code

# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(srtm)

# Convert sat basemap to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)

bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255

bmap_array %>%
  aperm(c(2,1,3)) %>%
  # Stretch contrast
  rescale(to=c(0,1)) -> bmap_array

2D view of the generated basemap with satellite image overlay.

srtm_water <- srtm_array
srtm_water[srtm_water < 353] <- 0

base_map_sat <- srtm_array %>%
  height_shade() %>%
  add_overlay(bmap_array) %>%
  add_shadow(ray_shade(srtm_array, zscale=90)) %>%
  add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))

plot_map(base_map_sat)

Figure 1: Hillshaded Basemap of Selingue Dam (Niger River basin)

That doesn’t look very clear, so instead we’ll create a basemap, not using the satellite image but a built-in texture.

base_map <- srtm_array %>% 
  height_shade() %>% 
  add_overlay(sphere_shade(srtm_array, texture="desert", 
    zscale=4, colorintensity=5), alphalayer=0.5) %>%
  add_shadow(lamb_shade(srtm_array, zscale=6), 0) %>%
  add_shadow(ambient_shade(srtm_array), 0) %>%
  add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0.1) %>%
  add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))

saveRDS(base_map, "./_data/base_map.rds")

Show code

base_map <- readRDS("./_data/base_map.rds")
plot_map(base_map)

Figure 2: Hillshaded Basemap of Selingue Dam

Looks better, so let’s acquire and overlay spatial features from OSM. Note that water-related features for Mali are scarcely available from OSM, but most waterways are recorded.

osm_roads <- opq(bbox(zoi)) %>% 
  add_osm_feature("highway") %>% 
  osmdata_sf()

osm_water = opq(bbox(zoi)) %>% 
  add_osm_feature("water") %>% 
  osmdata_sf()

osm_waterway = opq(bbox(zoi)) %>% 
  add_osm_feature("waterway") %>% 
  osmdata_sf()

osm_place = opq(bbox(zoi)) %>% 
  add_osm_feature("place") %>% 
  osmdata_sf()

road_layer <- generate_line_overlay(
  dplyr::filter(osm_roads$osm_lines, highway %in% c("primary", "secondary")),
  extent=ext, srtm_array, linewidth=5, color=pal["black"])

water_layer <- generate_line_overlay(
  osm_waterway$osm_lines, 
  extent=ext, srtm_array, linewidth=3, color=pal["blue"])

place_layer <- generate_label_overlay(
  dplyr::filter(osm_place$osm_points, !is.na(name) & nchar(name)<10), 
  extent=ext, heightmap=srtm_array,
  font=2, text_size=1.6, point_size=1.6, color=pal["black"],
  halo_color="white", halo_expand=2, halo_blur=1, halo_alpha=.9, seed=1,
  data_label_column="name")

scene <- base_map %>% 
#scene <- base_map_sat %>%   
  add_overlay(road_layer) %>%
  add_overlay(water_layer, alphalayer=1) %>%
  add_overlay(place_layer)

saveRDS(scene, "./_data/scene.rds")

Finally we’ll use WebGL to render this scene in 3D. We also test adding polygon (bar) annotations.

amb_layer <- ambient_shade(srtm_array, zscale=1/5)

scene2 <- srtm_array %>% 
  height_shade() %>%
  add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0) %>%
  add_shadow(amb_layer, 0) %>%
  add_overlay(road_layer) %>%
  add_overlay(water_layer, alphalayer=1) %>%
  add_overlay(place_layer)

saveRDS(scene2, "./_data/scene2.rds")

Show code

scene <- readRDS("./_data/scene.rds")
scene2 <- readRDS("./_data/scene2.rds")

scene %>% plot_3d(resize_matrix(srtm_array, 1/10), zscale=20, 
  theta=30, phi=20, fov=0, zoom=0.5,
  family="Roboto Condensed",
  shadow=TRUE, shadowcolor=pal["black"], solidcolor=pal["black"])

# Add polygon annotations
xyz <- sf::read_sf("~/Projects/WADashboard/shared/mli/srtm/xyz.shp")
render_polygons(xyz, ext, data_column_top="z",
  scale_data=1, color=alpha(pal["orange"], 0.8), 
  lit=F, light_intensity=0.01, clear_previous=T)

rglwidget()

Figure 3: Interactive 3D Scene of Selingue Dam

The generated scene (PNG image and associated JS code) is 3.5MB.

WebGL: Full Scene with Rayshader

Following the approach above, we experiment with a 3D scene for the entire Niger River basin area used in the WA+ analysis.

Data Acquisition

# Basin boundaries from IWMI
zoi <- shapefile("~/Projects/WADashboard/shared/mli/srtm/mli_basin.shp")
ext <- extent(zoi)
center <- coordinates(zoi)

# GAUL admin-2 boundaries
adm <- sf::st_read(
  "~/Projects/WADashboard/shared/mli/srtm/mli_adm2_lines.shp")
places <- sf::st_read(
  "~/Projects/WADashboard/shared/mli/srtm/mli_adm2_centroids.shp")

# CGIAR SRTM DEM (upsampled)
alt <- lapply(c("MLI", "GIN", "CIV"), 
  function(x) getData("alt", country=x, mask=FALSE) %>%
    crop(zoi) %>%
    mask(zoi)
)
alt <- mosaic(alt[[1]], alt[[2]], alt[[3]], fun=mean)
alt <- mask(crop(alt, zoi), zoi)
writeRaster(alt, "~/Projects/WADashboard/shared/mli/srtm/mli_alt.tif",
  overwrite=T)
alt

# ESA WorldCover (clipped in QGIS)
luc <- raster("~/Projects/WADashboard/shared/mli/srtm/ESA WorldCover clip.tiff")

# Note that I can't find a colormap for this raster, so switching to ESA CCI instead
luc <- raster("~/Projects/WADashboard/shared/mli/srtm/mli_esa_cci_300m.tif")

# And its colormap
pal.luc <- fread("~/Maps/ESA/ESA CCI Colormap.txt")
setnames(pal.luc, c("value", "R", "G", "B", "A", "label"))
pal.luc[, hex := rgb(R, G, B, maxColorValue=255)]

# We'll use GloRIC v1.0 stream network instead of OSM
# Features are filtered to 'Class_hydr > 12'
gloric <- sf::st_read("~/Projects/WADashboard/shared/mli/srtm/mli_gloric_filtered.shp")

# Country labels
osm_place <- opq(bbox(zoi)) %>% 
  add_osm_feature("place", "country") %>% 
  osmdata_sf()

save(zoi, ext, adm, places, alt, luc, pal.luc, gloric, osm_place,
  file="./_data/osm.RData")

load("./_data/osm.RData")

plot(terra::rast(alt),
  main="Niger River Basin (Mali) - SRTM 1km")
plot(adm, lty=3, col=pal["black"], lwd=1, add=T)
plot(gloric, col=pal["blue"], lwd=.4, add=T)
grid()

plot(terra::rast(luc), col=pal.luc$hex, 
  breaks=pal.luc$value, legend=F,
  main="Niger River Basin (Mali) - ESA CCI 300m")
plot(adm, lty=3, col=pal["black"], lwd=1, add=T)
grid()

Scene Rendering

Same configuration as above.

Show code

# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(alt)

bmap <- RGB(luc, col=pal.luc$hex, breaks=c(0,pal.luc$value))

# Convert land cover raster to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)

bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255

bmap_array <- bmap_array %>%
  aperm(c(2,1,3)) %>%
  # Stretch contrast
  rescale(to=c(0,1))

2D view of the generated basemap with land cover overlay.

base_map_luc <- srtm_array %>%
  height_shade() %>%
  add_overlay(bmap_array) %>%
  add_shadow(lamb_shade(srtm_array, zscale=6), 0.5) %>%
  add_shadow(ambient_shade(srtm_array), 0.5)

saveRDS(base_map_luc, "./_data/base_map_luc_basin.rds")

Show code

base_map_luc <- readRDS("./_data/base_map_luc_basin.rds")
plot_map(base_map_luc)

Figure 4: Hillshaded Land Cover - Niger River Basin

The choice of hillshading parameters below might need tweaking.

base_map <- srtm_array %>%
  sphere_shade(texture="desert") %>%
  #add_shadow(ray_shade(srtm_array, zscale=50))
  add_shadow(lamb_shade(srtm_array, zscale=6), 0.5) %>%
  # This step takes up to 5min on 4 cores / 8GB
  add_shadow(ambient_shade(srtm_array), 0.5)

saveRDS(base_map, "./_data/base_map_basin.rds")

2D view of the generated basemap with simple texture overlay.

Show code

base_map <- readRDS("./_data/base_map_basin.rds")
plot_map(base_map)

Figure 5: Hillshaded Basemap - Niger River Basin

As above, we choose to overlay place names, waterways, and level-2 administrative boundaries.

admin_layer <- generate_line_overlay(adm,
  extent=ext, srtm_array, linewidth=4, lty=1, color="white")

admin_layer2 <- generate_line_overlay(adm,
  extent=ext, srtm_array, linewidth=2, lty=1, color=pal["black"], 
  offset=c(0, 0.03))

water_layer <- generate_line_overlay(gloric, 
  extent=ext, srtm_array, linewidth=6, color=pal["light-blue"],
  data_column_width="Class_hydr")

place_layer <- generate_label_overlay(places,
  extent=ext, heightmap=srtm_array,
  font=2, text_size=1.6, point_size=0, color="white",
  halo_color=pal["black"], halo_expand=2, halo_blur=2, halo_alpha=.6,
  seed=1, data_label_column="ADM2_NAME")

country_layer <- generate_label_overlay(osm_place$osm_points,
  extent=ext, heightmap=srtm_array, offset=c(0.1, -0.5),
  font=3, text_size=3, point_size=0, color=pal["light"],
  halo_color=pal["black"], halo_expand=4, halo_blur=4, halo_alpha=.6,
  halo_offset=c(0.1, -0.5),
  seed=1, data_label_column="name")

scene <- base_map %>%
  add_overlay(admin_layer) %>%
  add_overlay(water_layer) %>%
  add_overlay(place_layer)

plot_map(scene)
saveRDS(scene, "./_data/scene_basin.rds")

scene_luc <- base_map_luc %>% 
  add_overlay(admin_layer2, alphalayer=.8) %>%
  add_overlay(admin_layer, alphalayer=.9) %>%
  add_overlay(country_layer) %>%
  add_overlay(water_layer) %>%
  add_overlay(place_layer)

plot_map(scene_luc)
saveRDS(scene_luc, "./_data/scene_luc_basin.rds")

Finally we’ll use WebGL to render these scenes in 3D.

Show code

scene <- readRDS("./_data/scene_basin.rds")
scene %>% plot_3d(resize_matrix(srtm_array, 1/5), zscale=40, 
  linewidth=0, theta=15, phi=45, fov=0, zoom=0.4, 
  family="Roboto Condensed",
  shadow=TRUE, shadowcolor=pal["black"], solidcolor=pal["black"])

rglwidget()

Figure 6: Interactive 3D Scene of Niger River Basin

Show code

rgl.close()

Show code

scene_luc <- readRDS("./_data/scene_luc_basin.rds")
scene_luc %>% plot_3d(resize_matrix(srtm_array, 1/5), zscale=40, 
  linewidth=0, theta=15, phi=45, fov=0, zoom=0.4, 
  family="Roboto Condensed",
  shadow=TRUE, shadowcolor=pal["black"], solidcolor=pal["black"])

rglwidget()

Figure 7: Interactive 3D Scene of Niger River Basin with Land Cover Basemap

Show code

rgl.close()

#save_3dprint("mli_basin.stl", maxwidth="120mm", clear=TRUE)
#save_obj("mli_basin.obj")

3D Processing with Three.js

For comparison’s sake we build a similar scene using Three.js libraries instead of WebGL. The scene is 12 MB in size (using a low resampling rate), it was generated using the Qgis2threejs plugin. The main difference is we’re showing ESA CCI land cover classification as an image overlay instead of a texture.

The widget is rendered here.

In Part 4 we experiment with another WebGL-derived Javascript library Deck.gl and 3D rendering of terrain layers.

---

1. https://www.rayshader.com/↩︎