• A static graph is two-dimensional.

• Showing a third dimension is a challenge.

• For statistical data three dimensions is not particularly special.

• But we can take advantage of our ability to perceive three-dimensional structure.

## Soil Resistivity Data

Data from Cleveland’s Visualizing Data book contains measurements of soil resistivity of an agricultural field along a roughly rectangular grid.

Using functions from the `lattice` package we can view the locations at which measurements were taken and a perspective plot of the 3D point cloud:

``````soil <- read.table("soil.dat")
p <- cloud(resistivity ~ easting * northing, pch = ".", data = soil)
s <- xyplot(northing ~ easting, pch = ".", aspect = 2.44, data = soil)
print(s, split = c(1, 1, 2, 1), more = TRUE)
print(p, split = c(2, 1, 2, 1))`````` • The data is quite noisy but there is some structure.

• The viewing angle and distance for the point cloud can be adjusted.

## Conditioning Plots (Coplots)

One way to try to get a handle on higher dimensional data is to try to fix values of some variables and visualize the values of others in 2D.

This can be done with

• interactive tools;

• small multiples with trellis displays or faceting.

A conditioning plot, or coplot:

• Shows a collection of plots of two variables for different settings of one or more additional variables, the conditioning variables.

• For ordered conditioning variables the plots are arranged in a way that reflects the order.

• When a conditioning variable is numeric, or ordered categorical with many levels, the values of the conditioning variable are grouped into bins.

• A useful strategy is to allow the bins to overlap somewhat.

• Trellis/lattice plots support this with the concept of a shingle.
• This does not seem easy to do with `ggplot`.

A lattice plot of `resistivity` against `northing` conditioned on `easting`:

``xyplot(resistivity ~ northing | equal.count(easting), data = soil)`` `equal.count` produces a shingle:

``````str(equal.count(soil\$easting))
## Class 'shingle'  atomic [1:8641] 0.016 0.0252 0.0345 0.0437 0.0529 0.0621 0.0714 0.0806 0.0898 0.0991 ...
##   ..- attr(*, "levels")=List of 6
##   .. ..\$ : num [1:2] -0.00405 0.39615
##   .. ..\$ : num [1:2] 0.235 0.584
##   .. ..\$ : num [1:2] 0.396 0.795
##   .. ..\$ : num [1:2] 0.584 1.027
##   .. ..\$ : num [1:2] 0.795 1.273
##   .. ..\$ : num [1:2] 1.03 1.56
##   .. ..- attr(*, "class")= chr "shingleLevel"``````
• The default number of intervals in 6.
• The intervals used are highlighted in the title bars.
• The intervals are chosen to contain approximately equal numbers of observations.
• By default, the intervals overlap by 50%.
• The intervals are computed by `co.intervals`; the amount of overlap is controlled by the `overlap` argument.

Using smaller points to reduce over-plotting:

``xyplot(resistivity ~ northing | equal.count(easting), data = soil, cex = 0.2)`` A larger number of intervals can be specified as a second argument to `equal.bins`:

``````xyplot(resistivity ~ northing | equal.count(easting, 12), data = soil,
cex = 0.2)`````` ### Coplot Variations

To help see the signal we can add a smooth, by default computed using the `loess` smoother. This can be done with a `type` specification that includes the string `"smooth"`:

``````xyplot(resistivity ~ northing | equal.count(easting, 12), data = soil,
cex = 0.2, type = c("p", "smooth"))`````` The smooth is hard to see. Some options:

• Omit the data and only show the smooth.
• Show the data in a less intense color, such as gray.
• Use a contrasting color for the smooth curves.
• Show the data using alpha blending.

Color adjustments can be accomplished with the `col` and `col.line` arguments; the line width can be adjusted with `lwd`. For a thicker red line:

``````xyplot(resistivity ~ northing | equal.count(easting, 12), data = soil,
cex = 0.2, type = c("p", "smooth"), col.line = "red", lwd = 2)`````` For gray points and a red smooth:

``````xyplot(resistivity ~ northing | equal.count(easting, 12), data = soil,
cex = 0.2, type = c("p", "smooth"),
col.line = "red", col = "gray", lwd = 2)`````` The `xyplot` function supports an `alpha` argument, but it affects both points and lines.

