mesh forms for mapdeck uber.gl

[mdsumner]

From @mdsumner on August 31, 2018 2:21

Here is my attempt to get XYZ triangles into mapdeck, so we have a shared example to run with.

key <- Sys.getenv("MAPBOX_TOKEN")

elevate_triangles <- function(g, r) {
    ## assuming geometry of POLYGON
    sf::st_polygon(list(cbind(g[[1]], raster::extract(r, g[[1]]))))
}

elevate_sf <- function(x, r, exag = 1) {
    sf::st_set_geometry(x, sf::st_sfc(purrr::map(sf::st_geometry(x), elevate_triangles, r = r * exag), crs = sf::st_crs(x)))
}

## get sf data (in longlat)
nc <- sf::read_sf(system.file("shape/nc.shp", package="sf"))
bb <- sf::st_bbox(nc)
bb <- bb + c(-2, -2, 2, 2)
## get topo data (also longlat)
topo <- marmap::as.raster(marmap::getNOAA.bathy(bb["xmax"], bb["xmin"], bb["ymax"], bb["ymin"], 
                                                                                                resolution = 2))
## sf POLYGON triangle for every element (a is *area* in native coordinates )
triangles <- sf::st_cast(sfdct::ct_triangulate(nc, a = .1))

library(mapdeck)
zfac <- 1
x <- elevate_sf(triangles, topo * zfac)

## my attempt at mapdeck XYZ triangle glory
mapdeck(token = key, style = 'mapbox://styles/mapbox/dark-v9', location = c(mean(bb[c("xmin", "xmax")]), mean(bb[c("ymin", "ymax")])), zoom = 5) %>% 
        add_polygon(data = x, layer_id = "polylayer") 

I use exaggeration (zfac) to make sure it's not just a scaling issue. So far all POLYGONs are totally flat!

Here's how I go about it with silicate/anglr (if installation gives you grief just let me know, It might not be setting up all deps properly).

## same but done with rgl
library(anglr)
plot3d(TRI(x))  ## triangulating is redundant here, but works with any sf polygon layer (use DEL for better/denser triangles with "max_area" arg 
rgl::aspect3d(1, 1, 1/10)
## if you in RStudio Server the output needs a widget (it's slow, much faster with local rgl device)
rgl::rglwidget()  ## also rerun this if aspect3d is updated ...

Copied from original issue: SymbolixAU/mapdeck#40

[mdsumner]

Works now! More soon :+1:

[mdsumner]

Here's a more direct example from a raster.

library(quadmesh)
data("hawaii", package = "marmap")
qm <- hawaii %>% 
  marmap::as.raster() %>% 
  raster::aggregate(fact = 2) %>% 
  quadmesh()

## z exag
qm$vb[3, ] <- qm$vb[3, ] * 20
qm2sf <- function(x, crs = NA_character_) {
  tri <- if(is.null(x$it)) quadmesh::triangulate_quads(x$ib) else x$it
  psfc <- sf::st_sfc(lapply(split(rbind(tri, tri[1, ]),
                          rep(seq_len(ncol(tri)), each = 4))
                    , 
                    function(a) sf::st_polygon(list(t(qm$vb[, a])))), crs = crs)
  sf::st_sf(geometry = psfc, a = .colMeans(matrix(x$vb[3, tri], 3), 3, ncol(tri)))
}
p <- qm2sf(qm)
library(mapdeck)
mapdeck(token = Sys.getenv("MAPBOX_TOKEN"), style = 'mapbox://styles/mapbox/dark-v9', 
        location = c(mean(qm$vb[1, ]), mean(qm$vb[2, ])), zoom = 5) %>% 
  add_polygon(data = p, layer_id = "polylayer", fill_colour = "a" ) 

Keen to discuss ways to get these mesh data in without having to construct sf (but also keen to discuss fast sf construction).

Impressive work!

[mdsumner]

I'll move this now, it's just me experimenting really - given that Z-from-sf now works. 💯

[mpadge]

That's great! @layik is currently adapting geoplumber to use deck.gl via mapdeck instead of leaflet, so he'll soon have loads of insight. That work has, however, just begun, and he'll likely be wrestling for .. no idea, but a while.

[mdsumner]

Another example, (because no tetures) so we can put colours explicitly onto polygons

elevate_triangles <- function(g, r) {
  ## assuming geometry of POLYGON
  sf::st_polygon(list(cbind(g[[1]], raster::extract(r, g[[1]]))))
}

