In r-gris/gris: A Flexible Geometry and Topology Model for Spatial Data in R

knitr::opts_chunk$set(echo = TRUE)

junk <- capture.output({
library(raster)
library(maptools)
library(rgl)
library(rglwidget)  ## allows embedding OpenGL vis in RMarkdown HTML
library(gris) 
  library(dismo)
  library(geosphere)

})
rgl.material(specular = "grey10")

The gris package aims to provide a flexible framework for building topological spatial data structures in R. The model is flexible as it allows for straightforward translation between many other such models. Where possible, everything is stored as relational tables which also allows for straightforward transport from and to actual databases.

Motivation

Features of the gris approach are:

• ability to store traditional GIS data including multi-part polygons, lines and points
• ability to convert part-based (or branch-based) objects to segment and triangle primitive forms required by graph and graphics tools
• no limitations on the attributes that can be stored on vertices, or on the branches / primitives that compose them
• stable links from objects to their component branches/primitives and again to vertices
• ability to store complex surfaces and lines in 3D, 4D and higher geometries, not just "2.5D" TINs
• Constrained Delaunay triangulation

GIS data structures are not well suited for generalization, and visualizations and models in 3D require pretty forceful and ad hoc approaches. By the time a GIS polygon layer is converted to arrays of triangulated vertices, there's often no natural link back to the original data. Here we describe several examples:

• convert polygons to 3D and visualize as a globe
• copy terrain values to polygon objects using internal vertices

First we must dig into a set of polygons, and make our first 3D plot in R.

GIS polygons R

The R package maptools contains an in-built data set called wrld_simpl, which is a basic (and out of date) set of polygons describing the land masses of the world by country. This code loads the data set and plots it with a basic grey-scale scheme for individual countries.

library(maptools)
data(wrld_simpl)
print(wrld_simpl)
plot(wrld_simpl, col = grey(sample(seq(0, 1, length = nrow(wrld_simpl)))))

We also include a print statement to get a description of the data set, this is a SpatialPolgyonsDataFrame which is basically a table of attributes with one row for each country, linked to a recursive data structure holding sets of arrays of coordinates for each individual piece of these complex polygons.

