whew, nearly

[mdsumner]

library(sf)
sfx <- silicate::inlandwaters[6, ]
library(geos)
library(wk)
tri <- wk_flatten(geos_constrained_delaunay_triangles(geos_densify(wk_flatten(geos_constrained_delaunay_triangles(geos_densify(wk_flatten(geos_constrained_delaunay_triangles(geos_densify(
wk_flatten(geos_constrained_delaunay_triangles(sfx)), 50000))), 10000))), 200)))


library(dsn)
library(vapour)
bb <- wk::wk_bbox(tri)
elev <- gdal_raster_data(dsn::gebco21(), target_ext =as.numeric(bb)[c(1, 3, 2, 4)], target_res = c(40, 40), target_crs = wk::wk_crs(bb)$wkt)

xy <- wk::wk_coords(tri)
tri_idx <- seq_len(nrow(xy))
tri_idx <- tri_idx[!(tri_idx %% 4 == 0)]
## so now we can get the triangles, their vertices, the elevation of those, and plot
cell <- vaster::cell_from_xy(attr(elev, "dimension"), attr(elev, "extent"), xy[, c("x", "y"), drop  =FALSE])
xy <- cbind(xy, z =  elev[[1]][cell])
##now plot!
library(rgl)
mesh <- tmesh3d(t(xy[, c("x", "y", "z"), drop = FALSE]), matrix(tri_idx, nrow = 3))
clear3d()
shade3d(mesh, col = "white")
aspect3d(1, 1, .1)
bg3d("darkgrey")