Examples
https://twitter.com/mdsumner/status/1312960193982996480/photo/1
image from https://gpm.nasa.gov/data/imerg
br <- raster::brick("imergert_1080p.now.png")
library(raster)
library(anglr)
extent(br) <- extent(-180, 180, -90, 90)
mesh <- as.mesh3d(br[[1]] * 0, image_texture = br, triangles = FALSE)
mesh$vb[1:3, ] <- t(quadmesh::llh2xyz(t(rbind(mesh$vb[1:2, ], 0))))
plot3d(mesh, specular = "black")
rgl::bg3d("black")
rgl::clear3d("bboxdeco")https://twitter.com/mdsumner/status/1311627020825055233/photo/1
library(string2path)
library(ggplot2)
# This TTF file is downloaded from https://ipafont.ipa.go.jp/.
# For installed fonts, you can use systemfonts::system_fonts()
# to lookup the path.
d <- string2path("広島はこの時期に美しい", "C:/temp/ipam.ttf")
d <- tibble::rowid_to_column(d)
p <- silicate::PATH0_from_df(d %>% mutate(path_ = id))
im <- ceramic::cc_location(cbind(132.4553, 34.3853), buffer = 30000)
el <- ceramic::cc_elevation(im)
plotRGB(im)
extent(im) <- extent(1, 100, 20, 100)
extent(el) <- extent(im)
plot3d(as.mesh3d(el, triangles = FALSE, image_texture = im), specular = "black")
rgl::clear3d("bboxdeco")
plot3d(DEL0(p), add = TRUE, color = "mauve")
d <- string2path("R", "C:/temp/ipam.ttf")
https://twitter.com/mdsumner/status/1311615172746768384?s=20
isoband_raster <- function(x, lo, hi, auto = FALSE) {
if (auto) {
breaks <- pretty(values(x))
lo <- head(breaks, -1)
hi <- tail(breaks, -1)
}
## OMG: note the [[1]] to avoid the case of as.matrix(brick) which is not helpful ...
b <- isoband::isobands(xFromCol(x), yFromRow(x), as.matrix(x[[1]]), levels_low = lo, levels_hi = hi)
sf::st_sf(lo = lo, hi = hi, geometry = sf::st_sfc(isoband::iso_to_sfg(b), crs = projection(x)))
}
library(SOmap)
d <- aggregate(SOmap::Bathy, fact = 4)
bks <- c(cellStats(d, "min"), c(-5000, -4000, -3000, -2500, -2000, -1500, -1000, -500, 0))
```polys <- isoband_raster(d, head(bks, -1), tail(bks, -1))
ramp2 <- grDevices::colorRampPalette(c("#54A3D1","#60B3EB","#78C8F0","#98D1F5","#B5DCFF","#BDE1F0","#CDEBFA","#D6EFFF","#EBFAFF","grey92","grey94","grey96", "white"))
bluepal<-ramp2(nrow(polys))
plot(polys[1], border = NA, col = bluepal)
bluepal<- head(ramp2(nrow(polys) + 3), -3)
plot(polys[1], border = NA, col = bluepal)
library(anglr)
mesh <- DEL0(polys, max_area = 1e10)
mesh$object$color_ <- bluepal
x <- as.mesh3d(mesh)
x$vb[3, ] <- extract(SOmap::Bathy, t(x$vb[1:2, ]))
plot3d(x)
rgl::aspect3d(1, 1, .01)
# Sakurajima
apologies to @researchremora
https://twitter.com/mdsumner/status/1309068366477496320?s=20
```R
library(ceramic)
library(anglr)
library(raster)
library(rgl)
#https://gbank.gsj.jp/volcano/Act_Vol/sakurajima/map/eng/volcmap01e.html
img <- raster::brick("map_l.jpg")
## this georeg is pretty sloppy, coordinates on the image margin are
## tl, br
#xy <- rbind(c(130 + 35/60 + 20/3600, 31 + 38/60 + 50/3600),
#c(130 + 44/60 + 40/3600, 31 + 32/60 + 20 /3600))
extent(img) <- extent(130.587, 130.7452, 31.53746, 31.64785)
projection(img) <- "+proj=longlat"0
el <- cc_location(cbind(130.6667171,31.5833836), 7500, type = "elevation-tiles-prod", zoom = 12)
mesh <- as.mesh3d(el, image_texture = img, triangles = F)
plot3d(mesh)
aspect3d("iso")
clear3d("bboxdeco")https://twitter.com/mdsumner/status/1294818713300942848?s=20
to_polygon <- function(x, y) {
list(cbind(x, y)[c(seq_along(x), 1), ])
}
cut <- sf::st_sfc(sf::st_multipolygon(list(
to_polygon(rep(c(80.01, 79.99), each = 91), c(-90:0, 0:-90)), # third cut bottom
to_polygon(rep(c(-19.99, -20.01), each = 91), c(-90:0, 0:-90)), # second cut bottom)
to_polygon(rep(c(-99.99, -100.01), each = 91), c(-90:0, 0:-90)), # first cut bottom)
to_polygon(rep(c(-40.01, -39.99), each = 91), c(90:0, 0:90) )
)
), crs = 4326)
data("wrld_simpl", package = "maptools")
## we have to buffer, else topology problems (other sources are not perfect either)
world_sf <- sf::st_buffer(sf::st_as_sf(wrld_simpl), 0)
library(sf)
x <- st_difference(world_sf$geometry,cut)
#plot(st_transform(x, "+proj=igh"))
gr <- graticule::graticule(seq(-180, 180, by = 10), seq(-90, 90, by = 10), nverts = 120, tiles = TRUE)
grc <- st_difference(st_as_sf(gr), cut)
#plot(st_transform(grc, "+proj=igh"))
## avoid sf preventing overplot
plot(st_geometry(st_transform(grc, "+proj=igh")), col = rgb(0.95, 0.95, 0.95, 0.5), border = "transparent")
plot(st_cast(st_transform(grc, "+proj=igh"), "LINESTRING"), add = TRUE, col = rgb(0.65, 0.65, 0.65))
plot(st_transform(x, "+proj=igh"), add = TRUE, col = "transparent", border = "black")
## now mesh and plot in rgl
library(anglr)
im <- ceramic::cc_location(cbind(0, 0), buffer = 6378137 * pi)
mesh <- anglr::as.mesh3d(st_transform(x, "+proj=igh"), image_texture = im)
#mesh_plot(mesh)
plot3d(mesh)
plot3d(st_transform(grc, "+proj=igh"), add = TRUE, color = "darkgrey")
rgl::clear3d("bboxdeco")
rgl::play3d(rgl::spin3d())- 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)
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)
https://twitter.com/mdsumner/status/1283161616779866112/photo/1
library(anglr)
library(raster)
library(ceramic)
ex <- extent(-82, -37, -66, 12)
#topo <- crop(gebco, ex)
topo <- cc_location(ex, type = "elevation-tiles-prod", max_tiles = 64)
img <- cc_location(topo, base_url = "https://stamen-tiles-a.a.ssl.fastly.net/terrain/{zoom}/{x}/{y}.png", zoom = 5)
mesh <- as.mesh3d(topo, image_texture = img)
## convert to geocentric (might need to do this manually given reproj-CRAN status yep # mesh <- globe(mesh)
mesh$vb[1:3, ] <- t(quadmesh::llh2xyz(reproj::reproj(t(mesh$vb[1:3, ]), source = projection(topo), target = 4326),
exag = 10))
plot3d(mesh)
## zap on the country boundaries
ne <- subset(rnaturalearth::ne_countries(), continent == "South America")
pts <- silicate::sc_vertex(ne)
rgl::points3d(quadmesh::llh2xyz(cbind(pts$x_, pts$y_, 200)), color = "hotpink", size = 5)
rgl::bg3d("black")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")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")
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")
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)
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")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")
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")Bit more of a serious crack at this
https://twitter.com/mdsumner/status/1281217016653733894/photo/1
dangerurl <- "https://planetarymaps.usgs.gov/mosaic/Mars/HRSC_MOLA_Blend/Mars_HRSC_MOLA_BlendDEM_Global_200mp_v2.tif"
library(lazyraster)
## note prefix with GDAL's URL indirection
mars <- lazyraster(file.path("/vsicurl", dangerurl))
#plateau <- crop(mars, raster::extent(-148, -15, -47, 38))
#plot(plateau, col = grey(seq(0, 1, length = 256)))
library(anglr)
r <- as_raster(mars, c(3600, 3600))
mesh <- as.mesh3d(r, triangles = F)
mesh$vb[1:3, ] <- t(quadmesh::llh2xyz(t(mesh$vb[1:3, ]), rad = 3389500, exag = 5))
plot3d(mesh, col = "white", specular = "black")
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)
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)