Chapter 3 3D and mesh forms of spatial data

3.1 Questions

  • What are meshes and topology?
  • what is the relationship of meshes to geospatial raster/vector data?

3.2 What is a mesh?

Key ideas!

  • Topology: the shape of things and their relationships.
  • Geometry: the where of things.
  • Indexing: the link between topology and geometry.

Traditional spatial data tends to confuse topology and geometry.

These concepts motivate my interest in these ideas.

  • Lossless reprojection
  • Topology fixes
  • Tracks and point clouds
  • Visualization

Topology vs. geometry

This line has 1-dimensional topology depicted in 3-dimensional geometry and the triangle has 2-dimensional topology depicted in 3-dimensional geometry.

##        x   y   z
## [1,] 0.0 0.0 0.0
## [2,] 0.5 0.5 0.0
## [3,] 1.0 1.0 0.8
##     [,1] [,2]
## .v0    1    2
## .v1    2    3
##     [,1]
## .v0    1
## .v1    2
## .v2    3

3.3 Geospatial data

  • raster
  • vector

These are traditionally kept separate, but in computer graphics the distinction starts to disappear.

What is a raster?

A layer of neighbouring rectangles?

Or a continuous fields between points?

(Lots of ways to infer the field, including this very poor one).

What is a polygon?

A series of grouped paths?

What’s in the middle?

## Warning in st_centroid.sf(tri): st_centroid assumes attributes are constant
## over geometries of x

The fill we see in traditional 2D graphics is a trick!!.

Search:

it’s not what you draw it’s what you not draw ~Paul Murrell

Technically the trick comes in two types, either the even-odd or winding rule, and this trick is not part of this workshop. The graphics engine uses this rule to draw a pixel if it has been encircled an even or odd number of times, or using a rule about in which direction it was encircled. It happens deep in the graphics.

Where it does matter is for the concept of orientation, and 3D graphics do care about the direction that triangles are wound (consider that reversing the direction is like flipping the triangle in place in terms of how some algorithms behave …).

What’s the fill?

In 3D, and to fill our polygons properly as data - we need primitives.

3.4 Primitives

Terminology alert! (This is my working definition of primitive for this topic, not everyone agrees.)

  • Point - a single coordinate is a 0-dimensional primitive (vertex, coordinate)
  • Line - a line segment between two coordinates is a 1-dimensional primitive (edge, segment)
  • Triangle - a triangle joining three coordinates is 2-dimensional primitive
  • Quad - a four-sided shape (not as core-primitive as a triangle, but commonly used)

Topology is not geometry

(This is topological dimension. Every one of these types of shape can be depicted within a geometric space that is equal to or higher than the topological dimension.)

We will have a matrix of vertices and a matrix of primitive indices. Quads and triangles are generally called faces, line segments are alternatively called edges. All are primitives in computer graphics, but we’ll also see the term finite element used.

Topology can be 3D (tetrahedron) - imagine volumetric fill versus planar faces bounding a volume. Geometry can be 4D, X, Y, Z, T - or any dimensionality we need.

To fill our polygon we need triangles.

## Warning in st_cast.sf(sfdct::ct_triangulate(minimal_mesh, a = 0.01, D =
## TRUE)): repeating attributes for all sub-geometries for which they may not
## be constant

Note that if we turn off the border, we don’t notice the difference.

No tricky winding or even-odd rule to worry about, but we have lost our boundary around each distinct shape - we could find them by finding edges within a shape that are not shared by two triangles …

## Warning in plot.sf(tri, border = NA, col = rainbow(10)): col is not of
## length 1 or nrow(x): colors will be recycled; use pal to specify a color
## palette

Raster and vector are not a strong distinction when it comes to meshes

A raster is a very simple version of a mesh. When we store a raster we need the equivalent of

  • number of columns and rows
  • the extent in the coordinate system used (xmin, xmax, ymin, ymax)
  • the coordinate system
  • the cell values! (ncols * nrows of them)

In in computer graphics we store

  • the corner coordinates ((ncols + 1) * (nrows + 1) of them)
  • an index, 4 indices for every quad specify the coordinates
  • groupings (or at least material-properties), what quads belong to which objects
  • the coordinate system (hopefully)

It’s the same for a triangular mesh.

  • the corner coordinates
  • an index, 3 indices for every triangle
  • groupings, what triangles belong to which objects
  • the coordinate system (hopefully)

And lines are

  • the end point coordinates of each line segment (or edge)
  • an index, 2 indices for every segment
  • groupings, what line segments belong to which objects (features)
  • the coordinate system (hopefully)

3.5 A raster is a mesh (implicitly)

The simplest kind of mesh is a basic raster. Consider the matrix from above.

On its own this matrix has absolutely nothing to do with spatial data, it is literally a collection of 9 numeric values in a given order, and by the magic of programming we’ve nominated a shape of 3x3. We can’t help but think about this shape spatially however, but there’s a problem. Does each element occupy space or should we consider them to be infinitesimal locations?

R provides either interpretation (to simplify this story we nominate locations for the rows and columns explicitly).

When considered as an image, each matrix element occupies a certain space in width and height, but when considered as a point set the numbers simply float at the given locations. Which is correct? (Spoiler: Both are correct, it simply depends what we are doing.)

The raster package defaults to the image interpretation and helpfully assumes the values are nominally at the centre points as shown above. We have to nominate the extent or we end up in 0,1 range, we also have to invert the order of the values because raster counts from the top of the page and R’s matrix uses column-major order.

## class      : RasterLayer 
## dimensions : 3, 3, 9  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 0, 3, 0, 3  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : memory
## names      : layer 
## values     : 0, 0.5  (min, max)

R’s image and rasters in general are so efficient because they only store this minimal amount of information: the actual data values, and the extent and dimensions of the space they occur in. If we had to store the centre coordinate of every cell, or worse the corner coordinates then the data storage goes up dramatically. Every software that deals well with these kinds of data has to treat these coordinates as implicit. We can easily expand the centre coordinates.

##     x   y layer
## 1 0.5 2.5 0.250
## 2 1.5 2.5 0.125
## 3 2.5 2.5 0.000
## 4 0.5 1.5 0.375
## 5 1.5 1.5 0.500
## 6 2.5 1.5 0.375
##     x   y layer
## 4 0.5 1.5 0.375
## 5 1.5 1.5 0.500
## 6 2.5 1.5 0.375
## 7 0.5 0.5 0.000
## 8 1.5 0.5 0.125
## 9 2.5 0.5 0.250

but to expand the corners we have to jump through some hoops and even then we get every instance of the corners, not only for each cell but to explicitly close the cell as a polygon.

## class       : SpatialPointsDataFrame 
## features    : 45 
## extent      : 0, 3, 0, 3  (xmin, xmax, ymin, ymax)
## crs         : NA 
## variables   : 4
## names       : layer, Lines.NR, Lines.ID, Line.NR 
## min values  :     0,        1,        1,       1 
## max values  :   0.5,        9,        9,       1

The polygon version of the raster, converted to points has 45 coordinates in it, although there are only 16 unique coordinates at the corners, which is where mesh3d comes in.

3.6 Summary

  • Topology and geometry are independently dimensional
  • Meshes include a matrix of vertices and a matrix of indices
  • In mesh-terms rasters and polygons aren’t so different