elevate_sf <- function(x, r, exag = 1) {
  sf::st_set_geometry(x, sf::st_sfc(purrr::map(sf::st_geometry(x), 
                      elevate_triangles, r = r * exag), crs = sf::st_crs(x)))
}

e <- new("Extent", xmin = -400389.695476286, xmax = -394710.48041841, 
         ymin = -45598.0485452332, ymax = -41690.5544242159)
d <- raster("../gis_vs_web3D/data/ELVIS_CLIP.tif")
p <- spex::polygonize(crop(d, e))
library(ceramic)
im <- cc_location(cbind(131, -25.35))
dz <- elevate_sf(st_transform(p, projection(d)), d, exag = 2)
dz <- st_transform(dz, 4326)
#dz$my_colour <- colourvalues::color_values(extract(d, st_coordinates(st_centroid(dz))[, 1:2]))
v <- extract(im, rgdal::project(st_coordinates(st_centroid(dz))[, 1:2], projection(im)))
v[is.na(v)] <- 0
dz$my_colour <- rgb(v[,1], v[,2], v[,3], maxColorValue = 255, alpha = 255)
mapdeck( ) %>%
  add_sf(
    data = dz
    , fill_colour = "my_colour", update_view = TRUE
  )

Still needs direct mesh to feed mapdeck, but it's so fast we can push this pretty far.

[SymbolixAU]

I like where this is going!

[mdsumner]

Here's another example, I've kept it pretty coarse because I'm working only within the browser. It took quite a bit of code debugging but it's a more straightforward example than above. (The image overlay is bad, but because it's not a texture, just a triangle colouring)

##remotes::install_github("hypertidy/quadmesh")
library(quadmesh)
##remotes::install_github("hypertidy/ceramic")
library(ceramic)
library(sf)
library(mapdeck)

## define a region by polygon
sd <- rnaturalearth::ne_states(returnclass = "sf") %>%
  dplyr::filter(name == "Tasmania")

## obtain imagery and elevation from Mapbox
im <- cc_location(st_coordinates(st_centroid(sd)), buffer = 3e5, zoom = 7,  debug = TRUE)
el <- cc_elevation(st_coordinates(st_centroid(sd)), buffer = 4e5, zoom = 5)

## convert polygon to mercator and mask the mapbox imagery
prsd <- st_transform(sd, raster::projection(im))
im <- raster::mask(im, prsd)
el <- raster::mask(el, prsd)

## convert to rgl mesh form (triangles)
tm  <- triangmesh(el)

## convert triangle mesh to sf Z triangles, removing NA vertices/primitives
tm_to_sf <- function(x, exag = 1) {
  ## convert bad z verts to bad triangles
  ibad <- which(is.na(x$vb[3, ]))
  ii <- x$it
  ii[ii %in% ibad] <- NA
  bad <- is.na(colSums(ii))
  x$vb[3, ibad] <- min(x$vb[3, ], na.rm = TRUE)
  x$vb[3, ] <- x$vb[3, ] * exag
  sfc <- vector("list", ncol(x$it))
  ## template sfg for update, so we don't check the closing coordinate every time
  sfg <- sf::st_polygon(list(cbind(c(1, 1, 2, 1), c(1, 2, 2, 1), 0)))
  for (i in seq_along(sfc)) {
    sfg[[1]] <- t(x$vb[1:3,x$it[c(1L, 2L, 3L, 1L),i]])
    sfc[[i]] <- sfg
  }
  sf::st_sf(a = seq_len(i), geometry = sf::st_sfc(sfc))[!bad, ]
}
z <- tm_to_sf(tm, exag = 10)
z <- st_set_crs(z, st_crs(im))

## now map colours on
v <- raster::extract(im, st_coordinates(st_centroid(z))[, 1:2])
v[is.na(v)] <- 0
z$my_colour <- rgb(v[,1], v[,2], v[,3], maxColorValue = 255, alpha = 255)
library(mapdeck)

## convert back to longlat and send to mapdeck! 
zll <- st_transform(z, 4326)
mapdeck( ) %>%
  add_sf(
    data = zll
    , fill_colour = "my_colour", update_view = TRUE
  )

[SymbolixAU]

I see your Tasmania, and I raise you Victoria

[mdsumner]

It's on!

[mdsumner]

This is covered by anglr::as.mesh3d -> mapdeck