---
title: "Topography with quadmesh"
author: "Michael D. Sumner"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Topography with quadmesh}
  %\VignetteEncoding{UTF-8}
  %\VignetteEngine{knitr::rmarkdown}
editor_options: 
  chunk_output_type: console
---

# 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. 

```{r}
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))
qmtopo <- quadmesh(etopo, etopo)
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. 

```{r}

library(proj4)
qmtopo$vb[1:2, ] <- t(proj4::project(t(qmtopo$vb[1:2, ]), "+proj=laea +datum=WGS84 +lat_0=-90"))
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))
shade3d(lltopo, col = grey(scl(qmtopo$vb[3,lltopo$ib])))

```


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

```{r eval=FALSE, include=TRUE}
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. 

```{r}
library(raster)
plot(etopo, col = palr::bathyDeepPal(20))
prj <- "+proj=stere +lat_ts=-71 +lat_0=-90 +lon_0=147 +datum=WGS84"
mesh_plot(etopo, crs = prj, asp = 1, colfun = palr::bathyDeepPal)
xy <- proj4::project(xymap, prj)
lines(xy)


```
