from Miles:

In "What makes polygons different to lines?" None of your arguments for segments made much sense > to me, since they contained a lot of domain knowledge. I think that's fine given the package audience, > but you might consider rewording a little for accessibility.

Probably remove this function, or at least split existing functionality into

• primitives_triangle, primitives_segment, primitives_point
  or
• primitives2d, primitives1d, primitives0d

See impetus here: #46

Ultimately we need to be able to have a model that has all of these things, which is what anglr might ultimately be able to do, a hyper-mapview. rgl has these mixed forms to a small degree, and we'd want to have lines and polygons in a model, sharing the same vertex pool - i.e. lines that follow terrain, and surfaces from contours. Can we sensibly have an everything-model that one function can build from any input, but also the family of functions that can add smaller parts to a larger model probably includes the functions above.

http://www.quantdec.com/misc/Avencia.pdf

first stab at mesh3d:

library(rgl)

x <- tetrahedron3d()
class(x)
str(x)
library(tibble)
fun <- function(x) {
  v <- setNames(as_tibble(t(x$vb[1:3, ])), c("x_", "y_", "z_"))
  v$vertex_ <- spbabel:::id_n(nrow(v))
  tXv <- tibble(vertex_ = v$vertex_[x$it], 
                triangle_ = spbabel:::id_n(ncol(x$it))[rep(seq(ncol(x$it)), each = 3)])
  tt <- tibble(triangle_ = tXv$triangle_[seq(ncol(x$it))], object_ = "1")
  o <- tibble(object_ = "1")
  structure(list(o = o, t = tt, tXv = tXv, v = v), class = "trimesh")
}

(Doesn't work)

i.e. open the device or check for one, and that it's not the null device

It should at least replace the source projection with the geocent string used.

And . . . the meta table could record every transformation made, and maybe keep track of what the current default attributes are . . . and their units . . .

These files choke at pslg:

SSRUs_East_Ant_Polys.shp

cst10_polygon.shp

Is the dedupe failing?

@kent37 thanks for the data, I've included it within the package silicate, which I hope is ok? There it is in PATH form, which you can create from an sf with silicate::PATH( ). There's still another step into anglr form while I figure out the best way forward.

If you are feeling brave, try installing anglr, there is one other Github package silicate, and rgl and a few other packages are required from CRAN. You'll likely need a pretty recent version of R, I'm using 3.4.1.

The rgl.clear() at the start of the plot and rglwidget() is a bit of overkill, but it's necessary for me working on RStudio Server.

This project is still very early and raw, but it's looking promising. Let me know how it goes, good or bad - it's helpful to know any feedback. We can put other features in the scene easily enough, including polygons with topography relief and so forth, but all of that takes a bit of knowledge of the structures at the moment.

## devtools::install_github("hypertidy/silicate")
## devtools::install_github("hypertidy/anglr")

library(anglr)
library(rgl)
data("flight_tracks", package = "silicate")
## convert to topological graph
aa <- anglr(flight_tracks)

rgl.clear()
## sets up a basic plot, with a colour for each feature
##  (M is ignored, but it's in the model - see aa$v)
plot(aa)
## this is a hack, z-metres in lon-lat
## play with the third value by orders of magnitude
## but re-run the rgl.clear to rglwidget funs 
## completely start again if needed
aspect3d(1, 1, 0.0001)
## this call is not needed unless you are on RSTudio Server, or have no OpenGL rgl device
rglwidget()

Finally, a specific thing if you can let me know, are you able to rotate the scene as well as zoom in and out? I'm having troubles, I don't understand the details of rglwidget yet. :)

http://blog.thematicmapping.org/2016/09/creating-tins-in-saga-gis.html

There's an example somewhere I used terrain() to do a simplistic curvature filter.

ah, no it's not building that is slow - quads are slow in rglwidget, why is that? Also noted back here:

http://mdsumner.github.io/2016/01/20/ManifoldRaster.html

There seems to be a limit in the size that can be handled

