Topographic example

When we plot global data sets in 3D it can be very weird, since the data is stored with longitude-latitude angular coordinates, but the data are in metres above mean sea level. To make this work we need an arbitrary scaling to reduce the distortion due to these different units.

library(quadmesh)
library(rgl)
scl <- function(x) (x - min(x))/diff(range(x))
data(etopo)
## very low resolution, simply to keep vignette size down
etopo <- raster::setExtent(raster::aggregate(etopo, fact = 10), raster::extent(etopo))

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum WGS_1984 in CRS definition,
##  but +towgs84= values preserved

qmtopo <- quadmesh(etopo, etopo)
open3d()

## null 
##   10

shade3d(qmtopo, col = grey(scl(qmtopo$vb[3,qmtopo$ib])), asp = c(1, 1, 0.0001))

aspect3d(1, 1, 0.1)


rglwidget()

It’s much more natural to work in a map projection, or on the surface of the sphere, but legacy spatial tools make this more difficult than it should be.

We can do what ever we like with quadmesh coordinates, without breaking the topology of the mesh at all.

qmtopo$vb[1:2, ] <- t(reproj::reproj(t(qmtopo$vb[1:2, ]), "+proj=laea +datum=WGS84 +lat_0=-90", source = "+proj=longlat +datum=WGS84")[,1:2, drop = FALSE])
open3d()

## null 
##   12

shade3d(qmtopo, col = grey(scl(qmtopo$vb[3,qmtopo$ib])))
aspect3d(1, 1, .1)

rglwidget()