library(anglr)
library(raster)

## Loading required package: sp

x <- raster::crop(gebco1, raster::extent(100, 160, -60, -20))
x <- aggregate(x, fact = 3)


## this decouples the order of the coordinates from the values, and can be done without any triangulation library
## just using indexing (as quadmesh does) but I just don't have the brain space for that today
tri <- RTriangle::triangulate(RTriangle::pslg(P = cbind(raster::colFromCell(x, seq_len(raster::ncell(x))),
                                                        nrow(x) - raster::rowFromCell(x, seq_len(raster::ncell(x))) + 1)))


prj <- "+proj=laea +lon_0=140 +lat_0=-40 +datum=WGS84"
coords <- rgdal::project(coordinates(x), prj)
o <- structure(list(vb = t(cbind(coords, raster::values(x), 1)),
                    it = t(tri$T),
                    primitivetype = "triangle",
                    material = list(specular = "black"),
                    normals = NULL,
                    texcoords = NULL), class = c("mesh3d", "shape3d"))

rgl::open3d()

## null 
##    1

rgl::clear3d()

cols <- viridis::viridis(56)[scales::rescale(o$vb[3,o$it], c(1, 56))]

rgl::shade3d(o, col = cols)
rgl::aspect3d(1, 1, .0001)  ## don't set z-exag to 0 ...  
qm <- quadmesh::quadmesh(x)
qm$vb[1:2, ] <- t(rgdal::project(t(qm$vb[1:2, ]), prj))
rgl::wire3d(qm, col = "red")

rgl::rglwidget()

## equivalent in 2D (quadmesh is not doing the right interpolation currently ...)
##quadmesh::mesh_plot(x, crs = prj)
plot(coords, type = "n", asp = 1)

XY <- coords[t(tri$T), ]
ID <- rep(1:nrow(tri$T), each = 3)
                                                     ## watch out here, we need a triangle vertex to get the colours in the right order
COL <- viridis::viridis(56)[scales::rescale(values(x)[tri$T[,1]], c(1, 56))]
isLL <- FALSE
#graphics::plot.new()
#graphics::plot.window(xlim = range(XY[,1], finite = TRUE), ylim = range(XY[,2], finite = TRUE), asp = if (isLL) 1/cos(mean(x$y, na.rm = TRUE) * pi/180) else 1  )
    
vps <- gridBase::baseViewports()
grid::pushViewport(vps$inner, vps$figure, vps$plot)
grid::grid.polygon(XY[,1], XY[,2], ID,  gp = grid::gpar(col = NA, fill = COL),
                   default.units = "native")
grid::popViewport(3)
## these lines are over-drawn
lines(t(qm$vb)[rbind(qm$ib, qm$ib[1, ], NA), ], col = rgb(0, 0, 0, 0.4))

## that actually looks bad, but it's because we only get discrete values on triangles in grid.polygon, but
## we get within-triangle colour interpolotion in rgl
## in real data it looks much better because the triangulation is denser: 

x <- raster::crop(gebco1, raster::extent(30, 180, -80, 20))
x

## class       : RasterLayer 
## dimensions  : 240, 360, 86400  (nrow, ncol, ncell)
## resolution  : 0.4166667, 0.4166667  (x, y)
## extent      : 30, 180, -80, 20  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : in memory
## names       : gebco1 
## values      : -8173, 3983  (min, max)

tri <- RTriangle::triangulate(RTriangle::pslg(P = cbind(raster::colFromCell(x, seq_len(raster::ncell(x))),
                                                        nrow(x) - raster::rowFromCell(x, seq_len(raster::ncell(x))) + 1)))


prj <- "+proj=laea +lon_0=75 +lat_0=-40 +datum=WGS84"
coords <- rgdal::project(coordinates(x), prj)

plot(coords, type = "n", asp = 1)

XY <- coords[t(tri$T), ]
ID <- rep(1:nrow(tri$T), each = 3)
## watch out here, we need a triangle vertex to get the colours in the right order
COL <- viridis::viridis(56)[scales::rescale(values(x)[tri$T[,1]], c(1, 56))]
isLL <- FALSE
#graphics::plot.new()
#graphics::plot.window(xlim = range(XY[,1], finite = TRUE), ylim = range(XY[,2], finite = TRUE), asp = if (isLL) 1/cos(mean(x$y, na.rm = TRUE) * pi/180) else 1  )

vps <- gridBase::baseViewports()
grid::pushViewport(vps$inner, vps$figure, vps$plot)
grid::grid.polygon(XY[,1], XY[,2], ID,  gp = grid::gpar(col = NA, fill = COL),
                   default.units = "native")
grid::popViewport(3)