library(marmap)
library(raster)
library(rangl)
data("hawaii", package = "marmap")
r <- as.raster(hawaii)
p <- rangl(aggregate(r, fact = 8))
xy <- rgdal::project(as.matrix(p$v[, 1:2]), "+proj=laea +lon_0=180 +lat_0=-0")
p$v$y_ <- xy[,2]
p$v$x_ <- xy[,1]
plot(p)

There's also some issue with caching in RSTudio I think, the widgets are piling up maybe ... it's important to reproject the vertices above to remove any uncertainty about the aspect ratio in the plot.

https://github.com/mdsumner/gris/blob/master/R/trianglePrimitives.R

Perhaps add this id to the object, or maybe this is just better done by putting material properties on rangl always

is this just a wave for CRAN, or do we need to add more information to the pkg?

NOTE
Package has a FOSS license but eventually depends on the following
package which restricts use:
  RTriangle

It's unlikely that anyone is going to replace RTriangle::triangulate.

Could use manifoldr for those systems.

Having many many lines is less efficient than it could be. But plotting segments is fast since you only need one plot call.

Is it also much faster in base lines() the same as in rgl? What about ggplot2?

• remove object-loop for creating primitives in rangl.SpatialLines
• always create primitives for lines and polygons, consider removing the branch model completely
• develop examples based on particles, time series, tracks

Examples

See AMPS particle trajectories . . .

See geolife example

Dated contours from time series rasters

library(raadtools)
icefile_db <- subset(icefiles(), 
                     date > as.POSIXct("2015-12-15") & date < as.POSIXct("2016-03-01"))
dummy <- readice()
pex <- extent(60, 100, -70, -55)
xyl <- projectExtent(raster(pex, crs = "+proj=longlat +ellps=WGS84"), projection(dummy))

icestack <- readice(icefile_db$date, hemisphere = "south", setNA = FALSE, 
                    xylim = xyl)
library(dplyr)
library(spbabel)
l <- vector('list', nlayers(icestack))
for (i in seq(nlayers(icestack))) {
  cl <- spbabel::sptable(spTransform(rasterToContour(icestack[[i]], lev = 15), CRS(projection(xyl))))
  cl <- cl %>% filter(branch_ == 1, object_ == 1) %>% mutate(object_ = as.character(i)) 
 l[[i]] <- cl

}

ice <-  spbabel::sp(bind_rows(l), tibble(date = icefile_db$date), crs = projection(icestack))

Geozoo has a bunch of interesting geometries, but not many have edges. You can use ball pivoting from vcg to some extent.

Good grist for rangl.

I downloaded a shapefile which sf read in as a MULTILINESTRING. anglr converts this to a linemesh, which I don't understand how to draw in 3D. What am I missing?

Edit: I've committed the shape file and related code in alexander_clark_park in our vr_vs_sp repo.

include ability to "copy down" object attributes to vertices, for e.g. contour lines

need to think about how maintaining polygon branches sits with decomposing lines generally (probably always do it?)

Is it worth using triangles as per polyggon to do this?

http://gis.stackexchange.com/questions/224035/how-to-create-a-crisp-topographical-terrain-map-with-ggplot2

title: "Structures for hierarchical data in R"
author: "Michael Sumner"
date: "14 December 2016"
output: html_document

knitr::opts_chunk$set(echo = TRUE)

Hierarchical data is generally hard to deal with since the layout doesn't often match the way tools must then work with the data. Physical structures are used to match the logical structure in data sets, and this is not necessarily an efficient or flexible way to deal with the data itself.

** Example of extracting stuff from a nested list**

library(repurrrsive)
got_chars

Tables

The tidy framework in R generally dictates that the long form of tabular storage is the best way to work with statistical models and graphics.

The tools for converting from untidy nested structures, to long-form tables is well established, but the capacity to deal with
multiple related tables is more difficult.

What does nesting represent? Nesting in a tree-like structure represents data as nodes in a graph, with the parent-node and child-nodes as a straightforward encoding of relations.

Tree example of a nested list

These relations can be encoded in "normal-form", where the proper entities in the data structure are stored each
in a table, and the tree is not stored explicitly but as 1-many and many-1 links in relations.

** Table example of a nested list **

When the nodes in the graph are themselves tidy, it's very easy to convert from the nested recursive structure to a set of tables.

** Bag of relations **