One option for changing only the alpha level for points is to use a panel function:

``````xyplot(resistivity ~ northing | equal.count(easting, 12), data = soil,
cex = 0.2, panel = function(x, y, ...) {
panel.xyplot(x, y, ...)
panel.loess(x, y, col = "red", lwd = 2)})``````

Another option is to create a color for the points with the desired alpha level included in its encoding:

``````blue_w_alpha <- rgb(0, 0, 1, 0.1)
blue_w_alpha
##  "#0000FF1A"``````

This can then be used with the `col` argument:

``````xyplot(resistivity ~ northing | equal.count(easting, 12), data = soil,
cex = 0.2, type = c("p", "smooth"),
col.line = "red", col = blue_w_alpha, lwd = 2)`````` After conditioning on one of two explanatory variables it is useful to also consider conditioning on the other:

``````xyplot(resistivity ~ easting | equal.count(northing, 12), data = soil,
cex = 0.2, type = c("p", "smooth"),
col.line = "red", col = blue_w_alpha, lwd = 2)`````` ### Copots using `ggplot`

`ggplot` faceting is analogous to trellis/lattice conditioning.

• The main missing feature is the possibility of overlap among group.
• Every data point has to belong to exactly one facet.
• This also makes showing a selection of narrow windows more challenging.

The basic coplot:

``````p <- ggplot(soil, aes(northing, resistivity))
p + geom_point() + facet_wrap(~ cut_number(easting, 6))`````` • `cut_number` is used to form non-overlapping groups with approximately equal numbers of observations.
• `cut_interval` can be used to divide the range in equal intervals.
• The facet ordering is top to bottom; lattice orders bottom to top.

Using 12 facets and smaller points:

``p + geom_point(size = 0.5) + facet_wrap(~ cut_number(easting, 12))`` The default smoothing method is currently “loess” for less than 1000 points:

``````p + geom_point(size = 0.5) + geom_smooth(color = "red") +
facet_wrap(~ cut_number(easting, 12))
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'`````` For larger data sets `method = "gam"` would be used:

``````p + geom_point(size = 0.5) + facet_wrap(~ cut_number(easting, 12)) +
geom_smooth(method = "gam", color = "red", formula = y ~ s(x))`````` ## Contour Plots for Surfaces

The `loess` local polynomial smoother can be used to estimate a smooth signal surface as a function of the two location variables.

The estimated surface level can be computed on a grid of points using the `predict` method of the fit.

``````eastseq <- seq(.15, 1.410, by = .015)
northseq <- seq(.150, 3.645, by = .015)
soi.grid <- expand.grid(easting = eastseq, northing = northseq)
##   easting northing
## 1   0.150     0.15
## 2   0.165     0.15
## 3   0.180     0.15
## 4   0.195     0.15
## 5   0.210     0.15
## 6   0.225     0.15
m <- loess(resistivity ~ easting * northing, span = 0.25,
degree = 2, data = soil)
soi.fit <- predict(m, soi.grid)``````

Contour plots are one approach to visualizing a surface based on values on a grid.

• Contour plots compute contours, or level curves, as polygons at a set of levels.

• Contour plots draw the level curves, often with a level annotation.

• Contour plots can also have their polygons filled in with colors representing the levels.

A `lattice` contour plot for the fitted surface:

``````asp <- diff(range(soi.grid\$n)) / diff(range(soi.grid\$e))
contourplot(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9, aspect = asp,
xlab = "Easting (km)",
ylab = "Northing (km)")`````` A basic contour plot in `ggplot` using `geom_contour`:

``````p <- ggplot(mutate(soi.grid, fit = as.numeric(soi.fit)),
aes(x = easting, y = northing, z = fit)) +
coord_fixed()
p + geom_contour()`````` • Neither `lattice` nor `ggplot` seem to make it easy to fill in the contours.

• The base function `filled.contour` is available for this:

``````cm.rev <- function(...) rev(cm.colors(...))
filled.contour(eastseq, northseq, soi.fit, asp = 1, color.palette = cm.rev)`````` ## Level Plots

A level plot colors a grid spanned by two variables by the color of a third variable.

• Level plots are also called image plots

