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()

lltopo <- quadmesh(etopo, etopo)
lltopo$vb[1:3, ] <- t(llh2xyz(t(lltopo$vb[1:3, ]), exag = 50))
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, ])))
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))

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

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

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

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))

## quadmesh::mesh_plot() is to be deprecated
##  a future release will use the 'anglr' package, soon to be on CRAN

xy <- proj4::project(xymap, prj)
lines(xy)