In hypertidy/rangl: Mesh Topology and Visualization for Spatial Data

source("setup.R")
knitr::opts_chunk$set(rgl.newwindow = TRUE)
set.seed(123)

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>")
rgl::par3d(windowRect = c(100, 100, 512 + 100, 512 +100))
library(raster)

Introduction

The 'anglr' package illustrates some generalizations of GIS-y and topology tasks in R with "tables". See the package project for more information.

Get some maps and plot in 3D - in plane view, or globe view.

NOTE: these plots are interactive, but performance and quality will be better (for now) if run locally.

Set up a simple world countries data set.

library(rgl)
library(raster)
library(sf)
## convert to triangles and plot
library(anglr)
library(silicate)

simpleworld <- rgeos::gBuffer(simpleworld, width = 0, byid = TRUE)

cmesh <- DEL(simpleworld)
plot(cmesh)
plot3d(cmesh)
rglwidget()

rgl::rgl.clear()

sids <- sf::read_sf(system.file("shape/nc.shp", package="sf"))

ex <- extent(sf::st_bbox(sids)[c(1, 3, 2, 4)]) + 5
gl <- graticule::graticule(seq(xmin(ex), xmax(ex), length = 15), 
                           seq(ymin(ex), ymax(ex), length = 8))


## convert to triangles, but wrap onto globe then plot
smesh <- DEL(sids)
plot3d(globe(smesh))
mgl <- SC(gl)
mgl$o$color_ <- "black"
plot3d(globe(mgl), lwd = 2, add = TRUE)
rgl::rglwidget()

Holes are (trivially) supported.

It's trivial to have holes, because there are no holes, we have a true surface, composed of 2D primitives (triangles).

rgl::rgl.clear()

library(spbabel)
data(holey)

## SpatialPolygonsDataFrame
sph <- sp(holey)

glh <- TRI(sph)
plot(glh)
plot3d(glh)
rgl::rglwidget()

Lines are supported.

rgl::rgl.clear()
library(dplyr)
library(sp)
linehouse <- as(sph, "SpatialLinesDataFrame")
plot3d(PATH(linehouse))
rgl::rglwidget()

Globe lines

rgl::rgl.clear()

lmesh <- SC(as(simpleworld, "SpatialLinesDataFrame"))
plot3d(globe(lmesh))
rgl::rglwidget()

rgl objects

Rgl mesh3d objects that use "triangle" primitives are supported.

(This is meant to be deprecated but we'll need a bit more work to create silicate models from rgl).

rgl::rgl.clear()

dod <- anglr(dodecahedron3d(col = "cyan"))
octo <- anglr(translate3d(octahedron3d(col = "blue"), 6, 0, 0))
plot(dod, col = viridis::viridis(5)[1], alpha = 0.3)
plot(octo, col = viridis::viridis(5)[5], alpha = 0.3)
bg3d("grey")
rgl::rglwidget()

To complete the support for these rgl objects we need quads, and to allows a mix of quads and triangles in one data set (that's what extrude3d uses). Extrusion is a bit limiting so it's unclear how useful this is (yes it is common, though). See rgl::extrude3d for the most readily available method in R.

Points

And points work! (Don't laugh, actually they don't work atm).

rgl::rgl.clear()

library(anglr)

mpts <- as(as(simpleworld, "SpatialLinesDataFrame"), "SpatialMultiPointsDataFrame")
plot3d(ARC(mpts))
rgl::view3d(theta = 25, phi = 3)
rgl::rglwidget()

Trips

The trip package includes a 'walrus818' data set courtesy of Anthony Fischbach. Zoom around and see if you can find them.

Here we also use the finite element aspect of our topology representation. We can forget about complicated algorithms to cut and splice lines and just throw away triangles south of a chosen latitude. (Still needs work, but we've proved it works well enough. )

Needs work:

rgl::rgl.clear()
library(trip)
library(anglr)
data(walrus818)

library(graticule)
prj <-"+proj=laea +lon_0=0 +lat_0=90 +ellps=WGS84"
gr <- graticule(lats = seq(40, 85, by = 5), ylim = c(35, 89.5), proj = prj)


w <- spTransform(subset(simpleworld, coordinates(simpleworld)[,2] > 10),  prj)
walrus <- spTransform(walrus818, prj)

gr$color_ <- "black"
rgl::par3d(windowRect = c(100, 100, 912 + 100, 912 +100))
plot3d(SC(gr))
w$color_ <- sample(viridis::inferno(nrow(w)))
wmap <- TRI(w)
tXv <- wmap$tXv %>% dplyr::inner_join(wmap$v)

keep <- tibble::tibble(triangle_ = unique(tXv$triangle_[rgdal::project(as.matrix(tXv[c("x_", "y_")]), prj, inv = TRUE)[,2] > 35]))
wmap$t <- wmap$t %>% dplyr::inner_join(keep)
## trim by primitives

## we mashed the graph and bit but it works, later we
## can propagate the filter ...
plot3d(wmap, specular = "black", add = TRUE)
plot3d(SC(walrus), add = TRUE)
um <- structure(c(0.934230506420135, 0.343760699033737, 0.0950899347662926, 
                  0, -0.302941381931305, 0.905495941638947, -0.297159105539322, 
                  0, -0.188255190849304, 0.24880850315094, 0.950081348419189, 0, 
                  0, 0, 0, 1), .Dim = c(4L, 4L))
par3d(userMatrix = um)
rgl::rglwidget()

A polygon layer may be overlaid with a raster layer using the z argument.

rgl::rgl.clear()
topo <- raster::raster(system.file("extdata", "gebco1.tif", package = "anglr"))
poly <- subset(simpleworld, sovereignt == "India")
prj <- "+proj=laea +lon_0=80 +lat_0=45 +datum=WGS84"
poly_z <- anglr(poly, z = topo * 68, max_area = 0.02)
poly_z$v[c("x_", "y_")] <- rgdal::project(as.matrix(poly_z$v[c("x_", "y_")]), prj)
poly_z$meta$proj <- prj
plot(poly_z)
rgl::rglwidget()

A discrete transfer of z values can be done directly by copy_down.

data("polymesh", package = "silicate")
rgl::rgl.clear()
polymesh$val <- abs(polymesh$layer - 8) * 0.01
plot3d(copy_down(TRI(polymesh), z = polymesh$val))
rgl::rglwidget()