• The term heat map is also used, in particular with a specific color scheme. But heat map often means a more complex visualization with an image plot at its core.

• Superimposing contours on a level plot is often helpful.

The `lattice` package provides the function `levelplot`.

A basic level plot:

``````levelplot(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9, aspect = asp,
xlab = "Easting (km)",
ylab = "Northing (km)")`````` A level plot with superimposed contours:

``````lv <- levelplot(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9, aspect = asp, contour = TRUE,
xlab = "Easting (km)",
ylab = "Northing (km)")
lv`````` `ggplot` provides `geom_tile` that can be used for a level plot:

``````p + geom_tile(aes(fill = fit)) + geom_contour() + • Level plots do not require computing contours, but are not not as smooth as filled contour plots.

• Visually, image plots and filled contour plots are very similar for fine grids, but image plots are less smooth for coarse ones.

• Lack of smoothness is less of an issue when the data values themselves are noisy.

The grid for the `volcano` data set is coarser and illustrates the lack of smoothness.

``````vd <- expand.grid(x = seq_len(nrow(volcano)), y = seq_len(ncol(volcano)))
vd\$z <- as.numeric(volcano)
levelplot(z ~ x * y, data = vd, cuts = 10)`````` A `filled.contour` plot looks like this:

``filled.contour(volcano, nlevels = 10, color.palette = cm.rev)`` • A coarse grid can be interpolated to a finer grid.

• Irregularly spaced data can also be interpolated to a grid.

• The `interp` function in the `akima` package is useful for this kind of interpolation.

## Wire Frame and Perspective Plots

A surface can also be visualized using a wire frame plot showing a 3D view of the surface from a particular viewpoint.

• A simple wire frame plot is often sufficient.

• Lighting and shading can be used to enhance the 3D effect.

A basic wire frame plot for the volcano data:

``wireframe(z ~ x * y, data = vd, aspect = c(61 / 89, 0.3))`` • Wire frame is a bit of a misnomer since surface panels in front occlude lines behind them.

• For a fine grid, as in the soil surface, the lines are too dense.

• The use of shading for the surfaces can help.

``wireframe(z ~ x * y, data = vd, aspect = c(61 / 89, 0.3), shade = TRUE)`` A wire frame plot with shading for the soil resistivity data:

``````wf <- wireframe(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
aspect = asp, shade = TRUE,
xlab = "Easting (km)",
ylab = "Northing (km)")
wf`````` Both ways of looking at a surface are useful:

``````print(lv, split = c(1, 1, 2, 1), more = TRUE)
print(wf, split = c(2, 1, 2, 1))`````` • The level plot/contour representation is useful for recognizing locations of key features.
• The wire frame view helps build a mental model of the 3D structure.
• Being able to interactively adjust the viewing position for a wire frame model greatly enhances our ability to understand the 3D structure.

## Interactive 3D Plots Using OpenGL

OpenGL is a standardized framework for high performance graphics.

• The `rgl` package provides an R interface to some of OpenGL’s capabilities.

• WebGL is a JavaScript framework for using OpenGL within a browser window.

• Most desktop browsers support WebGL; some mobile browsers do as well.

• In some cases support may be available but not enabled by default. You may be able to get help at https://get.webgl.org/.

• `knitr` and `rgl` provide support for embedding OpenGL images in web pages.

• It is also possible to embed OpenGL images in PDF files, but not all PDF viewers support this.

This code run in R will open a new window containing an interactive 3D scene (but this currently fails on FastX):

``````library(rgl)
bg3d(color = "white")
clear3d()
par3d(mouseMode="trackball")
surface3d(eastseq, northseq,
soi.fit / 100, color = rep("red", length(soi.fit)))``````

To embed an image in an HTML, document first set the `webgl` hook with a code chunk like this:

``````knitr::knit_hooks\$set(webgl = rgl::hook_webgl)
options(rgl.useNULL=TRUE)``````

Then a chunk with the option `webgl = TRUE` can produce an embedded OpenGL image:

``````library(rgl)
bg3d(color = "white")
clear3d()
par3d(mouseMode="trackball")
surface3d(eastseq, northseq,
soi.fit / 100, color = rep("red", length(soi.fit)))``````