volcano example

[mdsumner]

  ## according to mapview
maxll <- cbind(174.76448,-36.87812)

prj <- "+proj=nzmg +datum=WGS84"
maxnz <- rgdal::project(maxll, prj)
library(raster)
r <- raster::raster(volcano[nrow(volcano):1,ncol(volcano):1])
plot(r)
wm <- which.max(r)
xmn <- maxnz[1,1] - 10 * colFromCell(r, wm)
xmx <- maxnz[1,1] + 10 * (ncol(r) - colFromCell(r, wm))
ymn <- maxnz[1,2] - 10 * (nrow(r) - rowFromCell(r, wm))
ymx <- maxnz[1,2] + 10 * rowFromCell(r, wm)
rr <- setExtent(r, extent(xmn, xmx, ymn, ymx))

projection(rr) <- prj
plot(rr)
abline(v = maxnz[1,1], h = maxnz[1,2])
#mapview::mapview(rasterToContour(rr))
im <- ceramic::cc_location(rr, zoom = 18)  ## hardcode zoom
plotRGB(im)
plot(spTransform(rasterToContour(rr), projection(im)), add = TRUE, 
     lwd = 3, col = "lightgrey")




mesh <- anglr::as.mesh3d(rr, image_texture = im)
anglr::plot3d(mesh);rgl::aspect3d("iso")
cl <- rasterToContour(rr)
cl$level <- as.numeric(levels(cl$level)[cl$level]) + 1 #(add a fudge to elevate)
anglr::plot3d(anglr::copy_down(silicate::SC(cl), "level"), add = TRUE, lwd = 3)

With decent zoom from ceramic it looks pretty neat

not sure what to do about the texture problem, save the mesh3d and hard code it to look up the system.file?