% -*- fill-column: 80; comment-column: 50; -*-

\documentclass{article}

% \VignetteIndexEntry{Introduction to oce}
% \VignetteDepends{oce}
% \VignetteKeyword{oceanography}

\usepackage{url}
\usepackage{fullpage}
\usepackage{boxedminipage}
\usepackage{hyperref}
\usepackage{makeidx}
\usepackage{titlesec}
\makeindex

\parskip=1.5ex plus 1.5ex minus 1.25ex

\titleformat{\section}[block]{\normalfont\large\bfseries}{\thesection}{1em}{}
\titlespacing{\section}{0em}{2em plus 3em minus 0.5em}{0.15em plus 0.15em minus 0.125em}
\titleformat{\subsection}[block]{\normalfont\itshape}{\thesubsection}{1em}{}
\titlespacing{\subsection}{0em}{1em plus 2em minus 0.5em}{-0.15em plus 0.15em minus 0.125em}

\newcommand{\workedexercise}[2]{
	\vspace{2ex plus 2ex minus 1ex}
	\begin{boxedminipage}[c]{0.95\linewidth}
		{\textbf{Exercise #1}.\hspace{1em}#2}
	\end{boxedminipage}
	\vspace{2ex plus 2ex minus 1ex}
}
\newcommand{\workedanswer}[2]{
\goodbreak
\vskip 1.5ex plus 0.5ex minus 0.5ex
\noindent\textbf{Exercise #1 -- #2.}
}


\begin{document}

\title{The OCE package}
\author{Dan E. Kelley}
\maketitle


\begin{abstract}

The \verb@oce@ package makes it easy to read, summarize and plot data from a
variety of Oceanographic instruments, isolating the researcher from the quirky
data formats that are common in this field. It also provides functions for
working with basic seawater properties such as the equation of state, and with
derived quantities such as the buoyancy frequency.  Although simple enough to be
used in a teaching context, \verb@oce@ is powerful enough for a research
setting.  These things are illustrated here, in the context of some practical
examples.  Worked examples are provided, in order to help readers take early
steps towards using the \verb@oce@ package in their research.

\end{abstract}

\section{Introduction}

Oceanographers must deal with measurements made by a wide variety of
instruments, a task that is complicated by the delight instrument manufacturers
seem to take in inventing new data formats. The manufacturers often provide
software for scanning the data files and producing some standard plots, but this
software is of limited use to researchers who work with several instrument types
at the same time, and who need to carry the analysis beyond the first steps,
e.g. moving beyond engineering plots to scientific plots and to statistical
analysis.

\begin{figure}[ht]
  \begin{center}
    \includegraphics[width=0.5\hsize]{ctd-object}
  \end{center}
  \caption{\label{f:ctdobject}Sketch of the contents of a \texttt{ctd}
    object. The elements \texttt{data}, \texttt{metadata}, and
    \texttt{processing.log} are shared by all \texttt{oce} objects, although
    their contents vary from object to object.}
\end{figure}

The need to scan diverse data files was one motivation for the creation of
\verb@oce@, but an equal goal was to make it easy to work with the data once
they are in the system.  This was accomplished partly by the provision of
functions to work with the data, and partly by developing a uniform object
design that lets users reach inside without guesswork.

At the core of the \verb@oce@ design process is a policy of adding features
according to the priorities of practical research. As a result, \verb@oce@ is a
fairly comfortable tool today, and it should remain so as it grows.


\section{Object design}

As illustrated in Figure~\ref{f:ctdobject}, each \texttt{oce} object is a list
containing three elements: (a) \texttt{data}, a list or a data frame containing
the actual data (e.g., for a CTD object, this will contain pressure,
temperature, etc.), (b) \texttt{metadata}, a list containing data such things as
file headers, the location of a CTD cast, etc., and (c) \texttt{processing.log},
a list that documents how the file was created (often by a \texttt{read} or
\texttt{as} method) and how it was changed thereafter (e.g. by decimating a CTD
cast).

The uniformity of the various \texttt{oce} objects makes it easy to build skill
in examining and modifying objects.  A simple S3-based object inheritance scheme
permits somewhat generic processing of \texttt{oce} objects, e.g. \verb@plot(x)@
produces different sorts of plots for different sorts of \verb@x@ objects.


\section{Calculations of seawater properties}

\index{calculation!seawater properties}
\index{seawater properties, calculations of}

The \verb@oce@ package provides many functions for dealing with seawater
properties. Perhaps the most used is \verb@sw.rho(S,T,p)@, which computes
seawater density $\rho=\rho(S, T, p)$, where $S$ (PSU) is salinity, $T$
($^\circ$C) is \emph{in-situ} temperature, and $p$ (dbar) is seawater pressure,
expressed as an increment over standard atmospheric pressure.  (This
  and similar functions starts with the letters \verb@sw.@ to designate that
  they relate to seawater properties; a future version of \verb@oce@ may  provide
  air properties, with functions names starting with \verb@air.@) The result is a
number in the order of $1000$\,kg/m$^3$.  For many purposes, Oceanographers
prefer to use the density anomaly $\sigma=\rho-1000$\,kg/m$^3$, provided with
\verb@sw.sigma(S,t,p)@, or its adiabatic cousin $\sigma_\theta$, provided with
\verb@sw.sigma.theta(S,t,p)@.

\looseness=1
\index{dynamic height}
Most of the functions use the UNESCO formulations of seawater properties, but new formulations
may be added as they come into use in the literature.
A partial list of seawater functions is as follows:
\verb@sw.dynamic.height@ (dynamic height),
\verb@sw.N2@ (buoyancy frequency),
\verb@sw.S.C.T.p@ (salinity $S$ from conductivity $C$, $T$ and $p$),
\verb@sw.S.T.rho@ ($S$ from $T$ and $\rho$),
\verb@sw.T.S.rho@ ($T$ from $S$ and $\rho$),
\verb@sw.T.freeze@ (freezing temperature),
\verb@sw.alpha@ (thermal expansion coefficient $\alpha=-\rho_0^{-1}\partial\rho/\partial T$),
\verb@sw.beta@ (haline compression coefficient $\beta=\rho_0^{-1}\partial\rho/\partial S$),
\verb@sw.alpha.over.beta@ ($\alpha/\beta$),
\verb@sw.conductivity@ (conductivity from $S$, $T$ and $p$),
\verb@sw.depth@ (depth from $p$ and latitude),
\verb@sw.lapse.rate@ (adiabatic lapse rate),
\verb@sw.rho@ (density $\rho$ from $S$, $T$ and $p$),
\verb@sw.sigma@ ($\rho-1000$\,kg/m$^3$),
\verb@sw.sigma.t@ ($\sigma$ with $p$ set to zero and temperature unaltered),
\verb@sw.sigma.theta@ ($\sigma$ with $p$ set to zero and temperature altered adiabatically),
\verb@sw.sound.speed@ (speed of sound),
\verb@sw.specific.heat@ (specific heat),
\verb@sw.spice@ (a quantity used in double-diffusive research),
\verb@sw.theta@ (potential temperature $\theta$),
and
\verb@sw.viscosity@ (viscosity).
Details and examples are provided in the documentation of these functions.


\workedexercise{1}{(a) What is the density of a seawater parcel at pressure
$100$\,dbar, with salinity $34$\,PSU and temperature 10$^\circ$C?
(b) What temperature would the parcel have if raised adiabatically to the surface?
(c) What density would it have if raised adiabatically to the surface?
(d) What density would it have if lowered about 100m, increasing the pressure to $200$dbar?
(e) Draw a blank $T$-$S$ diagram with $S$ from $30$ to $40$\,PSU and $T$ from $-2$ to 20$^\circ$C.
(Answers are provided at the end of this document.)
}


\section{CTD data}
\subsection{Example with pre-trimmed data}

\index{data!ctd profile \verb@ctd@}

To get you started with CTD data, \verb@oce@ provides a sample data set that has
been trimmed to just the downcast portion of the sampling.  (See the next
section to learn how to do this trimming.).  The commands
<<label=ctdplot,eval=FALSE>>=
library(oce)
data(ctd)
plot(ctd)
@
produce Figure~\ref{fig:ctd}. You may also get a summary of the data with
<<ctdsummary,eval=false>>=
summary(ctd)
@

The object used to hold CTD data stores not just the data, but also the raw
header sequence, and whatever has been discovered about the dataset by parsing
the header; use
<<ctdnames,eval=false>>=
names(ctd)
@
to learn about the data use
<<ctdnames2,eval=false>>=
names(ctd$data)
@
to find out what sensors were attached to the instrument, thus providing data columns.

Of course, you may apply any \verb@R@ techniques to the data in \verb@oce@
objects, e.g. \verb@hist(ctd$data$temperature)@ would produce a histogram of
temperature for the \verb@ctd@ object.  It is always worth checking, though, to
see if \verb@oce@ has already defined a function that you may be applying,
e.g. \verb@plot.TS@ will produce a lovely temperature-salinity diagram, with
isopycnals and proper units on the axes.

The package provides facilities for some common operations with oceanographic
data, such as trimming CTD profiles with \verb@ctd.trim()@, but of course you
may do that sort of work by acting on the data directly, if necessary.  Just
make sure you realize that the metadata will not be altered if you do that.
Also, it is a good idea to add log entries to any objects that you change, by
using the \verb@processing.log.append()@ function.  (You can see an example of
this in action with \verb@?section@.)

\begin{figure}
\begin{center}
\setkeys{Gin}{width=\textwidth}
<<label=ctdfig, width=9, height=6, fig=TRUE, echo=FALSE>>=
<<ctdplot>>
@
\end{center}
\caption{Overview graph of the sample CTD dataset \texttt{ctd}, acquired in
the St Lawrence Estuary Internal Wave Experiment.  (This dataset has been
trimmed to the downcast; see the text and Figure~\ref{fig:ctdraw}.)}
\label{fig:ctd}
\end{figure}

\workedexercise{2}{Plot a profile of $\sigma_\theta$ and $N^2$, for just the data in the pycnocline.  Do this with both \texttt{ctd.trim} and \texttt{subset}.}

\subsection{Example with raw data}

Practicing Oceanographers may be wondering how the CTD cast used in the
preceding section was trimmed of equilibration-phase and upcast-phase data. Data
from these sections are spurious and must be trimmed as a first step in
processing. For example, consider the following code.
<<label=ctdrawplot,eval=FALSE>>=
data(ctd.raw)
plot.ctd.scan(ctd.raw)
@
\begin{figure}
\begin{center}
\setkeys{Gin}{width=0.55\textwidth}
<<label=ctdrawfig,fig=TRUE,width=6, height=4, echo=FALSE>>=
<<ctdrawplot>>
@
\end{center}
\caption{Scanwise plot of the \texttt{ctd.raw} sample data set.  Note the wild
spike at the start, the equilibration phase before the downcast, and the
spurious freshening signal near the start of the upcast.  See the text for a
discussion of how inspection of such graphs can help in trimming CTD data.}
\label{fig:ctdraw}
\end{figure}

\noindent This produces a two-panel plot (Figure~\ref{fig:ctdraw}) of the data
as a time-series, revealing not just the (useful) downcast, but also the
subsequent upcast sequence.  The x-axis in this plot is the scan number, which
is a convenient index for extraction of the downcast portion of the profile by
an essentially manual method, e.g. proceeding with a sequence of commands such
as
<<ctdscaneg,eval=false>>=
plot.ctd.scan(ctd.trim(ctd.raw, "scan", c(140,250)))
plot.ctd.scan(ctd.trim(ctd.raw, "scan", c(150,250)))
@
This is the ``gold standard'' method, which is recommended for detailed
work. However, for quick work, you may find that the automatic downcast
detection scheme works adequately, e.g.
<<ctdtrimeg,eval=false>>=
ctd.trimmed <- ctd.trim(ctd.raw)
@

It should be noted that \verb@ctd.trim@ inserts entries into the object's log
file, so that you (or anyone else with whom you share the object) will be able to
see the details of how the trimming was done.

\index{reading!ctd profile}

Once the profile has been trimmed, you may wish to use \texttt{ctd.decimate()}
to smooth the data and interpolate the smoothed results to uniformly-spaced
pressure values. For example, a quick examination of a file might be done by the
following:
<<ctddectrim,eval=false>>=
plot(ctd.decimate(ctd.trim(read.ctd("stn123.cnv"))))
@

\subsection{Example with WOCE archive data}

The package has a harder time scanning the headers of data files in the WOCE
archive format than it does in the Seabird format illustrated in the previous
examples. This is mainly because front-line researchers tend to work in the
Seabird format, and partly because the WOCE format is odd. For example, the
first line of a WOCE file is of the form \texttt{CTD,20060609WHPOSIODAM} (or
\texttt{BOTTLE,...}).  Scanning the item to the left of the comma is not
difficult (although there are variants to the two shown, e.g. \texttt{CTDO}
sometimes occurs).  The part to the right of the comma is more difficult.  The
first part is a date (yyyymmdd) so that's no problem.  But then things start to
get tricky.  In the example provided, this text contains the division of the
institute (WHPO), the institute itself (SIO), and initial of the investigator
(DAM). The problem is that no dividers separate these items, and that there seem
to be no standards for the item lengths. Rather than spend a great deal of time
coding special cases (e.g. scanning to see if the string \texttt{WHOI} occurs in
the header line), the approach taken with \verb@oce@ is to ignore such issues
relating to quirky headers.  This frees up time to work on more important
things, such as plotting the data.

Of course, R provides access to object constituents, so that you are
free to do such things as
<<ctdfix1,eval=false>>=
x <- read.ctd("nnsa_00934_00001_ct1.csv", type="WOCE")
x$metadata$institute <- "SIO" # bad way
x$metadata$scientist <- "DAM" # bad way
@
but it is \emph{bad practice} to alter metadata in such a way, because
doing so alters the object without altering its documentation.  The
preferred scheme is to do as follows.
\index{changelogs in oce objects}
\index{oce object changelog}
\index{logged changes to oce objects}
\index{editing oce objects}
\index{oce.edit}
<<ctdfix2,eval=false>>=
x <- read.ctd("nnsa_00934_00001_ct1.csv", type="WOCE")
x <- oce.edit(x, "institute", "SIO") # better way
x <- oce.edit(x, "scientist", "DAM") # better way
@
which will store a note about the changes in the object's log.  Even
better, provide a reason for the change, and sign the change with your
name:
<<ctdfix3,eval=false>>=
x <- read.ctd("nnsa_00934_00001_ct1.csv", type="WOCE")
x <- oce.edit(x, "institute", "SIO", "human-parsed","Dan Kelley") # best way
x <- oce.edit(x, "scientist", "DAM", "human-parsed","Dan Kelley") # best way
@

For a real-world example (with warts!), visit
\url{http://cchdo.ucsd.edu/data_access?ExpoCode=58JH199410} and download the zip
file containing the Arctic section called ``CARINA'', measured in 1994. Expand
the zip file, enter the directory, and run the code below.
<<arcticeg,eval=false>>=
library(oce)
# Source: http://cchdo.ucsd.edu/data_access?ExpoCode=58JH199410
files <- system("ls *.csv", intern=TRUE)
for (i in 1:length(files)) {
    cat(files[i], "\n")
    x <- read.ctd(files[i])
    if (i == 1) {
        plot.TS(x, xlim=c(31, 35.5), ylim=c(-1, 10), type='l', col="red")
    } else {
        lines(x$data$salinity, x$data$temperature, col="red")
    }
}
@

What you'll see is an overall $T$-$S$ diagram for the entire dataset. It may take
a while, since the dataset contains over 90,000 observations. You may note that,
even though this is an official, quality-controlled dataset, it is not without
problems. The graph that is produced by this code has several spurious lines
oriented horizontally (indicating spurious salinity) and vertically (indicating
spurious temperature). One way to find such values is to put the lines
<<rangeg,eval=false>>=
print(range(x$data$temperature))
print(range(x$data$salinity))
@
after the \verb@read.ctd()@ command. One thing you'll find is that station 987
has a minimum salinity range of 0.0009 to 987. These values are clearly in
error, as are the temperatures at this spot in the file. (It is perhaps
revealing that the spurious salinity is equal to the station number.) Indeed, at
this spot in the file it can be seen that the pressure jumps from 1342 to 0, and
then starts increasing again; the file contains two profiles, or the same
profile twice. This is not the only flaw that is revealed by the graph, and by
\verb@range@ commands; a generous user would spend a week tracking down such
issues, and would then contact the data provider (or the chief scientist of the
field work) with specific suggestions for correcting the files. The point here
is to highlight how this package can be used with real-world data.

\subsection{Section plots}
The commands
<<label=sectionplot,eval=FALSE>>=
data(section)
data(coastline.hal)
plot(section, coastline=coastline.hal)
@
will plot a summary diagram containing sections of
$T$, $S$, and $\sigma_\theta$, along with a chart indicating station
locations. In addition to such overview diagrams, \verb@plot@ can also create
individual plots of individual properties.

\workedexercise{3}{Draw a $T$-$S$ diagram for the section data, colour-coded by station}

\index{reading!ctd section}
\index{data!section \verb@a03@, North Atlantic along 36$^\circ$N}

\index{Gulf Stream}

The Halifax section is supplied in a pre-gridded format, but some datasets have
different pressure levels at each station.  For such cases, the
\verb@section.grid@ function may be used, e.g.
<<label=sectionplota03,eval=FALSE>>=
# File source -- http://cchdo.ucsd.edu/data_access?ExpoCode=90CT40_1
# a03 <- read.section("a03_hy1.csv")
data(a03)
Gulf.Stream <- subset(a03, indices=124:102)
Gulf.Stream.gridded <- section.grid(Gulf.Stream, p=seq(0, 1600, 25))
data(coastline.world)
plot(Gulf.Stream.gridded, coastline=coastline.world, map.xlim=c(-80,-60))
@
produces Figure~\ref{fig:sectiona03}.  The ship doing the sampling was
travelling westward from the Mediterranean towards North America, taking 124
stations in total; the \verb@station.indices@ value selects the last few
stations of the section, during which the ship heading was changed to run in a
northwesterly direction, to cross isobaths (and perhaps, the Gulf Stream) at
right angles.

\begin{figure}
\begin{center}
\setkeys{Gin}{width=1\textwidth}
<<label=sectionfiga03,fig=TRUE,width=9, height=5, echo=FALSE>>=
<<sectionplota03>>
@
\end{center}
\caption{\label{fig:sectiona03}
Portion of the CTD section designated A03, showing the Gulf Sream region.
}
\end{figure}

\workedexercise{4}{Plot dynamic height and geostrophic velocity across the Gulf Stream.}

\section{Coastline and topographic data}

\index{data!coastline}
\index{data!topography}

Coastline data are available from a variety of sources. The NOAA site
\url{http://www.ngdc.noaa.gov/mgg_coastline/} is particularly popular, and it
has the advantage of providing data in Splus format.  The function
\verb@read.coastline@ can handle reading that format (plus some other formats),
and \verb@plot@ on the resulting object will produce a simple coastline map.
The only real advantage over plotting things yourself is that latitude and
longitude are scaled to give natural shapes near the centre of the plot.

Bathymetric charts, or more generally topographic maps, can be produced easily.
A sample data set is provided, so that
<<label=topoplot,eval=FALSE>>=
library(oce)
data(topo.maritimes)
plot(topo.maritimes, xlim=c(-66, -58), ylim=c(44, 50),
     water.z=c(  -50, -100, -150, -200, -300, -400, -500, -1000, -2000),
     water.lwd=c(  1,    1,    1,    1,    1,   1,   1.5,   1.5,   1.5))
@
will produce a chart of the waters embraced by the Maritime Provinces of Canada
(Figure~\ref{fig:topo}).
\begin{figure}
\begin{center}
\setkeys{Gin}{width=0.7\textwidth}
<<label=topofig,fig=TRUE,width=7, height=5, echo=FALSE>>=
<<topoplot>>
@
\end{center}
\caption{Topography of eastern Canada, centred on Les \^{i}les de la Madeleine.}
\label{fig:topo}
\end{figure}

\section{Sea-level data}

\subsection{Time-domain analysis}

\index{sea level!during Hurricane Juan}
\index{Hurricane Juan!surge seen in time-series of sea level}

The commands
<<label=sealevelplot,eval=FALSE>>=
library(oce)
#sealevel <- read.oce("../../tests/h275a96.dat")
data(sealevel.hal)
plot(sealevel.hal)
@
load and graph a build-in dataset of sea-level timeseries. The result, shown in
Figure~\ref{fig:sealevel}, is a four-panel plot. The top panel is a timeseries
view that provides an overview of the entire data set. The second panel is
narrowed to the most recent month, which should reveal spring-neap cycles if the
tide is mixed. The third panel is a spectrum, with a few tidal constituents
indicated. At the bottom is a cumulative spectrum, which makes these
narrow-banded constituents quite visible.

\begin{figure}
\begin{center}
\setkeys{Gin}{width=\textwidth}
<<label=sealevelfig,fig=TRUE,width=7,height=5,echo=FALSE>>=
<<sealevelplot>>
@
\end{center}
\caption{Sea-level timeseries measured in 2003 in Halifax Harbour.  (The spike
in September is the storm surge associated with Hurricane Juan, regarded by the
Canadian Hurricane Centre to be one of the most powerful and damaging hurricanes
to ever hit Canada.)}
\label{fig:sealevel}
\end{figure}

\workedexercise{5}{Illustrate Halifax sealevel variations during Hurricane Juan.}

\workedexercise{6}{Draw a spectrum of sea-level variation, with the M2 tidal component indicated.}

\subsection{Tidal analysis}

In a future version, tidal analysis will be provided, along the lines of the
t-tide package in Matlab.  A preliminary version of tidal analysis is provided
by the \verb@tidem@ function provided in this version of the package, but
readers are cautioned that the results are certain to change in a future
version.  (The problems involve phase, and the inference of satellite nodes.)

%%% \index{sea level!storm surge during Hurricane Juan}
%%% \index{Hurricane Juan!surge detection by detiding}
%%%
%%% <<label=tideplot,eval=FALSE>>=
%%% library(oce)
%%% data(sealevel.hal)
%%% tide <- tidem(sealevel.hal)
%%% days <- 15 # focus on these days at end of series
%%% n <- length(sealevel.hal$data$elevation)
%%% i <- seq(n-24*days, n)
%%% t   <- sealevel.hal$data$t[i]
%%% eta <- sealevel.hal$data$elevation[i]
%%% eta.pred <- predict(tide)[i]
%%% plot(t, eta, type="l", ylim=c(-0.5,3))
%%% abline(h=0, col="pink")
%%% lines(t, eta - eta.pred, col="red")
%%% @
%%% produce Figure~\ref{fig:tide}, illustrating the use of tidal analysis to detect
%%% storm surges.  Note the use of the \verb@model@ that is returned by the tide
%%% analysis.
%%%
%%% (It should be noted that the tidal analysis portion of \verb@oce@ is still in
%%% active development, and the features are subject to change.  Use
%%% \verb@?tidem@ to learn more about this.)
%%%
%%% \begin{figure}
%%% 	\begin{center}
%%% \setkeys{Gin}{width=0.8\textwidth}
%%% <<label=tidefig,fig=TRUE,width=8,height=5,echo=FALSE>>=
%%% <<tideplot>>
%%% @
%%% \end{center}
%%% \caption{Sea-level variation in Halifax Harbour in 2003 (black line) and the residual
%%% after removal of tides using \texttt{tidem()}.  The spike on September 28th
%%% is a storm surge caused by Hurricane Juan, which caused a great deal of damage
%%% to the region.}
%%% \label{fig:tide}
%%% \end{figure}
%%%
%%% \goodbreak

\section{Temperature-Depth Recorder data}
The commands
<<label=ptplot,eval=FALSE>>=
library(oce)
data(pt)
plot(pt)
@
produce a plot (Figure~\ref{fig:pt}) of pressure-temperature recorder
(pt) measurements made during St Lawrence Estuary Internal Wave
Experiment, a CFCAS-funded research program in which the \texttt{oce}
author was a principal investigator.  The device was mounted on
a mooring that tilted significantly with the tide, yielding
strong tidal signals in pressure as well as temperature, thus
providing a rough indication of the temperature  profile.

This graph is trimmed of the in-air data; access to the raw data is
part of a worked example.

\begin{figure}
\begin{center}
\setkeys{Gin}{width=\textwidth}
<<label=ptfig ,fig=TRUE, width=7, height=5, echo=FALSE>>=
<<ptplot>>
@
\end{center}
\caption{Measurements made with a pressure-temperature recorder in the
  St Lawrence Estuary.  Note that the record contains in-air
  measurements made during an hour prior to lowering into the water,
  and during several days of transit from the field site back to
  Dalhousie University.  The text suggests a procedure for isolating
  in-water data, and removing the approximately 10 dbar atmospheric
  component of pressure.}
\label{fig:pt}
\end{figure}

\workedexercise{7}{Trim the in-air measurements from this \texttt{pt} record, and
then remove the atmospheric pressure from the signal.}


\section{Acoustic Doppler Current Profiler data}
\index{ADP}
The commands
<<label=adpplot,eval=FALSE>>=
library(oce)
data(adp)
e <- adp.xyz2enu(adp.beam2xyz(adp))
plot(e, which=1:2,
     adorn=expression({
       lines(x$data$ts$time, x$data$ts$pressure, lwd=3, col='green')
     }))
@
produce a plot (Figure~\ref{fig:adp}) of ADP measurements made
during St Lawrence Estuary Internal Wave Experiment.
\begin{figure}
\begin{center}
\setkeys{Gin}{width=\textwidth}
<<label=adpfig ,fig=TRUE, width=7, height=7, echo=FALSE>>=
<<adpplot>>
@
\end{center}
\caption{Measurements made with an ADP in the
  St Lawrence Estuary.  The data were recorded in beam format, so a two-step
  process is used to convert to enu coordinates.  The green curve indicates
  the water surface, inferred from the pressure measured by the ADP.}
\label{fig:adp}
\end{figure}



\section{The future of oce}

The present version of \verb@oce@ can only handle data types that the author has
been using lately in his research. New data types will be added as the need
arises in that work, but the author would also be happy to add other data types
that are likely to prove useful to the Oceanographic community.  (The data types
need not be restricted to Physical Oceanography, but the author will need some
help in dealing with other types of data, given his research focus.)

As for algorithms, there are plenty of gaps in \verb@oce@.  Dealing
with measurements of turbulence is a high priority for the author, for
example, as is improving the functions for processing data from
pressure-temperature sensors and acoustic Doppler current profilers.

Two principles will guide the addition of data types and functions:
(a) the need, as perceived by the author or by other contributors and
(b) the ease with which the additions can be made.

\section{Development website}

The site \url{http://code.google.com/p/r-oce/} provides a window on the
development that goes on between the CRAN releases of the package. Please visit
the site to report bugs, to suggest new features, or just to see how \verb@oce@
development is coming along.

Oce versions are only made official, i.e. released to the CRAN package system,
when they are considered to be of a high quality.  However, you may also get the
bleeding-edge version from the development website, perhaps to see whether the
author has successfully addressed a deficiency that you reported.  You will need
to have subversion installed on you machine, and then can do e.g.

\verb@svn checkout http://r-oce.googlecode.com/svn/trunk/ -r 835 r-oce-read-only@

\noindent to retrieve the intermediate version number 835. Naturally, you should
rename the directory from \verb@r-oce-read-only@ to \verb@oce@ before attempting
to check, build, or install the package.

%\begin{center}
%\vspace{2cm}
%
%------
%
%\vspace{2cm}
%
%See the end of the document for answers to exercises.
%
%\end{center}
%
%\newpage


\goodbreak

\section*{Answers to exercises}

\workedanswer{1}{Seawater properties}
<<>>=
library(oce)
sw.rho(S=34, t=10, p=100)
sw.theta(S=34, t=10, p=100)
sw.rho(S=34, t=sw.theta(S=34, t=10, p=100), p=0)
sw.rho(S=34, t=sw.theta(S=34, t=10, p=100, pref=200), p=200)
plot.TS(as.ctd(c(30,40),c(-2,20),rep(0,2)), grid=TRUE, col="white")
@

\workedanswer{2}{Profile plots}
Although one may argue as to the limits of the pycnocline, for illustration let us say it is in 5bar to 12dbar range.  One way to do this is
<<>>=
library(oce)
data(ctd)
pycnocline <- ctd.trim(ctd, "pressure", c(5,12))
plot.profile(pycnocline, type="density+N2")
@
and another is
\index{subset!extracting a portion of a profile}
<<>>=
library(oce)
data(ctd)
pycnocline <- subset(ctd, 5<=pressure & pressure<=12)
plot.profile(pycnocline, type="density+N2")
@

\index{section!extracting profile data from}
\workedanswer{3}{TS diagram for section data}
<<>>=
library(oce)
data(section)
SS <- TT <- pp <- id <- NULL
n <- length(section$data$station)
for (i in 1:n) {
	stn <- section$data$station[[i]] # save typing
	SS <- c(SS, stn$data$salinity)
	TT <- c(TT, stn$data$temperature)
	pp <- c(pp, stn$data$pressure)
	id <- c(id, rep(i, length(stn$data$pressure)))
}
ctd <- as.ctd(SS, TT, pp)
plot.TS(ctd, col=hsv(0.7*id/n), cex=2, pch=21)
@

\workedanswer{4}{Gulf Stream}
(Try \verb@?sw.dynamic.height@ for hints on smoothing.)
\index{Gulf Stream!geostrophic calculation}
\index{subset!extracting stations from a section}
<<>>=
library(oce)
data(a03)
Gulf.Stream <- subset(a03, indices=124:102)
dh <- sw.dynamic.height(Gulf.Stream)
par(mfrow=c(2,1))
plot(dh$distance, dh$height, type='b', xlab="", ylab="Dyn. Height [m]")
grid()
# 1e3 metres per kilometre
f <- coriolis(Gulf.Stream$data$station[[1]]$metadata$latitude)
g <- gravity(Gulf.Stream$data$station[[1]]$metadata$latitude)
v <- diff(dh$height)/diff(dh$distance) * g / f / 1e3
plot(dh$distance[-1], v, type='l', col="blue", xlab="Distance [km]", ylab="Velocity [m/s]")
grid()
abline(h=0)
@
\workedanswer{5}{Halifax sealevel during Hurricane Juan.}
\index{Hurricane Juan!worked example of sea-level plot}
\index{sea level!during Hurricane Juan}
A web search will tell you that Hurricane Juan hit about midnight, 2003-sep-28.
The author can verify that the strongest winds occurred a bit after midnight -- that was the time
he moved to a room without windows, in fear of flying glass.
<<>>=
library(oce)
data(sealevel.hal)
# Focus on 2003-Sep-28 to 29th, the time when Hurricane Juan caused flooding
plot(sealevel.hal,which=1,xlim=as.POSIXct(c("2003-09-24","2003-10-05"), tz="UTC"))
abline(v=as.POSIXct("2003-09-29 04:00:00", tz="UTC"), col="red")
mtext("Hurricane\nJuan", at=as.POSIXct("2003-09-29 04:00:00", tz="UTC"), col="red")
@

\workedanswer{6}{Sealevel spectrum}
Notice the use of \verb@(object)$data$(item)@ here.  All \verb@oce@ objects are lists,
and all of them contain a \verb@data@ element of a similar form to this.
<<>>=
library(oce)
data(sealevel.hal)
spectrum(sealevel.hal$data$elevation, spans=c(3,7))
abline(v=1/12.42)
mtext("M2",at=1/12.42,side=3)
@

\workedanswer{7}{Pressure-temperature Recorder plot}
\index{subset!using with pt data}
\index{data!TDR}

The accuracy of automatic removal of the atmospheric component done with
\verb@pt.trim()@ is often quite good, but data are so valuable that it always
makes sense to investigate the data by eye, in order to select a time window.
Manual selection of a time window is also important when a researcher is dealing
with a set of instruments whose measuremnts are to be combined and so need a
uniform time base (e.g. to make a matrix to be contoured).

In the code, the automatic method is used first, and then a series of time
indicators is drawn on the plot, to determine proper cutoff times.

Hint: be careful on timezones!
<<>>=
library(oce)
data(pt)
# First, get an idea for the data
plot(pt, which=c(1,3))

pt.trim.auto <- pt.trim(pt)
plot(pt.trim.auto)
print(range(pt.trim.auto$data$time)) # gives a starting point

abline(v=as.POSIXct("2008-06-25 00:00:00"),col="red")
abline(v=as.POSIXct("2008-06-26 00:00:00"),col="red")
abline(v=as.POSIXct("2008-06-25 15:00:00"),col="red")

abline(v=as.POSIXct("2008-07-04 00:00:00"),col="blue")
abline(v=as.POSIXct("2008-07-05 10:00:00"),col="blue")
abline(v=as.POSIXct("2008-07-04 12:00:00"),col="blue")
pt.trim <- subset(pt, as.POSIXct("2008-06-25 15:00:00") <= time &
                  time <= as.POSIXct("2008-07-04 12:00:00"))
print(pt.patm(pt)) # take 10.2dbar from this
pt.trim <- oce.edit(pt.trim, action="x$data$pressure <- x$data$pressure - 10.2", reason="remove atm pressure", person="DK")
plot(pt.trim)
@

\printindex
\end{document}
