hypertidy / anglr

Mesh creation and topology for spatial data (and not just geographic)
https://hypertidy.github.io/anglr/
83 stars 10 forks source link

add walls for `copy_down(x, "colname")` #102

Open mdsumner opened 4 years ago

mdsumner commented 4 years ago

probably use sc and quads for paths, but for the sides of a matrix it might be easiest to consider the margin and the corners as a polygon to triangulate

mdsumner commented 4 years ago

Walls: https://github.com/hypertidy/anglr/wiki/Examples#walls

mdsumner commented 4 years ago

done with quads on a tmesh: https://gist.github.com/mdsumner/da52ff4b476f25595c31c9341a90d568

mdsumner commented 4 years ago

A bit more of a go, you need to mesh each polygon separately, so it's easier to do on each sf row (but you could loop over the TRI0 indexes)

library(anglr)
library(ozmaps)

sfx <- rmapshaper::ms_simplify(abs_lga)

sfx$z <- rpois(nrow(sfx), 2)

tri <-copy_down( TRI0(sfx[1, ]), "z")

mesh0 <- as.mesh3d(tri, meshColor = "faces")  

add_walls <- function(mesh, Zfloor = 0) {                

## find non-recurring edges (they are the boundary)
edges <- t(cbind(mesh$it[1:2, ], mesh$it[2:3, ], mesh$it[c(3, 1), ]))
native <- edges[,1] < edges[,2]
## by sorting each row
edges1 <- cbind(pmin(edges[,1], edges[,2]), pmax(edges[,1], edges[,2]))

## and find duplicates, keep only those that occur once (they are mesh boundary)
tib <- tibble::tibble(.vx0 = edges1[,1], .vx1 = edges1[,2])
tib <- dplyr::ungroup(dplyr::mutate(dplyr::group_by(tib, .vx0, .vx1), n = dplyr::n() ))
tib <- dplyr::filter(tib, n < 2)
boundary <- as.matrix(tib[c(".vx0", ".vx1")])

## add the quads with a constant z
nv <- rbind(mesh$vb[1:2, t(boundary[,2:1])], z = Zfloor, h = 1)
## or a relative z
#nv <- rbind(mesh$vb[1:2, t(boundary[,2:1])], z = mesh$vb[3, rep(boundary[,1], each = 2)] - 30, h = 1)
ncurr <- dim(mesh$vb)[2L]
nnew <- nrow(boundary) * 2
quads <- cbind(boundary, matrix(ncurr + seq_len(nnew), ncol = 2, byrow = TRUE))
mesh$ib <- t(quads)            ## the quads are all new
vb <- cbind(mesh$vb, nv)  ## but the geometry is part new
mesh$vb <- vb 

mesh$material$color <- c("#A020F0","#FFFFFF")[(mesh$vb[3, ] > 0) + 1]

mesh 
}

mesh_list <- lapply(seq_len(nrow(sfx)), function(i) {

  tri <-copy_down( TRI0(sfx[i, ]), "z")
  tri$object$color_ <- "grey"
  mesh0 <- as.mesh3d(tri, meshColor = "vertices")  

  add_walls(mesh0)  
})

library(sf)
plane <- as.mesh3d(st_sfc(st_polygon(list(cbind(c(0, 3600, 3600, 0, 0), 
                                             c(-90, -90, 9000, 9000, -90))))))
clear3d()
lapply(mesh_list, plot3d, add = TRUE)
# par3d(ignoreExtent = T)
# plot3d(plane, add = TRUE, color = "darkgrey")
bg3d("#333333")

rgl::aspect3d(1, 1, .05)

image