These structures are quite complicated, involving nested lists of matrices with X-Y coordinates. First I will extract all of the coordinates from this structure in order to make our first 3D plot. I can use class coercion from polygons, to lines, then to points as the most straightforward way of obtaining every XY coordinate by dropping the recursive hierarchy structure to get at every single vertex in one matrix. (There are other methods to obtain all coordinates while retaining information about the country objects and their component "pieces", but I'm ignoring that for now.)

allcoords <- coordinates(as(as(wrld_simpl, "SpatialLines"), "SpatialPoints"))
dim(allcoords)
head(allcoords)  ## print top few rows

We need to put these "X/Y" coordinates in 3D so I simply add another column filled with zeroes.

allcoords <- cbind(allcoords, 0)
head(allcoords)

OpenGL in R

In R we have access to 3D visualizations in OpenGL via the rgl package, but the model for data representation is very different to that of the polygons so I first plot the vertices of the wrld_simpl layer as points only.

library(rgl)
library(rglwidget)  ## allows embedding OpenGL vis in RMarkdown HTML
plot3d(allcoords, xlab = "", ylab = "") 
subid <- currentSubscene3d()
rglwidget(elementId="plot3d_wrld")
rgl.close()

Plotting in the plane is one thing, but more striking is to convert the vertices from planar longitude-latitude to Cartesizan XYZ. Define an R function to take "longitude-latitude-height" and return spherical coordinates (we can leave WGS84 for another day).

## deploy our custom function on the longitude-latitude values
xyzcoords <- geocentric(allcoords)

Now we can visualize these XYZ coordinates in a more natural setting, and even add a blue sphere for visual effect.

plot3d(xyzcoords, xlab = "", ylab = "")
spheres3d(0, 0, 0, radius = 6370000, col = "lightblue")
subid <- currentSubscene3d()
rglwidget(elementId="plot3d_wrldxyz")
rgl.close()

This is still not very exciting, since our plot knows nothing about the connectivity between vertices. Before continuing we have to revisit how our polygons are organized.

Organization of polygons

The in-development R package gris provides a way to represent spatial objects as a set of relational tables. In short a gris object has tables "o" (objects), "b" (for branches), "bXv" (links between branches and vertices) and "v" the vertices.

It can also have tables 'oXt' (links between objects and triangles) and 'tXv' (links between triangles and vertices) - this is not yet generalized to deal with line segments or multi-points rather than triangles so we ignore it for now.

These table structures are described in more detail in the "gris-tables" vignette.

library(maptools)
data(wrld_simpl)
gobject <- gris(wrld_simpl)

More 3D already!

Why go to all this effort just for a few polygons? The structure of the gris objects gives us much more flexibility, so I can for example store the XYZ Cartesian coordinates right on the same data set. I don't need to recursively visit nested objects, it's just a straightforward calculation and update - although we're only making a simple point, this could be generalized a lot more for user code.

gobject$v$zlonlat <- 0
do_xyz <- function(table) {
  xyz <- llh2xyz(dplyr::select(table, x, y, zlonlat))
  table$X <- xyz[,1]
  table$Y <- xyz[,2]
  table$Z <- xyz[,3]
  table
}

gobject$v <- do_xyz(gobject$v)

gobject$v

I now have XYZ coordinates for my data set, and so for example I will extract out a few nearby countries and plot them. First plot in the traditional way, using the default "evenodd" rule used by R's polypath() function.

localarea <- gobject[gobject$o$NAME %in% c("Australia", "New Zealand"), ]
## plot in traditional 2d
plot(localarea, col = c("dodgerblue", "firebrick"))

The plot is a bit crazy since parts of NZ that are over the 180 meridian skews everything; while we could fix that easily by modifiying the vertex values for longitude it is naturally more sensible in 3D.

plot3d(localarea$v$X, localarea$v$Y, localarea$v$Z, xlab = "", ylab = "")
subid <- currentSubscene3d()
rglwidget(elementId="plot3d_wrldoznz")
rgl.close()

The R package RTriangle wraps Richard Shewchuk's Triangle library, allowing constrained Delaunay triangulations. To run this we need to make a Planar Straight Line Graph from the polygons, but this is fairly straightforward by tracing through paired vertices in the data set. The key parts of the PSLG are the vertices P and the segment indexes S defining paired vertices for each line segment. This is a "structural" index where the index values are bound to the actual size and shape of the vertices, as opposed to a more general but perhaps less efficient relational index.

Maybe put this detail in another doc??

pslgraph <- mkpslg(gobject)
dim(pslgraph$P)
range(pslgraph$S)
head(pslgraph$P)
head(pslgraph$S)

The PSLG is what we need for the triangulation, but the triangulat.gris method handles this work for us.

tobject <- triangulate(gobject)

The triangulation vertices (long-lat) can be converted to XYZ, and plotted.

plot3d(tobject, globe = TRUE)
subid <- currentSubscene3d()
rglwidget(elementId="plot3d_wrld_polyxyz")
rgl.close()

These are very ugly polygons since there's no internal vertices to carry the curvature of this sphere. This is the same problem we'd face if we tried to drape these polygons over topography, as some point we need internal structure.

Luckily Triangle can set a minimum triangle size. We set a constant minimum area, which means no individual triangle can be larger in area than so many "square degrees". This gives a lot more internal structure so the polygons are more elegantly draped around the surface of the sphere. (There's not really enough internal structure added with this minimum area, but I've kept it simpler to make the size of this document more manageable).

tobject2 <- triangulate(gobject, a = 5)  ## a (area) is in degrees, same as our vertices
plot3d(tobject2)
subid <- currentSubscene3d()
rglwidget(elementId="plot3d_wrld_surfacexyz")
rgl.close()

Github, sfr

There is ongoing work to upgrade R's Spatial support, but this is still wedded to the basic flat and dumb polygonal model used by open source projects.

Image textures

## marmap has some topography
library(marmap)
data("hawaii", package = "marmap")
r <- marmap::as.raster(hawaii)
## but is otherwise undesirable so we unload it
unloadNamespace("marmap")
library(dismo)
gm <- gmap(r, type = "satellite", scale = 2)
#gm <- raster()
## this is our texture
paletteim2RGB <- function(x) {
  ct <-  col2rgb(x@legend@colortable[values(x) + 1])

  setValues(brick(x, x, x), t(ct))
}
library(rgdal)
texture <- writeGDAL(as(paletteim2RGB(gm), "SpatialGridDataFrame"), "texture.png", drivername = "PNG", type = "Byte")

library(gris)
library(rgl)
library(rglwidget)
b <- quadmeshFromRaster(r, r)
tcoords <- xyFromCell(setExtent(gm, c(0, 1, 0, 1)), cellFromXY(gm, project(t(b$vb[1:2, ]), projection(gm))))
aspect3d(1, 1, 0.0001)
#shade3d(b, col = "white", texture= texture, texcoords = tcoords[b$ib, ])

quads3d(t(b$vb[1:3, b$ib]), col = "white", texture= texture, texcoords = tcoords[b$ib, ])
subid <- currentSubscene3d()
rglwidget(elementId="plot3d_texture")
rgl.close()

Polygon pathologies

Orthodromes, loxodromes.

Densifying vertices

Projection boundaries

D3 approach to curvature

Real world topography

A simplistic layer of country polygons wrapped around a sphere is not very practical. Here we show the techniques applied above to build a relief map.

Extract a polygon map for Tasmania.

oz <- getData(name = "GADM", country = "AUS", level = 1)
tas <- disaggregate(subset(oz, NAME_1 == "Tasmania"))[-c(1, 2), ] ## drop MI

gtas<- gris(tas)

topo <- getData("SRTM", lon = mean(gtas$v$x), lat = mean(gtas$v$y))

f <- "ftp://xftp.jrc.it/pub/srtmV4/tiff/srtm_66_21.zip"
download.file(f, basename(f), mode = "wb")
unzip(basename(f))
topo <- raster("srtm_66_21.tif")

prj <- "+proj=laea +lon_0=147 +lat_0=-42 +ellps=WGS84"
xy <- rgdal::project(gtas$v %>% select(x, y) %>% as.matrix, prj)
gtas$v$x <- xy[,1]; gtas$v$y <- xy[,2]
tri <- RTriangle::triangulate(mkpslg(gtas), a = 5e6)

z <- extract(topo, rgdal::project(tri$P, prj, inv = TRUE))
triangles3d(cbind(tri$P, z)[t(tri$T), ])
aspect3d(1, 1, 2)

## getData('ISO3') 
nz <- getData("GADM", country = "NZL", level = 0)

## explode into separate pieces
nz <- disaggregate(nz)[order(areaPolygon(nz), decreasing = TRUE)[1:6], ]
## grisify and modify all vertices that are too far west in -180,180
gnz <- gris(nz[2,])
gnz$v <- gnz$v %>% mutate(x = ifelse(x < 0, x + 360, x))
#topo <- getData("SRTM", lon = mean(gnz$v$x), lat = mean(gnz$v$y))

Colophon? I use the programming environment R for the data manipulation and the creation of this document via several extensions (packages) to base R.

r-gris/gris documentation built on May 14, 2019, 12:57 a.m.