A list of tables doesn't seem especially helpful, and it's not - at the user-level - but for a developer it's a very straightforward data structure to work with.

Tools to build a map of nested structure

Applying the map to extract tidy tables

Collapsing to entity tables

library(repurrrsive)

library(purrr)
library(tibble)
library(dplyr)
library(tidyr)
map_df(wesanderson, function(x) tibble(color = x), .id = "palette")

map(got_chars[[1]], length)


library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf"))
map(st_geometry(nc), length)

nc %>% select(CNTY_ID, AREA, geometry) %>% mutate(len = at_depth(geometry, 1, length)) %>% filter(len > 1)

shp <- as(nc, "Spatial")
library(spbabel)
library(sp)
map_table(nc)

## break down nc and rebuild as tibble
tnc <- nc %>% select(-geometry) %>% as_tibble()
tnc$geometry <- st_geometry(nc)
st_as_sf(tnc) %>% by_row(map_table)
map_table.tbl_df <- function(x, ...) {
  map_table(st_as_sf(x))
}

x <- st_as_sf(tnc) %>% by_row(map_table) %>% select(.out)

Examples with spatial data

** Normal form for simple features **

** Primitives form (topology) **

** Recording aesthetics as variables with primitives **

** Recording variables generally on entities **

** Primitives from ggvis **

Try example with trip or similar, an auto rgl animation akin to the flying light ball example.

Is it worth piggy-backing on tourr at all?

This doesn't work, a certain combination of polygons and max_area - the linking with vindex$row_n is getting out of whack somehow, is it the order of the verts or is there somehow a too-small triangle throwing it out?

## problem
p <- sp::spTransform(wrld_simpl, "+proj=laea +ellps=WGS84")

edge <- rangl(p[seq(2, nrow(p), by= 8), ], max_area = 5e8)
plot(edges$v$x_, edges$v$y_, type = "n")
## the rangl tables are ... not fine
## because this plot is ... not fine
lapply(split(edge$tXv, edge$tXv$triangle_), function(x) lines((inner_join(x, edge$v) )[, c("x_", "y_")], col = "firebrick"))
#Joining, by = "vertex_"
#Error: not compatible with requested type 

## the rgl constructs are now bogus
## (the plot is true, the indexes are out of whack some how)
a <- plot(edge)
g <- plot(globe(edge))
ncol(a$vb)
ncol(a$it)

Probably this means we just have one function "mesh" - or "rangl()"?

Build notes

• Travis needs libproj-dev installed, via .travis.yml (CRAN does not)
• R-hub needs SystemRequirements: PROJ.4
• All need export RGL_USE_NULL=TRUE
• plot in examples is fine, but there's a warning from R CMD check in interactive mode (??)

BGM needs

• mixed quads and triangles (use extrude3d though later use RTriangle)
• bgm data structure probably should be rangl from the start

I've seen this somewhere else too

cl <- contourLines(volcano)
d <- dplyr::bind_rows(lapply(cl, as_tibble))
library(dplyr)
library(spbabel)
x <- sp(transmute(d, x_ = x, y_ = y, branch_ = level, object_ = 1, order_ = seq(nrow(d))))

Warning: stack imbalance in '.Call', 168 then 167
Warning: stack imbalance in '<-', 166 then 165
Warning: stack imbalance in '{', 163 then 162
Warning: stack imbalance in '.Call', 167 then 166
Warning: stack imbalance in '<-', 165 then 164
Warning: stack imbalance in '{', 162 then 161
Warning: stack imbalance in '.Call', 166 then 165
Warning: stack imbalance in '<-', 164 then 163
Warning: stack imbalance in '{', 161 then 160
Warning: stack imbalance in '.Call', 165 then 164
Warning: stack imbalance in '<-', 163 then 162
Warning: stack imbalance in '{', 160 then 159

This is a simple d admin map with 2 very detailed polygons with upwards of 100K and 200K vertices each, due to a baked-in coastline. The 14th polygon is fine, and comes out of rangl fine - the 17th takes a lot of memory (more than 20Gb) at last check to run the constrained Delaunay algo (without a max area value). This is not a reasonable test

https://github.com/mikedh/trimesh

