## Other 3D Viewing Options

Other options for viewing a 3D scene include:

A stereogram

• presents each eye with an image from a slightly different viewing angle;
• the brain fuses the images into a 3D view in a process called stereopsis.

This approach is or was used You can create your own viewer with paper or cardboard tubes, for example.

A stereo image for the soil data surface:

``````tp <- trellis.par.get()
trellis.par.set(theme=col.whitebg())
ocol <- trellis.par.get("axis.line")\$col
oclip <- trellis.par.get("clip")\$panel
trellis.par.set(list(axis.line = list(col = "transparent"),
clip = list(panel = FALSE)))
print(wireframe(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9,
aspect=diff(range(soi.grid\$n))/diff(range(soi.grid\$e)),
xlab = "Easting (km)",
screen = list(z = 40, x = -70, y = 3)),
split = c(1,1,2,1), more = TRUE)

print(wireframe(soi.fit ~ soi.grid\$easting * soi.grid\$northing,
cuts = 9,
aspect=diff(range(soi.grid\$n))/diff(range(soi.grid\$e)),
xlab = "Easting (km)",
screen = list(z = 40, x = -70, y = 0)),
split = c(2,1,2,1))`````` ``````trellis.par.set(list(axis.line = list(col = ocol), clip = list(panel = oclip)))
trellis.par.set(tp)``````

## Some Calibration Experiments

### Simple Linear Stucture

Linear structure with normally distributed data:

``````n <- 1000
d1 <- data.frame(x = rnorm(n), y = rnorm(n))
d1 <- mutate(d1, z = x + y)``````
``splom(d1)`` ``````
xyplot(z ~ x | equal.count(y), data = d1)`````` ``xyplot(z ~ x | equal.count(y, 12), data = d1)`` ``xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d1)`` ``splom(d1, group = abs(x) < 0.1)`` ``````clear3d()
points3d(d1)``````

Linear structure with uniformly distributed data:

``````d2 <- data.frame(x = runif(n, -1, 1), y = runif(n, -1, 1))
d2 <- mutate(d2, z = x + y)``````
``splom(d2)`` ``xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d2)`` ``splom(d2, group = abs(x) < 0.1)`` ``````clear3d()
points3d(d2)``````

### Nonlinearity and Noise

Normal data with non-linearity and noise:

``d3 <- mutate(d1, z = cos(x) + cos(y) + 0.5 * (x + y) + rnorm(x, 0, .2))``
``splom(d3)`` ``xyplot(z ~ x | equal.count(y, overlap = -0.8), data = d3)`` ``splom(d3, group = abs(x) < 0.1)`` ``splom(d3, group = abs(x + 1) < 0.1)`` Uniform data with non-linearity and noise:

``d4 <- mutate(d2,  z = cos(x) + cos(y) + 0.5 * (x + y) + rnorm(x, 0, .1))``
``xyplot(z ~ x | equal.count(y, 12, overlap = -0.8), data = d4)`` ``splom(d4, group = abs(x) < 0.1)`` Full data in the background for context:

``````xyplot(z ~ x | equal.count(y, 12, overlap = -0.8), data = d4,
panel = function(x, y, ...) {
panel.xyplot(d4\$x, d4\$z, col = "grey")
panel.xyplot(x, y,...)
})`````` ``````clear3d()
points3d(d4)``````

### More Variables

Two-dimensional views, conditioning on two variables:

``````n <- 10000
d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
d <- mutate(d,
x4 = as.numeric(scale(x1 + cos(x2) + sin(x3) + rnorm(n, 0, 0.1))))``````
``splom(d)`` ``splom(~d, data = d, group = abs(x1) < 0.1)`` ``splom(~d, data = d, group = abs(x1) < 0.1 & abs(x2) < 0.1)`` ``````xyplot(x4 ~ x3 |
equal.count(x1, 4, overlap = -0.8) +
equal.count(x2, 4, overlap = -0.8),
data = d, asp = 1)`````` Three-dimensional views, conditioning on one variable:

``dd <- filter(d, abs(x1) < 0.1)[2:4]``
``````clear3d()
points3d(dd)``````

``````clear3d()

``````clear3d()
spheres3d(dd, radius = 0.05, col = ifelse(abs(dd\$x2) < 0.1, "red", "blue"))``````

``````dd1 <- filter(d, abs(x1 - 1) < 0.1)[2:4]
dd2 <- filter(d, abs(x1 + 1) < 0.1)[2:4]
clear3d()
spheres3d(dd1, radius = 0.05, color = rep("red", ncol(dd1)))
spheres3d(dd2, radius = 0.05, color = rep("blue", ncol(dd2)))``````

``````ddx <- mutate(d,
v = ifelse(abs(x1 - 1) < 0.1,
1,
ifelse(abs(x1 + 1) < 0.1, 2, NA)))
ddx <- filter(ddx, ! is.na(v))
cloud(x4 ~ x2 * x3, group = v, data = ddx, auto.key = TRUE)`````` ``cloud(x4 ~ x2 * x3 | v, data = ddx)`` ## Interactive Explorations

