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.

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.

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`

.

- Trellis/lattice plots support this with the concept of a

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)
```

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
## [1] "#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)
```

`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))
```

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)
head(soi.grid)
## 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)
```

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() +
scale_fill_gradientn(colors = rev(cm.colors(100)))
```

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.

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.

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)))
```