http://artax.karlin.mff.cuni.cz/r-help/library/geometry/html/trimesh.html

also matlab uses this name

http://blog.ryanwalker.us/2014/03/going-viral-with-rs-igraph-package.html?m=1

Something like this, index is not right yet.

Also, quad_mesh is using a structural index, which is baad.

quadToTriangle <- function(x) {
  v <- x$v
  v$vertex_ <- seq(nrow(v))
  meta <- x$meta
  tab <- x$qXv
  n4 <- nrow(tab) / 4L
  tXv <- tibble(vertex_ = rep(c(1, 2, 3, 2, 3, 4), n4) + rep(seq(0, length = n4, by = 3), 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
}

Why doesn't $t[7, ] get filtered out by over?

library(sf)
library(dplyr)
library(rangl)
vpal <- function(n) {
  if (!is.null(dim(n))) {
    n <- nrow(n)
  }
  viridis::viridis(n)
}
p_self <- st_sfc(st_polygon(list(cbind(c(0, 0, 1, -0.1, -0.1, 1.5, 0), c(0, 1, 1, 0.2, 1.2, 1.2, 0)))))
ri <- rangl(as(p_self, "Spatial"))

p <- do.call(cbind, structure(list(x = c(-0.205701357466063, -0.0357229982293921, 
                                         -0.0697186700767263, 0.312732638205784, 0.4855439700964, 0.343895337399174, 
                                         0.471379106826677, 1.62156600432815), y = c(0.424072398190045, 
                                                                                     0.160605941373205, 0.336250245917765, 1.09831988982884, 0.54022427700177, 
                                                                                     0.732866417469998, -0.0801967342120795, 1.71024198308086)), .Names = c("x", 
                                                                                                                                                            "y")))
library(sp)
inp <- over(SpatialPoints(p), as(p_self, "Spatial"))
ploti <- function(x, col = NULL) {
  xx <- x %>% inner_join(pi$tXv) %>% inner_join(pi$v) %>% select(x_, y_) 
  #print(xx)
  polypath(xx, col = col)
  #$print(apply(as.matrix(xx), 2, mean))
  text(matrix(apply(as.matrix(xx), 2, mean), ncol = 2), lab = i)
}
plot(p)
plot(p_self, add = TRUE)
text(p, lab = sprintf("pt%i", seq_len(nrow(p))))
for (i in seq_len(nrow(pi$t))) ploti(pi$t[i, ], col =  scales::alpha("red", i/10))
data.frame(pt = seq_len(nrow(p)), pip = over(SpatialPoints(p), as(p_self, "Spatial")))

    pt pip
1  1  NA
2  2  NA
3  3   1
4  4   1
5  5   1
6  6  NA
7  7  NA
8  8  NA    

break and re-mesh by combining two together

We need to differentiate

• continuous: the vertex "v" table
• discrete: the primitives or the object "o" table

Examples are the colour of a country in countriesHigh (discrete) vs. the colour scale across triangles (continuous).

Are there lessons from ggplot2 or scales for this?

Could we have "fill_cts" and "fill_dsc" args, mutually exclusive?

Copy down attributes from objects

• to vertices, they are unique
• to branches/primitives, potentially splitting what were shared (x.y) vertices

This allows columnar polygons or ribbon lines

this might be spbabel, but the indexes are extremely heavy - and this is really useful but - perhaps they just get created when merging or subsetting operations are happening, and we always resort to a counting index (though why have it at all ...)

We need

• RTriangle triangulate for sf types, especially MULTIPOINTZ
• questions regarding triangulation within features, so this is more like the model of not-sharing vertices to connect multi-point geoms (it fits rangl but hasn't really come up before)
• etc.

here's a raw example

library(lidR)
library(sf)
LASfile <- system.file("extdata", "Topography.laz", package="lidR")
lidar <- readLAS(LASfile)
ground <- lasfilterground(lidar)
library(dplyr)
pts <- local({
points <- st_sfc(st_multipoint(as.matrix(ground@data[, c("X", "Y", "Z")])))
d <- data.frame(file = "Topography.laz", method = "filterground")
d[["geometry"]] <- points
st_as_sf(d)
})
#sp_tri <- st_triangulate(pts)

rg_triangulate <- function(x, ...) UseMethod("rg_triangulate")

rg_triangulate.MULTIPOINT <- function(x, ...) {
    m <- as.matrix(x)
    x <- RTriangle::triangulate(RTriangle::pslg(P = m[, 1:2, drop = FALSE], PA = m[,3,drop = FALSE]), ...)
    a <- rgl::tetrahedron3d()
    a$vb <- t(cbind(x$P, x$PA[, 1, drop = FALSE], 1))
    a$it <- t(x$T)
    a
}

shade3d(rg_triangulate(pts$geometry[[1]], a = 5), col = "white")
points3d(as.matrix(lidar@data[, c("X", "Y", "Z")]))

#46

these ones: https://hypertidy.github.io/anglr/articles/spatial-normal-forms.html

Hawaii example

library(rangl)
library(raster)
library(marmap) ## has topo data
data(hawaii)
library(rgeos)
haw <- setExtent(raster(unclass(hawaii)), extent(unlist(lapply(attr(hawaii, "dimnames"), function(x) as.numeric(c(head(x, 1), tail(x, 1)))))))
projection(haw) <- "+proj=longlat +ellps=WGS84"

hpoly <- rasterToPolygons(cut(haw, c(cellStats(haw, min), 0, cellStats(haw, max))), dissolve = TRUE)

## just some polygons
hpoly <- rgeos::gPolygonize(rasterToContour(cut(haw, c(cellStats(haw, min), 0, cellStats(haw, max))), lev = 1.5))

plot(globe(rangl(haw * 280)))
rh <- rangl(hpoly)
rh$v$z_ <- 4e6
plot(globe(rh))

• TIN from raster http://blog.thematicmapping.org/2016/09/creating-tins-in-saga-gis.html
• Whuber's manifesto: http://www.quantdec.com/misc/Avencia.pdf
• igraph http://blog.ryanwalker.us/2014/03/going-viral-with-rs-igraph-package.html?m=1

See rgl examples in mesh, the colours and other stuff can be maintained and carried through. (Joy to the world).

need functions for raster vertexIndex_FromCell, or

ll_FromCell(raster, cells) 
lr_FromCell/tl_FromCell/tr_FromCell

where the four indices are implicit on the transform, and sparsely generated only for populated cells. Then generalize for L3bin, to tie together tabularaster and quadmesh.

Does this work by

• translating centre cell index to ll/lr/tl/tr index
• then just doing mod and %/% x, y to get the values

So the full set of coords is never generated, maybe never generated until transformed or visualized by rgl?

https://www.linkedin.com/pulse/hermit-crab-gets-3d-printed-designer-shell-modeled-using-jose-i-rey?trk=prof-post

It would be neat to have "selected" behaviour, say to switch discrete / continuous aesthetics as per #3

A brushing mode could switch between these.

http://gmsh.info/doc/texinfo/gmsh.html

• rgl needs ascii read for STL, see demos/sphere-surf.stl
• explore .geo format

Everything is filtered out as being a hole, why? What about sc?

library(sf)
f <- "maps.nyc.gov/download/3dmodel/DA_WISE_Multipatch/DA_WISE_Multipatch/DA10_3D_Buildings_Multipatch.gdb"
d <- read_sf(f)

r <- rangl(d[1:10, ])

Sigh

sos::ff <- findFn("mesh")

AT.FLUKA.read.USRBIN.mesh
closemeshKD
distmesh2d
distmeshnd
mesh
mesh.2D.rectangular
mesh.2D.simple
mesh.dcircle
mesh.diff
mesh.drectangle
mesh.dsphere
mesh.hunif
mesh2grey
mesh2ply
mesh3d
meshcube
meshDist
meshgrid
meshintegrity
meshPlaneIntersect
meshres
mvmesh-geom
mvmesh-methods
mvmesh-package
mvmeshmisc
ply2mesh
quad2trimesh
quadmesh
ray2mesh
rmvmesh
rotmesh.onto
scalemesh
submesh.grid
submesh.mesh
tetramesh
theme_mesh
tri.mesh
trimesh
vcgUniformRemesh