Topography with quadmesh

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)
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))
qmtopo <- quadmesh(etopo, etopo)

If you have the rgl package installed you can run the following code to produce an interactive plot.

library(rgl)
open3d()
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.

If you have the rgl package installed you can run the following code to produce an interactive plot.

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()
shade3d(qmtopo, col = grey(scl(qmtopo$vb[3,qmtopo$ib])))
aspect3d(1, 1, .1)

#rglwidget()
lltopo <- quadmesh(etopo, etopo)
lltopo$vb[1:3, ] <- t(llh2xyz(t(lltopo$vb[1:3, ]), exag = 50))

If you have the rgl package installed you can run the following code to produce an interactive plot.

shade3d(lltopo, col = grey(scl(qmtopo$vb[3,lltopo$ib])))
#rglwidget()

Try a local area, this is not rendered in order to keep the vignette size down.

data(etopo)
lt <- raster::crop(etopo, raster::extent(100, 168, -58, -40))
ltt <- ltt0 <- quadmesh(lt, lt)
ltt$vb[1:3, ] <- t(llh2xyz(t(ltt$vb[1:3, ])))

If you have the rgl package installed you can run the following code to produce an interactive plot.


open3d()
shade3d(ltt, col = grey(scl(ltt0$vb[3,ltt$ib])))
#rglwidget()

The package also includes a mesh_plot() function for straightforward plotting of a raster to an arbitrary map projection.

library(raster)
plot(etopo, col = palr::bathy_deep_pal(20))
prj <- "+proj=stere +lat_ts=-71 +lat_0=-90 +lon_0=147 +datum=WGS84"
mesh_plot(etopo, crs = prj, asp = 1, col = palr::bathy_deep_pal(20))
xy <- reproj::reproj(xymap, prj, source = 4326)
lines(xy)