Skip to content

Examples

Jump to bottom
Michael Sumner edited this page Jul 9, 2020 · 24 revisions

Bake polygons into a raster mesh

  • in real world coords
  • raster is a triangle mesh (but could be quads)
  • lines added by first densifying to the DEM resolution, then copy down the polygon boundaries

library(tigris)
library(raster)
nv <- counties(state = "Nevada", cb = TRUE)
library(ceramic)
elv <- cc_elevation(nv, max_tiles = 36)

p <- sf::st_transform(sf::st_as_sf(nv), raster::projection(elv))
l <- sf::st_segmentize(p, mean(res(elv)))
p$id <- 1:nrow(p)
img <- fasterize::fasterize(p, elv, "id")
## where it's NA, make it blue
blue <- cellStats(img, max) + 1
img[is.na(img)] <- blue
cols <- sample(viridis::viridis(nrow(p)))
cols <- c(cols, "#FFC0CBFF")
cp <- col2rgb(cols)

## create rgb raster
rrggbb <- setValues(brick(img, img, img), t(cp[, values(img) ]))
library(anglr)
mesh <- as.mesh3d(elv, image_texture = rrggbb)
mesh$material$color <- "#FFFFFFFF"
plot3d(mesh); rgl::aspect3d(1, 1, .02)
plot3d(copy_down(SC0(l), elv + 10), add = TRUE, col = "grey", lwd = 4)

Walls

https://twitter.com/mdsumner/status/1257742410190225408/photo/1

library(anglr)
library(ozmaps)
lga <- ozmap_data("abs_lga") 
x <- dplyr::filter(lga, sf::st_coordinates(sf::st_centroid(lga))[,2] < -40)
#x <- sf::read_sf(system.file("gpkg/nc.gpkg", package = "sf", mustWork = TRUE))
## get the boundaries as segments
s <- SC0(sf::st_union(x))
segs <- as.matrix(do.call(rbind, lapply(s$object$topology_, 
                      function(x) x[c(".vx0", ".vx1")])))
tri <- TRI0(x)
mesh <- as.mesh3d(tri)
nv <- cbind(rbind(t(as.matrix(s$vertex)), 
                  z = 0, h = 1), rbind(t(as.matrix(s$vertex)), 
              z = -1, h = 1))
n <- dim(mesh$vb)[2]
quads <- cbind(segs+n, .vx2 = segs[,2] + n + nrow(segs), .vx3 = segs[,1] + n + nrow(segs))

mesh$ib <- t(quads)
mesh$vb <- cbind(mesh$vb, nv)
rl_primitives <- c(unlist(lapply(tri$object$topology_, nrow)), nrow(quads))
plot3d(mesh, color = rep(viridis::viridis(length(rl_primitives)), rl_primitives), 
       alpha = 0.6); rgl::aspect3d(1, 1, .1)

Uluru

https://twitter.com/mdsumner/status/1251800819495694337/photo/1

cc <- ceramic::cc_location(cbind(131.037, -25.3456562),
                           buffer = c(1800, 1200), 
                           zoom = 17)
el <- ceramic::cc_elevation(cc)
library(anglr)
mesh <- as.mesh3d(el, image_texture = cc)
clear3d(type= "lights")
clear3d(type= "bboxdeco")

light3d(45, 75)
plot3d(mesh, specular = "black")
rgl::aspect3d("iso")
rgl::bg3d("#111111")

Hawaii

https://twitter.com/mdsumner/status/1251807688398233600

cc <- ceramic::cc_location(cbind(-155.45, 19.61),
                           buffer = c(85000, 85000), zoom = 11)
#raster::plotRGB(cc)
el <- ceramic::cc_elevation(cc, zoom = 10)
library(anglr)
mesh <- as.mesh3d(el, image_texture = cc)
library(rgl)
plot3d(mesh, specular = "black", axes = FALSE, b)
clear3d(type= "lights")
clear3d(type= "bboxdeco")
light3d(45, 75)
rgl::aspect3d(1, 1, .1)
rgl::bg3d("#111111")

Precipitous Bluff