### Tk Sliders Using `tkrplot`

A simple animation function using the tkrplot` package:

``````library(tkrplot)

tkranim <- function(fun, v, lo, hi, res = 0.1, title, label) {
tt <- tktoplevel()
if (! missing(title))
tkwm.title(tt, title)
draw <- function() fun(v)
img <-tkrplot(tt, draw)
tv <- tclVar(init = v)
f<-function(...) {
vv <- as.numeric(tclvalue(tv))
if (vv != v) {
v <<- vv
tkrreplot(img)
}
}
s <- tkscale(tt, command=f, from=lo, to=hi, variable=tv,
showvalue=FALSE, resolution=res, orient="horiz")
if (! missing(label))
tkpack(img, tklabel(tt, text=label), s)
else tkpack(img,s)
}``````

Interactive conditioning for one of the calibration examples:

``````fd4 <- function(v)
print(splom(~ d4, group = abs(x - v) < 0.1, data = d4))

tkranim(fd4, 0, -1, 1, 0.05)``````

Interactive slicing for the soil resistivity surface:

``````f <- function(i) {
e <- eastseq[i]
idx <- which(sf\$easting == e)
r <- range(sf\$fit)
print(xyplot(fit ~ northing, data = sf[idx,],
ylim = r, type = "l", main = sprintf("Easting = %f", e)))
}

tkranim(f, 1, 1, 85, 1)``````

Interactive conditioning for the soil resistivity data:

``````binwidth <- 0.1 ## bad idea to use a global variable
g <- function(i) {
e <- eastseq[i]
idx <- which(abs(soil\$easting - e) <= binwidth)
ry <- range(soil\$resistivity)
rx <- range(soil\$northing)
print(xyplot(resistivity ~ northing, data = soil[idx,],
xlim = rx, ylim = ry, main = sprintf("Easting = %f", e)))
}

tkranim(g, 1, 1, 85, 1)``````

### Shiny Sliders

A version of the animation function using the `shiny` package:

``````shinyAnim <- function(fun, v, lo, hi, step = NULL, label = "Frame") {
require(shiny)
shinyApp(
ui = fluidPage(
sidebarLayout(
sidebarPanel(sliderInput("n", label, lo, hi, v, step,
animate = animationOptions(100))),
mainPanel(plotOutput("plot"))
)
),
server = function(input, output, session) {
session\$onSessionEnded(function() shiny::stopApp())
output\$plot <- renderPlot(fun(input\$n))
}
)
}``````
``shinyAnim(fd4, 0, -1, 1, 0.05)``
``shinyAnim(f, 1, 1, 85)``
``shinyAnim(g, 1, 1, 85, 1)``

### Animations as Movies

To show how a plot changes with a single parameter it is possible to pre-compute the images and assemble them into a movie.

This is easy to do in Rmarkdown, but requires the `ffmpeg` software to be installed.

Currently, `ffmpeg` is not available on the CLAS Linux systems.

Linked plots and brushing were pioneered at Bell Laboratories in the early 1980s.

The `iplots` package provides a set of basic plot that are all linked by default.

``````library(iplots)
with(soil, iplot(northing, resistivity))``````
``with(soil, iplot(easting,northing))``

The `rggobi` package provides another framework for linked plots and brusing within R.

``````library(rggobi)
ggobi(soil)``````

Several other tools with linking/brushing support and good integration with R are under development.

• One of the furthest along seems to be `rbokeh`.
• `plotly` also provides some support for linked brushing.

## Example: Ethanol Data

The `ethanol` data frame in the `lattice` package contains data from an experiment efficiency and emissions in small one-cylinder engines.

The data frame contains 88 observations on three variables:

• `NOx`: Concentration of nitrogen oxides (NO and NO2) in micrograms.

• `C` Compression ratio of the engine.

• `E` Equivalence ratio, a measure of the richness of the air and ethanol fuel mixture.

A scatterplot matrix:

``````library(lattice)
splom(ethanol)`````` Conditional plots of `NOx` against `E` for each level of `C`:

``xyplot(NOx ~ E | C, data = ethanol)`` ``````
ggplot(ethanol, aes(x = E, y = NOx)) + geom_point() + facet_wrap(~C)`````` A single plot with color to separate levels of `C`:

``xyplot(NOx ~ E, group = C, data = ethanol, auto.key = TRUE)`` ``````
ggplot(ethanol, aes(x = E, y = NOx, color = factor(C))) + geom_point()`````` Augmented with smooths:

``````xyplot(NOx ~ E, group = C, type = c("p", "smooth"), data = ethanol,
auto.key = TRUE)`````` ``````
ggplot(ethanol, aes(x = E, y = NOx, color = factor(C))) +
geom_point() + geom_smooth(se = FALSE)
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'`````` Conditional plots of `NOx` against `C` given levels of `E`:

``````with(ethanol, {
Equivalence.Ratio <- equal.count(E, number = 9, overlap = 0.25)
xyplot(NOx ~ C | Equivalence.Ratio,
panel = function(x, y) {
panel.xyplot(x, y)
panel.loess(x, y, span = 1)
},
aspect = 2.5,
layout = c(5, 2),
xlab = "Compression Ratio",
ylab = "NOx (micrograms/J)")
})`````` An interactive 3D view:

``````library(rgl)
clear3d()
with(ethanol,points3d(E, C / 5, NOx / 5,
size = 4, col = rep("black", nrow(ethanol))))``````

## Example: Earth Quake Locations and Magnitudes

The `quakes` data frame contains data on locations of seismic events of magnitude 4.0 or larger in a region near Fiji.

The time frame is from 1964 to perhaps 2000. More recent data is available from a number of sources on the web.

Some scatter plot matrices:

``````library(lattice)
splom(quakes)`````` ``splom(quakes, cex = 0.1)`` ``splom(quakes, alpha = 0.1)`` ``````splom(quakes, panel = panel.smoothScatter) Where the quakes are located:

``````library(maps)
map("world2", c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = "red", cex = 0.1))`````` ``````
map("world2", c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = heat.colors(nrow(quakes))[rank(depth)],
cex = 0.1))`````` ``````
map("world2",c("Fiji", "Tonga", "New Zealand"))
with(quakes, points(long, lat, col = heat.colors(nrow(quakes))[rank(depth)]))`````` Showing `depth` with grouping:

``````quakes2 <- mutate(quakes, depth_cut = cut_number(depth, 3))
xyplot(lat ~ long, aspect = "iso", groups = depth_cut,
data = quakes2, auto.key = list(columns = 3, title = "Depth"))`````` ``````
ggplot(quakes2, aes(x = long, y = lat, color = depth_cut) ) +
geom_point() + theme_bw() + coord_map()`````` ``````
ggplot(quakes2, aes(x = long, y = lat, colour = depth_cut, size = mag)) +
geom_point() + theme_bw() + coord_map()`````` Showing `depth` with a trellis display:

``xyplot(lat ~ long | depth_cut, aspect = "iso", data = quakes2)`` ``````
xyplot(lat ~ long | depth_cut, aspect = "iso", type=c("p","g"),
data = quakes2, auto.key = TRUE)`````` ``````
ggplot(quakes2, aes(x = long, y = lat)) +
geom_point(color = "blue") + facet_wrap(~depth_cut, nrow = 1) +
theme_bw() + coord_map()`````` Adding the full data for background context:

``````xyplot(lat ~ long | cut(depth, 3), aspect = "iso", type=c("p","g"),
panel = function(x, y, ...) {
panel.xyplot(quakes\$long, quakes\$lat, col = "grey")
panel.xyplot(x, y, ...)
},
data = quakes)`````` ``````
## quakes does not contain the depth_cut variable used in the facet
ggplot(quakes2, aes(x = long, y = lat)) +
geom_point(data = quakes, color = "gray") +
geom_point(color = "blue") + facet_wrap(~depth_cut, nrow = 1) +
theme_bw() + coord_map()`````` 3D views of `lat`, `long`, and `depth`:

``cloud(-depth ~ long * lat, data = quakes)`` ``````library(rgl)
bg3d(col = "black")
clear3d()
with(quakes, points3d(long, lat, -depth / 50, size = 2,
col = rep("white", nrow(quakes))))``````

``````clear3d()
with(quakes, points3d(long, lat, -depth / 50, size = 2,
col=heat.colors(nrow(quakes))[rank(mag)]))``````

## Visualizing 2D and 3D Density Estimates

### 2D Density Estimates

We have already seen some examples of plotting 2D density contours on a scatter plot using `geom_density2d`.

Here is a version that explicitly creates the density estimate:

``````data(geyser, package = "MASS")
geyser <- mutate(geyser, duration = lag(duration))
geyser <- filter(geyser, ! is.na(duration))
plot(waiting ~ duration, data = geyser)`````` ``````sc <- function(x) as.numeric(scale(x))
de <- with(geyser, MASS::kde2d(sc(duration), sc(waiting), n = 50,  h = 0.5))
contour(de)

with(geyser, points(sc(duration), sc(waiting)))`````` The density estimate can also be viewed with a perspective plot or an interactive 3D plot:

``persp(de)`` ``````clear3d()
bg3d(col = "white")
surface3d(unique(de\$x), unique(de\$y), de\$z, col = "red")
box3d()
axes3d()``````

### 3D Density Estimates

A 3D normal mixture:

``````library(rgl)
library(misc3d)

x <- matrix(rnorm(9000), ncol = 3)
x[1001 : 2000, 2] <- x[1001 : 2000, 2] + 4
x[2001 : 3000, 3] <- x[2001 : 3000, 3] + 4``````
``````clear3d()
bg3d(col = "white")
points3d(x = x[,1], y = x[,2], z = x[,3], size = 2, color = "black")``````