example code - raster to quadmesh to trimesh

[mdsumner]

Here's some functions that might be of interest, to convert raster directly to a trimesh, probably you'd crop first but I was happy to see that it worked on the nt_raster example. There is a little more detail in the raster data, and it's hard to see that without comparing like this.

quadmesh_to_trimesh <- function(x, ...) {
  n4 <- ncol(x$ib)
  idx <- rep(c(1, 2, 3, 1, 3, 4), n4) + rep(seq(1, length = n4, by = 4)-1, each = 6)
  x$it <- array(x$ib[idx], c(3, length(idx)/3))
  x$primitivetype <- "triangle"
  x$ib <- NULL  ## remove quad interpretation
  x
}

raster_to_trimesh <- function(x, ...) {
  quadmesh_to_trimesh(quadmesh::quadmesh(x), ...)
}

tri <- raster_to_trimesh(nt_raster)
wire3d(tri)
rglwidget()

(This doesn't allow constraining to a boundary, though the two edge-sets can be merged together and re-triangulated in order to do that - I'm not totally au fait with that work flow yet).

It turns out I added it to the scene that already had uluru_trimesh wire3d, so it was a good comparison. It's great to have examples like this that are readily convertible, primitives make these changes quite simple, if a bit abstract!

[mdsumner]

I've also extended from that example above to obtain a visible image, and do the (many!) necessary conversions and tweaks to texture map that on.

First up, I would definitely crop nt_raster before starting here, otherwise it's a bit big if using rglwidget - though it did work and it looks good.

## there are better ways to get imagery, but this is a good old faithful detault 
im <- dismo::gmap(nt_raster, type = "satellite", scale = 2)

## so this is a bit head-stretching, because multiple coordinate systems in play

## we want the [0,1,0,1] coordinates of the image in terms of the geographic mesh
## 1) mesh is lcc, backwards to longlat WGS84
xyl <- rgdal::project(t(tri$vb[1:2, ]), projection(nt_raster), inv = TRUE)
## 2) image is merc, forwards to that
xym <- rgdal::project(xyl, projection(im))
## 3) the cell (in [0,1,0,1]) for our mesh coordinates
cell <- cellFromXY(im, xym)
## 4) use that cell to get native from a rescaled im
xyim <- xyFromCell(setExtent(im, extent(0, 1, 0, 1)), cell)

## 5) convert to RGB from palette (might be a raster fun for this...)
imrgb <- setValues(brick(im, im, im), t(col2rgb(im@legend@colortable[values(im) + 1])))

## 6) write to PNG (that's the only way we can texture map)
texfile <- sprintf("%s.png", tempfile())
rgdal::writeGDAL(as(imrgb, "SpatialGridDataFrame"), texfile, driver = "PNG")

## clear plot if need be
## rgl.clear()

## finally, plot the mesh as filled triangles and specify the texture mapping
## (it has to not be col = "black", which is the default )
shade3d(tri, texcoords = xyim[tri$it, ], texture = texfile, col = "white")
rglwidget()

Finally, I think steps 3,4 would be better if we interpolated the position rather than just the cell lookup, since it might be non-aligned with the cell (like your triangulation in the boundary). For that we can put the x/y values in copies of the raster and extract(, method = "bilinear") or something but I can't quite manage it right now.

[MilesMcBain]

:scream: :clap: :clap: :clap:

This is freaking awesome! Thanks Mike!