https://twitter.com/mdsumner/status/1251810265848442880/photo/1

cc <- ceramic::cc_location(cbind(146.57, -43.48),
                           buffer = c(8000, 11000), zoom = 14)
raster::plotRGB(cc)
el <- ceramic::cc_elevation(cc, zoom = 13)
library(anglr)
mesh <- as.mesh3d(el, image_texture = cc)
library(rgl)
plot3d(mesh, specular = "black", axes = FALSE)
clear3d(type= "lights")
clear3d(type= "bboxdeco")
light3d(45, 75)
rgl::aspect3d("iso")
rgl::bg3d("#111111")

Tasmania

This examples a bit more involved, with a wire3d mesh DEM and then polygons with the image texture hovering over the mesh - the polygons have a white outline.

library(sf)
library(ceramic)
library(anglr) ## remotes::install_github("hypertidy/anglr")
mp <- ozmaps::abs_lga
mp <- mp[st_coordinates(st_centroid(mp))[,2] < -39, ]
plot(mp)
cc <- ceramic::cc_location(raster::extent(mp), zoom  = 8)
#raster::plotRGB(cc)
el <- ceramic::cc_elevation(cc, zoom = 5)
library(anglr)
mesh <- as.mesh3d(el)
library(rgl)
mesh$material$color <- "grey"
wire3d(mesh, specular = "black", axes = FALSE)
im <- silicate::TRI0(mp)
clear3d(type= "bboxdeco")
# clear3d(type= "lights")
# light3d(0, 75)
rgl::aspect3d(1, 1, .08)

rgl::bg3d("#111111")
sc <- reproj::reproj(DEL0(mp), 
                     raster::projection(cc))
ln <- reproj::reproj(SC0(mp), 
                     raster::projection(cc))
ln$vertex$z_ <- 1700
sc$vertex$z_ <- 1700
sc$object$color_ <- sample(viridis::viridis(nrow(sc$object)))
scmesh <- as.mesh3d(sc, image_texture = cc)
scmesh$material$color <- "white"
plot3d(scmesh, add = TRUE, alpha = 0.95)
plot3d(ln, add = TRUE, color= "white", size = 2)

Heard Island

https://twitter.com/mdsumner/status/1251847551818412032/photo/1

pt <- cbind(73 + 31/60, -(53 + 6/60))

#pt <- cbind(78 + 15/60, -(68 + 33/60))
cc <- ceramic::cc_location(pt,
                           buffer = c(45000, 45000))

el <- ceramic::cc_elevation(cc)
library(anglr)
mesh <- as.mesh3d(el, image_texture = cc)
raster::plotRGB(cc)
plot3d(mesh, specular = "black")
clear3d(type= "lights")
clear3d(type= "bboxdeco")
light3d(45, 75)

rgl::aspect3d("iso")
rgl::bg3d("#111111")

Geoscience Australian Topography on globe

This now with adaptive-density triangulation thanks to terrainmeshr/hmm

https://twitter.com/mdsumner/status/1255702228171546624/photo/1

Original quadmesh version: https://twitter.com/mdsumner/status/1195281139096621056/photo/1

library(ceramic)
library(anglr) ## remotes::install_github("hypertidy/anglr")
library(wmts) ## remotes::install_github("mdsumner/wmts")
library(rgl)

dem <- cc_elevation(raster::extent(110, 155, -45, -10), 
                    max_tiles = 36)

url <- "http://gaservices.ga.gov.au/site_7/rest/services/Topographic_Base_Map_WM/MapServer/WMTS/1.0.0/WMTSCapabilities.xml"
dsn <- glue::glue("WMTS:{url},layer=Topographic_Base_Map_WM,tilematrixset=default028mm")
img <- wmts(dsn, loc = dem, max_tiles = 64)
x <- DEL0(dem * 18, max_triangles = 60000)

qm <- as.mesh3d(x, image_texture = img)
qm$vb[1:3, ] <- t(reproj::reproj(t(qm$vb[1:3, ]), "+proj=geocent", source = crsmeta::crs_proj(img)))
plot3d(qm, specular = "darkgrey"); bg3d("black")

