raster-as-tri-mesh.r

quadToTriangle <- function(x) {
  v <- x$v
  v$vertex_ <- seq(nrow(v))
  meta <- x$meta
  tab <- x$qXv
  n4 <- nrow(tab) / 4L
  tXv <- tibble(vertex_ = x$qXv$vertex_[rep(c(1, 2, 3, 1, 3, 4), n4) + rep(seq(1, length = n4, by = 4)-1, each = 6)])
  tXv$triangle_ <- rep(seq(nrow(tXv)/3), each = 3)
  x <- list(o = tibble(object_ = "1"), t = tibble(triangle_ = seq(nrow(tXv)/3), object_ = "1"),
            tXv = tXv, v = v, meta = meta)
  class(x) <- "trimesh"
  x
}

library(rangl)
library(raster)

## Loading required package: sp

library(marmap) ## has topo data

## 
## Attaching package: 'marmap'

## The following object is masked from 'package:raster':
## 
##     as.raster

## The following object is masked from 'package:grDevices':
## 
##     as.raster

data(hawaii)
library(rgeos)

## rgeos version: 0.3-21, (SVN revision 540)
##  GEOS runtime version: 3.5.0-CAPI-1.9.0 r4084 
##  Linking to sp version: 1.2-3 
##  Polygon checking: TRUE

haw <- setExtent(flip(t(raster(unclass(hawaii))), "y"), extent(unlist(lapply(attr(hawaii, "dimnames"), function(x) as.numeric(c(head(x, 1), tail(x, 1)))))))
projection(haw) <- "+proj=longlat +ellps=WGS84"
qmesh <- globe(rangl(haw * 280))
library(tibble)
mesh <- plot(quadToTriangle(qmesh))

## Joining, by = "object_"

## Joining, by = "triangle_"

library(plotly)

## Loading required package: ggplot2

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:raster':
## 
##     select

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

# plot point cloud
x <- mesh$vb[1,]
y <- mesh$vb[2, ]
z <- mesh$vb[3,]
m <- matrix(c(x,y,z), ncol=3, dimnames=list(NULL,c("x","y","z")))

# now figure out the colormap
zmean <- apply(t(mesh$it),MARGIN=1,function(row){mean(m[row,3])})

library(scales)
facecolor = colour_ramp(
  brewer_pal(palette="RdBu")(9)
)(rescale(x=zmean))

plot_ly(
  x = x, y = y, z = z,
  i = mesh$it[1,]-1, j = mesh$it[2,]-1, k = mesh$it[3,]-1,
  facecolor = facecolor,
  type = "mesh3d"
)