\documentclass[article,nojss]{jss}
\DeclareGraphicsExtensions{.pdf,.eps,.png}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Add-on packages and fonts


\newcommand{\noun}[1]{\textsc{#1}}
%% Bold symbol macro for standard LaTeX users
\providecommand{\boldsymbol}[1]{\mbox{\boldmath $#1$}}

%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\newcommand{\p}{\textbf{\textsf{plot3D }}}
\newcommand{\R}{\proglang{R}}
\title{
  \p: Tools for plotting 3-D and 2-D data.
}
\Plaintitle{plot3D}

\Keywords{
  plot,  persp, image, 2-D, 3-D, scatter plots, surface plots, slice plots, 
  oceanographic data, \proglang{R}
}

\Plainkeywords{
  plot,  persp, image, 2-D, 3-D, scatter plots, surface plots, slice plots, 
  oceanographic data, R
}


\author{Karline Soetaert\\
NIOZ-Yerseke\\
The Netherlands
}

\Plainauthor{Karline Soetaert}

\Abstract{
  \R{ }package \p \citep{plot3D} contains functions for plotting multi-dimensional
  data. Many functions are derived from the \code{persp} function, other
  functions start from the \code{image} or \code{contour} function.

  Two related packages are: 
  \begin{itemize}
   \item   \pkg{plot3Drgl} \citep{plot3Drgl}, that plots multidimensional data
    using openGL graphics (and using package \pkg{rgl} \citep{rgl}).
   \item  \pkg{OceanView} \citep{OceanView} that contains functions for visualing
    oceanographic data.  
  \end{itemize}
  
  A graphical gallery using one of  \pkg{plot3D}, \pkg{plot3Drgl} or \pkg{OceanView}
  is in \url{http://www.rforscience.com/rpackages/visualisation/oceanview/}
  and \url{http://www.rforscience.com/rpackages/visualisation/plot3d/}
}

%% The address of (at least) one author should be given
%% in the following format:
\Address{
  Karline Soetaert\\
  Royal Netherlands Institute of Sea Research (NIOZ)\\
  4401 NT Yerseke, Netherlands\\
  E-mail: \email{karline.soetaert@nioz.nl}\\
  URL: \url{http://http://www.nioz.nl/}\\
}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% R/Sweave specific LaTeX commands.
%% need no \usepackage{Sweave}
%\VignetteIndexEntry{plot3D: functions for plotting 3-D and 2-D data}
%\VignetteKeywords{plot, persp, image, 2-D, 3-D, scatter, surface plots, slice plots}
%\VignettePackage{plot3D}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Begin of the document
\begin{document}
\SweaveOpts{engine=R,eps=FALSE,resolution=80}
\SweaveOpts{keep.source=TRUE}

<<preliminaries,echo=FALSE,results=hide>>=
library(plot3D)
options(prompt = " ")
options(continue = "  ")
options(width=75)
@

\maketitle

\section{Introduction}

\R{ } package \p provides functions for plotting 2-D and 3-D data, and 
that are either extensions of R's \code{persp} function
or of R's \code{image} and \code{contour} function.

The main extensions to these functions are:
\begin{itemize}
 \item In addition to the x, y (and z) values, an additional data dimension 
 can be represented by a color variable (argument \code{colvar}). 
 \item A color key (argument \code{colkey}) can be written next to the figure.  
   It is possible to log-transform the color key, rescale it, adjust its position, ...   
 \item The resolution of a figure can be increased (argument \code{resfac}). 
 \item Either the \code{facets} can be colored, just the border, or both.    
\end{itemize}

Package \p contains:
\begin{itemize}
 \item Functions that are based on the \code{persp} function, for visualising 3-D data:
   \begin{itemize}
     \item \code{persp3D}: an extended version of the \code{persp} function.
     \item \code{ribbon3D}: perspective plots as ribbons.
     \item \code{hist3D}: 3-D histograms.
     \item \code{scatter3D}, \code{points3D}, \code{lines3D}, \code{text3D}: 
       scatter plots in 3-D, points, lines, labels.
     \item \code{surf3D}: 3-D shapes (or surfaces).
     \item \code{slice3D}, \code{slicecont3D}, \code{isosurf3D}, \code{voxel3D}: 
     slices, isosurfaces and voxels from a full 3-D data set.
     \item \code{arrows3D}: arrows in 3D.
     \item \code{contour3D}, \code{image3D}: contours and images in 3D.
     \item \code{segments3D}, \code{polygon3D}, \code{rect3D}, \code{border3D}, \code{box3D}: 
       line segments, polygons, rectangles, boxes in 3D.
   \end{itemize}
 \item Functions defined on the \code{image} or \code{contour} function:
   \begin{itemize}
     \item \code{image2D}, \code{contour2D}, for an extended version of these functions to visualise 2-D (or 3-D) data.
     \item \code{ImageOcean}, for an image of the ocean's bathymetry.
   \end{itemize}
 \item Other functions
  \item \code{scatter2D}: colored points, lines, ... in 2-D. 
  \item \code{text2D}, \code{arrows2D}, \code{segments2D}, \code{rect2D}, 
    \code{polygon2D} for other 2D functions, comparable to R's base graphics but
    that have a color key.
 \item Colors and colorkeys:
   \begin{itemize}
     \item \code{colkey}: color legends.
     \item \code{jet.col}, \code{jet2.col}, \code{gg.col}, \code{ramp.col}: suitable color palettes.
   \end{itemize}
 \item Utility functions:
   \begin{itemize}
     \item \code{mesh}: generating rectangular (2D) or (3D) meshes.
     \item \code{plotdev}: plotting on the current device.
   \end{itemize}
 \item Data sets:
   \begin{itemize}
     \item \code{Oxsat}: a (rather large) 3-D data set with the 
      ocean's oxygen saturation values. 
     \item \code{Hypsometry}: a 2-D data set with the worlds elevation and 
       the ocean's depth.
   \end{itemize}
\end{itemize}

This vignette contains some examples; more can be found in the package's help files.
To run all examples:
\begin{verbatim}
 example(persp3D)  
 example(surf3D)  
 example(slice3D)
 example(scatter3D)
 example(segments3D)
 example(image2D)
 example(image3D)
 example(contour3D)
 example(colkey) 
 example(jet.col) 
 example(perspbox)
 example(mesh)
 example(trans3D)
 example(plot.plist)
 example(ImageOcean)
 example(Oxsat)
 example(legendplot)
\end{verbatim}

\section{Functions image2D and persp3D}
\code{image2D} and \code{persp3D} are extensions of R's \code{image} and \code{persp} 
functions.
The arguments of \code{persp3D} are (see the help file for what they mean):
<<>>=
args(persp3D)
@

Many examples of the use of \code{image2D} and \code{persp3D} are in vignette \code{volcano}. 

The \code{Hypsometry} data set is depicted first as an \code{image}, 
with 0 m contour lines added. 
Slight shading gives the plot a perspective view. 
The zoomed region (used in next figure) is then added.
<<label=hyps,include=FALSE>>=
image2D(Hypsometry, xlab = "longitude", ylab = "latitude", 
      contour = list(levels = 0, col = "black", lwd = 2),
      shade = 0.1, main = "Hypsometry data set", clab = "m")
rect(-50, 10, -20, 40, lwd = 3)
@

\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=hyps,fig=TRUE,echo=FALSE, width = 10, height = 6, pdf = FALSE, png = TRUE>>=
<<hyps>>
@
\end{center}
\caption{Hypsometry data set}
\label{fig:hyps}
\end{figure}


<<>>=
ii <- which(Hypsometry$x > -50 & Hypsometry$x < -20)
jj <- which(Hypsometry$y >  10 & Hypsometry$y <  40)
zlim <- c(-10000, 0)
@
The perspective figure is made with black side-panels (\code{bty}). 
Grey contour lines are added on the bottom panel (\code{"zmin"}) and on 
the persp plot itself (\code{"z"}).
The resolution is increased (\code{resfac}) to make smoother images.
A color key (\code{colkey}) is added on the first margin (\code{side})
<<label=ocean2,include=FALSE>>=
par(mfrow = c(1, 1))

# Actual bathymetry, 4 times increased resolution, with contours
persp3D(z = Hypsometry$z[ii,jj], xlab = "longitude", bty = "bl2",
        ylab = "latitude", zlab = "depth", clab = "depth, m", 
        expand = 0.5, d = 2, phi = 20, theta = 30, resfac = 2,  
        contour = list(col = "grey", side = c("zmin", "z")),
         zlim = zlim, colkey = list(side = 1, length = 0.5))
@

\setkeys{Gin}{width=0.6\textwidth}
\begin{figure}
\begin{center}
<<label=ocean2,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<ocean2>>
@
\end{center}
\caption{Bathymetry of a part of the ocean}
\label{fig:ocean2}
\end{figure}

\newpage
\section{slices and isosurfaces}
Function \code{slice3D} draws slices from volumetric (3D) data, 
function \code{isosurf3D} creates and plots isosurfaces. 
It makes use of a function from package \code{misc3d} \citep{misc3d}.

<<>>=
args(slice3D)
@
Function \code{mesh} is used to generate a full rectangular 3-D mesh. 
This is used to generate the volumetric data (\code{p}) that defines the coloration.
The data are visualised by one slice in x (\code{xs}) and 3 slices in y direction
(\code{ys}).
Function \code{isosurf3D} plots the data for p-values that are equal to \code{0}.

<<label=slice,include=FALSE>>=
par(mfrow = c(1, 2))
x <- y <- z <- seq(-4, 4, by = 0.2)
M <- mesh(x, y, z)

R <- with (M, sqrt(x^2 + y^2 +z^2))
p <- sin(2*R)/(R+1e-3)

slice3D(x, y, z, colvar = p, 
        xs = 0, ys = c(-4, 0, 4), zs = NULL)

isosurf3D(x, y, z, colvar = p, level = 0, col = "red")
@

\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}
\begin{center}
<<label=slice,fig=TRUE,echo=FALSE, width = 10, height = 6, pdf = FALSE, png = TRUE>>=
<<slice>>
@
\end{center}
\caption{Slices and isosurfaces from volumetric data}
\label{fig:slice}
\end{figure}

\newpage
\section{surf3D}
Function \code{surf3D} creates 3-D surface plots.
<<>>=
args(surf3D)
@
Here are 4 applications, showing the different options of coloration.
<<label=surf,include=FALSE>>=
  par(mfrow = c(2, 2), mar = c(0, 0, 0, 0)) 

 # Shape 1
  M  <- mesh(seq(0,  6*pi, length.out = 80), 
             seq(pi/3, pi, length.out = 80))
  u  <- M$x ; v <- M$y

  x <- u/2 * sin(v) * cos(u)
  y <- u/2 * sin(v) * sin(u)
  z <- u/2 * cos(v)

  surf3D(x, y, z, colvar = z, colkey = FALSE, box = FALSE)


 # Shape 2: add border
  M  <- mesh(seq(0, 2*pi, length.out = 80), 
             seq(0, 2*pi, length.out = 80))
  u  <- M$x ; v  <- M$y

  x  <- sin(u)
  y  <- sin(v)
  z  <- sin(u + v)

  surf3D(x, y, z, colvar = z, border = "black", colkey = FALSE)

 # shape 3: uses same mesh, white facets
  x <- (3 + cos(v/2)*sin(u) - sin(v/2)*sin(2*u))*cos(v)
  y <- (3 + cos(v/2)*sin(u) - sin(v/2)*sin(2*u))*sin(v)
  z <- sin(v/2)*sin(u) + cos(v/2)*sin(2*u)
 
  surf3D(x, y, z, colvar = z, colkey = FALSE, facets = FALSE)

 # shape 4: more complex colvar
  M  <- mesh(seq(-13.2, 13.2, length.out = 50), 
             seq(-37.4, 37.4, length.out = 50))
  u  <- M$x   ; v <- M$y

  b <- 0.4; r <- 1 - b^2; w <- sqrt(r)
  D <- b*((w*cosh(b*u))^2 + (b*sin(w*v))^2)
  x <- -u + (2*r*cosh(b*u)*sinh(b*u)) / D
  y <- (2*w*cosh(b*u)*(-(w*cos(v)*cos(w*v)) - sin(v)*sin(w*v))) / D
  z <- (2*w*cosh(b*u)*(-(w*sin(v)*cos(w*v)) + cos(v)*sin(w*v))) / D

  surf3D(x, y, z, colvar = sqrt(x + 8.3), colkey = FALSE, 
         border = "black", box = FALSE)
@

\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}
\begin{center}
<<label=surf,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<surf>>
@
\end{center}
\caption{Surface plots}
\label{fig:surf}
\end{figure}

\subsection{scatter2D and scatter3D}
Functions \code{scatter2D} and \code{scatter3D} draw scatterplots.

<<>>=
args(scatter2D)
args(scatter3D)
@
The dataset \code{quakes} is plotted using function \code{scatter3D}.
Before the 3-D quakes data are drawn, small dots are added on the bottom 
and on the depth plane (\code{panelfirst}).

<<label=scatter,include=FALSE>>=
par(mfrow = c(1, 1))
  panelfirst <- function(pmat) {
    zmin <- min(-quakes$depth)
    XY <- trans3D(quakes$long, quakes$lat, 
                  z = rep(zmin, nrow(quakes)), pmat = pmat)
    scatter2D(XY$x, XY$y, colvar = quakes$mag, pch = ".", 
            cex = 2, add = TRUE, colkey = FALSE)

 
    xmin <- min(quakes$long)
    XY <- trans3D(x = rep(xmin, nrow(quakes)), y = quakes$lat, 
                  z = -quakes$depth, pmat = pmat)
    scatter2D(XY$x, XY$y, colvar = quakes$mag, pch = ".", 
            cex = 2, add = TRUE, colkey = FALSE)
  }

   with(quakes, scatter3D(x = long, y = lat, z = -depth, colvar = mag, 
       pch = 16, cex = 1.5, xlab = "longitude", ylab = "latitude", 
       zlab = "depth, km", clab = c("Richter","Magnitude"),
       main = "Earthquakes off Fiji", ticktype = "detailed", 
       panel.first = panelfirst, theta = 10, d = 2, 
       colkey = list(length = 0.5, width = 0.5, cex.clab = 0.75))
       )
@

\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=scatter,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<scatter>>
@
\end{center}
\caption{Scatter plot}
\label{fig:scatter}
\end{figure}

\newpage
\subsection{arrows3D, arrows2D}
Functions \code{arrows2D} and \code{arrows3D} extend R function \code{arrows} with
a color variable.
<<label=arrows,include=FALSE>>=
par (mfrow = c(1, 2))
arrows2D(x0 = runif(10), y0 = runif(10), 
         x1 = runif(10), y1 = runif(10), colvar = 1:10, 
         code = 3, main = "arrows2D")

arrows3D(x0 = runif(10), y0 = runif(10), z0 = runif(10),
         x1 = runif(10), y1 = runif(10), z1 = runif(10),
         colvar = 1:10, code = 1:3, main = "arrows3D", colkey = FALSE)
@
\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=arrows,fig=TRUE,echo=FALSE, width = 10, height = 6, pdf = FALSE, png = TRUE>>=
<<arrows>>
@
\end{center}
\caption{arrows}
\label{fig:arrows}
\end{figure}

\newpage
\section{Functions based on image}
The \code{image2D} function is an extended version of \code{image}.
It has two S3 methods:
\begin{verbatim}
image2D(z =, ...)

image2D.matrix(z, x = NULL, y = NULL, ...,
             col = jet.col(100), NAcol = "white", facets = TRUE, 
             contour = FALSE, colkey = NULL, resfac = 1, 
             clab = NULL, theta = 0, border = NA) 
image2D.array(z, margin = c(1, 2), subset, ask = NULL, ...) 

\end{verbatim}

The data set \code{Oxsat} has oxygen saturation values in the ocean, at 2dg 
horizontal resolution, and for 33 depth intervals.
<<>>=
names(Oxsat)
dim(Oxsat$val)
@
Function \code{image2D.array} plots several depth intervals at once, looping over
the first and second margin. The color key is added in a separate figure.
<<label=image2D,include=FALSE>>=
sub <- c(1, 5, 9) 
image2D(z = Oxsat$val, subset = sub, 
        x = Oxsat$lon, y = Oxsat$lat,
        margin = c(1, 2), NAcol = "black", colkey = FALSE,
        xlab = "longitude", ylab = "latitude", 
        main = paste("depth ", Oxsat$depth[sub], " m"),
        clim = c(0, 115), mfrow = c(2, 2))
colkey(clim = c(0, 115), clab = c("O2 saturation", "percent"))        
@

\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=image2D,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<image2D>>
@
\end{center}
\caption{image2D function}
\label{fig:image2D}
\end{figure}

\newpage
\section{Composite figures}
It is also possible to make a composite figure combining several functions.
<<label=Composite,include=FALSE>>=
 persp3D(z = volcano, zlim = c(-60, 200), phi = 20,    
    colkey = list(length = 0.2, width = 0.4, shift = 0.15,
      cex.axis = 0.8, cex.clab = 0.85), lighting = TRUE, lphi = 90,
    clab = c("","height","m"), bty = "f", plot = FALSE)

# create gradient in x-direction
 Vx <- volcano[-1, ] - volcano[-nrow(volcano), ]

# add as image with own color key, at bottom 
 image3D(z = -60, colvar = Vx/10, add = TRUE, 
    colkey = list(length = 0.2, width = 0.4, shift = -0.15,
      cex.axis = 0.8, cex.clab = 0.85),
   clab = c("","gradient","m/m"), plot = FALSE)

# add contour  
 contour3D(z = -60+0.01, colvar = Vx/10, add = TRUE, 
    col = "black", plot = TRUE)
@
\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=Composite,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<Composite>>
@
\end{center}
\caption{Several color keys in composite figure}
\label{fig:Composite}
\end{figure}

\newpage
\section{plotting with legends and colorkeys}
As from version 1.4, a new feature has been added, allowing plotting functions to have colorkeys, or legends that are positioned outside of the plotting region.
See ?legendplot. 

For instance, to put a legend outside of a boxplot:

<<label=legendboxplot,include=FALSE>>=
pm <- par(mar = c(4,3,4,2))

legend.plt(formula = len ~ dose:supp, data = ToothGrowth,
          boxwex = 0.5, col = c("orange", "yellow"),
          main = "Guinea Pigs' Tooth Growth",
          xlab = "Vitamin C dose mg", ylab = "tooth length",
          sep = ":", lex.order = TRUE, ylim = c(0, 35), yaxs = "i", 
          method = "boxplot", legend.side = 2, 
          legend = list(legend = c("Ascorbic acid", "Orange juice"),
                          fill = c("yellow", "orange")))
@

\setkeys{Gin}{width=0.8\textwidth}
\begin{figure}
\begin{center}
<<label=legendboxplot,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<legendboxplot>>
@
\end{center}
\caption{adding a legend to boxplot}
\label{fig:legendboxplot}
\end{figure}

\newpage
To add a colorkey 

<<label=colorkeypie,include=FALSE>>=
n <- 100
colorkey.plt(method = "pie", x = rep(1, n), labels = "", 
             col = rainbow(n), border = NA,
             main = "colorkeyplot with 'pie'", 
             colorkey = list(col = rainbow(n), clim = c(1,n)))
@

\setkeys{Gin}{width=0.5\textwidth}
\begin{figure}
\begin{center}
<<label=colorkeypie,fig=TRUE,echo=FALSE, pdf = FALSE, png = TRUE>>=
<<colorkeypie>>
@
\end{center}
\caption{adding a colorkey}
\label{fig:colorkeypie}
\end{figure}


\newpage
\section{Issues}
\subsection{Specifying axes limits}
There are two ways in which the axes limits can be set, either allowing the 3D data
to overflow the axes, or to be clipped.
Consider the following code, which is based on a demo from the \pkg{rgl} package \citep{rgl}.
In this code, the \code{iris} data set is fitted with two models, and the 
data plotted using \code{scatter3D} and the fitted surfaces using \code{persp3D}.

The z-axis limits are set to \code{c(1, 9)}; plotting is postponed until all
objects have been added (this avoids clogging - see next subsection).
Then the result is plotted using \code{plotdev} (this plots the graph on any device
that has been opened). The first time, without specifying the axes limits, the 
second time setting the axes limits.
<<label=Fit,include=FALSE>>=
nout <- 30
xout <-  with(iris, seq(min(Sepal.Length), max(Sepal.Length), length = nout))
yout <-  with(iris, seq(min(Sepal.Width) , max(Sepal.Width),  length = nout))

xy  <- expand.grid(Sepal.Length = xout, Sepal.Width = yout)

# Fit two models, linear and quadratic
mod   <- with(iris, lm(Petal.Length ~Sepal.Length  + Sepal.Width))
mod2  <- with(iris, lm(Petal.Length ~Sepal.Length  + Sepal.Width +
                       I(Sepal.Length^2) + I(Sepal.Width^2) +
                       I(Sepal.Length*Sepal.Width)))
# prodict at new values
zpred.1 <- matrix(
   predict(mod, newdata = xy), nrow = nout, ncol = nout)
zpred.2 <- matrix(
   predict(mod2, newdata = xy), nrow = nout, ncol = nout)

# make graph, postpone plotting till the end
par(mfrow = c(1, 2))

with(iris,
  scatter3D(Sepal.Length, Sepal.Width, Petal.Length,
            colvar = as.numeric(Species), colkey = FALSE,
            col = c("blue", "red", "gold"), bty = "b",
            xlab = 'SL', ylab = 'PL', zlab = 'SW', zlim = c(1, 9),
            pch = 16, cex = 2, theta = 0, plot = FALSE))

persp3D(x = xout, y = yout, z = zpred.1, facets = NA,
  add = TRUE, col = "blue", plot = FALSE)

persp3D(x = xout, y = yout, z = zpred.2,
  add = TRUE, col = "green", plot = FALSE)

# plot using traditional device
plotdev(theta = -50, alpha = 0.5)
plotdev(theta = -50, alpha = 0.5, zlim = c(1, 9))

# if you want to see this in rgl:
# library(plot3Drgl)
#plotrgl(alpha = 0.5)
@
\setkeys{Gin}{width=1.0\textwidth}
\begin{figure}
\begin{center}
<<label=Fit,fig=TRUE,echo=FALSE, width = 10, height = 6, pdf = FALSE, png = TRUE>>=
<<Fit>>
@
\end{center}
\caption{calling plotdev() with axes limits causes the figure to be clipped}
\label{fig:hyps}
\end{figure}
\subsection{Preventing clogging of the figures}
When a lot of objects are \code{add}ed then the ultimate figure may appear clogged.
This is because each time something is added, all the objects are redrawn on top of 
what was already there - the only thing that is not redrawn are the axes and titles.
To create `slim' figures, use \code{plotdev()} (e.g. see previous section).



\section{Finally}

This vignette was made with Sweave \citep{Leisch02}.

\clearpage
\bibliography{vignettes}

\end{document}