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

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

A conditioning plot, or coplot:

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"

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
## [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)

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

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

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.

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

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)

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 basic wire frame plot for the volcano data:

wireframe(z ~ x * y, data = vd, aspect = c(61 / 89, 0.3))

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

Interactive 3D Plots Using OpenGL

OpenGL is a standardized framework for high performance graphics.

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