TRI for raster with hmm

[mdsumner]

Notes for terrainmeshr use:

• use sentinels for missing values, not handled atm
• you only get row and col coordinates

Example

r <- anglr::gebco
heightmap <- raster::as.matrix(r)
## if missing data you have to sentinelize them
#heightmap[is.na(heightmap)] <- -10000
m <- terrainmeshr::triangulate_matrix(heightmap, maxTriangles = 256)
image(1:nrow(heightmap), 1:ncol(heightmap), heightmap)
plot(r)
## respatialize
cell <- cellFromRowCol(r, m[,1], m[,3])
d <- tibble::tibble(x = xFromCell(r, cell), y = yFromCell(r, cell), 
                    z = m[,2])
idx <- cbind(matrix(seq_len(nrow(d)), ncol = 3, byrow = TRUE), NA)
polygon(d[head(as.vector(t(idx)), -1), c("x", "y")])

[mdsumner]

A mesh3d example

library(raster)
r <- anglr::gebco
heightmap <- raster::as.matrix(r)
## if missing data you have to sentinelize them
#heightmap[is.na(heightmap)] <- -10000
m <- terrainmeshr::triangulate_matrix(heightmap, maxTriangles = 40000, y_up = TRUE)
#image(1:nrow(heightmap), 1:ncol(heightmap), heightmap)
#plot(r)
## respatialize
cell <- cellFromRowCol(r, m[,1], m[,3])
d <- tibble::tibble(x = 180 - xFromCell(r, cell), y = yFromCell(r, cell), 
                    z = m[,2])
idx <- cbind(matrix(seq_len(nrow(d)), ncol = 3, byrow = TRUE), NA)
#polygon(d[head(as.vector(t(idx)), -1), c("x", "y")])

library(rgl)
xyz <- as.matrix(d)

mesh <- tmesh3d(rbind(t(as.matrix(d)), 1), t(idx[,1:3]), 
                material = list(color = colourvalues::colour_values(d$z)))
#triangles3d(d$x, d$y, d$z)
mesh$vb[1:3, ] <- t(quadmesh::llh2xyz(t(mesh$vb[1:3, ]), rad = 637813))
plot3d(mesh, specular = "darkgrey")

[mdsumner]

Oops, note the 180- hack in there to get the orientation right, have to fix that up

[mdsumner]

This works, needs some thought:

• normalize vertices (see ?rgl::as.mesh3d)
• put in anglr, and use DEL instead? (or hack around raster and avoid import)

#' @name TRI0
#' @export
TRI0.BasicRaster <- function(x, ..., max_triangles = NULL) {
  heightmap <- t(raster::as.matrix(x))[,nrow(x):1]
  if (is.null(max_triangles)) max_triangles <- prod(dim(heightmap))/20
  ## if missing data you have to sentinelize them
  sentinel <- min(heightmap, na.rm  = TRUE) -100
  heightmap[is.na(heightmap)] <- sentinel
  hmm <- terrainmeshr::triangulate_matrix(heightmap, maxTriangles = max_triangles)
  ## respatialize
  cell <- raster::cellFromRowCol(r, hmm[,3], hmm[,1])
#  n_triangles <- dim(hmm)[1]/3

  topology <- tibble::tibble(.vx0 = seq(1, dim(hmm)[1], by = 3L),
                             .vx1 = .vx0 + 1L,
                             .vx2 = .vx1 + 1L)
  meta <- tibble(crs = raster::projection(x), ctime = Sys.time())
  structure(list(object = tibble(a = 1, topology_ = list(topology)),
                 vertex = tibble::tibble(x_ = hmm[,1], y_ = hmm[,3], z_ = hmm[,2]),
                 meta = meta),
            class = c("TRI0", "sc"))
}

[mdsumner]

Now with missing value handling

#' @name TRI0
#' @export
TRI0.BasicRaster <- function(x, ..., max_triangles = NULL) {
  heightmap <- t(raster::as.matrix(x))[,nrow(x):1]
  if (is.null(max_triangles)) max_triangles <- prod(dim(heightmap))/20
  ## if missing data you have to sentinelize them

  dosentinel <- FALSE
  if (anyNA(heightmap)) {
    dosentinel <- TRUE
    sentinel <- as.integer(min(heightmap, na.rm  = TRUE) -100)
    heightmap[is.na(heightmap)] <- sentinel
  }
  hmm <- terrainmeshr::triangulate_matrix(heightmap, maxTriangles = max_triangles)
  if (dosentinel) {
    #browser()
    hmm <- hmm[hmm[,2] > sentinel, ]
    bad <- table(hmm[,4]) < 3
    if (any(bad)) {
      hmm <- hmm[!hmm[,4] %in% as.integer(names(bad[which(bad)])), ]
    }
  }
  ## respatialize
  cell <- raster::cellFromRowCol(r, hmm[,3], hmm[,1])

  topology <- tibble::tibble(.vx0 = seq(1, dim(hmm)[1], by = 3L),
                             .vx1 = .vx0 + 1L,
                             .vx2 = .vx1 + 1L)
  meta <- tibble(crs = raster::projection(x), ctime = Sys.time())
  structure(list(object = tibble(a = 1, topology_ = list(topology)),
                 vertex = tibble::tibble(x_ = hmm[,1], y_ = hmm[,3], z_ = hmm[,2]),
                 meta = meta),
            class = c("TRI0", "sc"))
}

[mdsumner]

Done in anglr::DEL0, it's Delaunay so it makes sense.

https://github.com/hypertidy/anglr/tree/hmm