Mars

https://twitter.com/mdsumner/status/1265263093052932096/photo/1

dangerurl <- "/vsicurl/https://planetarymaps.usgs.gov/mosaic/Mars/HRSC_MOLA_Blend/Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif"

library(lazyraster)
library(anglr)
lr <- lazyraster(dangerurl)
ar <- as_raster(lr, c(1024, 512) * 2)
d <- as.mesh3d(ar)
 
d$vb[1:3, ] <- t(anglr:::.ll_to_globe(t(d$vb[1:3, ]), exag = 20))
plot3d(d, col = "white")

Stacked plot

Adapted from https://urbandemographics.blogspot.com/2016/04/creating-tilted-and-stacked-maps-in-r.html

See https://twitter.com/mdsumner/status/1273947106525364232?s=20

# load libraries
library(rgeos)
library(UScensus2000tract)
library(dplyr)
library(RColorBrewer)

# load data
data("oregon.tract")


library(anglr)

## triangles from the polygons
tri <- TRI0(oregon.tract)

## base layer
bounds <- PATH0(oregon.tract)
mesh <- as.mesh3d(tri)

## population
cols <- RColorBrewer::brewer.pal(11, "RdBu")
tri2 <- TRI0(copy_down(tri, 1.5))## elevate arbitrarily
tri3 <- TRI0(copy_down(tri, 3))

## not exact, and bit awkward but it's something
tricolors3 <- rep(colourvalues::colour_values(tri$object$pop2000 - tri$object$white, palette = t(col2rgb(rev(cols)))), 
                 unlist(lapply(tri2$object$topology_, nrow)))

tricolors2 <- rep(colourvalues::colour_values(tri$object$white, palette = t(col2rgb(rev(cols)))), 
                  unlist(lapply(tri2$object$topology_, nrow)))

plot3d(mesh, col = "white")
plot3d(bounds, add = TRUE, color = "black")

plot3d(tri2, add =TRUE, color = tricolors2)
plot3d(tri3, add =TRUE, color = tricolors3)

Stacked walls plot

https://twitter.com/mdsumner/status/1273985235814764544/photo/1

# load libraries
library(rgeos)
library(UScensus2000tract)
library(dplyr)
library(RColorBrewer)

# load data
data("oregon.tract")


library(anglr)

## triangles from the polygons
tri <- TRI0(oregon.tract)

## base layer
bounds <- PATH0(oregon.tract)
mesh <- as.mesh3d(tri)

## population
cols <- RColorBrewer::brewer.pal(11, "RdBu")
scale2 <- ggplot2::scale_fill_distiller(palette = "RdBu", limits = c(range(oregon.tract$pop2000)))

## not exact, and bit awkward but it's something
idx <- unlist(lapply(tri2$object$topology_, nrow))
tricolors2 <- rep(scale2$map(oregon.tract$pop2000), idx)
tricolors3 <- rep(scale2$map(oregon.tract$pop2000 - oregon.tract$white), idx)
flip <- 
function(x) {
  xx <- x - min(x)
  max(xx) - xx
}


scale2 <- ggplot2::scale_fill_distiller(palette = "RdBu", limits = c(range(oregon.tract$pop2000)))

tri2 <- tri3 <- tri

tri2$vertex <- 
  tri2$vertex %>% dplyr::mutate(x0 = x_, y0 = y_) %>% 
  transmute(x_ =  min(x0), y_= scales::rescale(x0, to = range(y_)), z_ = scales::rescale(y0))

tri3$vertex <- 
  tri3$vertex %>% dplyr::mutate(x0 = x_, y0 = y_) %>% 
  transmute(x_ =  x_, y_= max(y_), z_ = scales::rescale(y0)) 

plot3d(mesh, col = "white")
plot3d(bounds, add = TRUE, color = "black")

plot3d(tri2, add =TRUE, color = tricolors2)

plot3d(tri3, add =TRUE, color = tricolors3)

rgl::aspect3d(1, 1, .5)
Clone this wiki locally