From 11bc6913edd46e125f8a7def27c9f053617d2f60 Mon Sep 17 00:00:00 2001 From: dankelley Date: Thu, 7 Jun 2018 09:17:30 -0300 Subject: [PATCH] rebuild vignettes --- inst/doc/adp.html | 494 +++++++---------------- inst/doc/flags.html | 410 +++++-------------- inst/doc/oce.html | 944 +++++++++++++++++--------------------------- 3 files changed, 609 insertions(+), 1239 deletions(-) diff --git a/inst/doc/adp.html b/inst/doc/adp.html index f49c5375a4..635ba9b127 100644 --- a/inst/doc/adp.html +++ b/inst/doc/adp.html @@ -12,7 +12,7 @@ - + Acoustic Doppler profilers @@ -20,256 +20,46 @@ - + @@ -280,7 +70,7 @@

Acoustic Doppler profilers

Dan Kelley

-

2018-06-03

+

2018-06-07

@@ -307,32 +97,32 @@

2018-06-03

Abstract. This vignette explains the basics of using the oce package to handle measurements made by acoustic Doppler profilers. As in the related vignette about working with data-quality flags, the material is organized to reflect typical workflow. See the main oce vignette for general information about oce. The sample code provided here relies on the oce package, loaded as follows.

- +
library(oce)

1 Reading adp files

Acoustic Doppler profiling instruments are called ADCPs by some manufacturers and ADPs by others. The class name adp is used in oce, to handle both types. It is important to realize that there are many data formats for data of this type, and that the files are normally in a dense binary format that cannot be read with a text editor. It is necessary to use software to work with such files. Although oce can handle many file formats, and although its authors do their best to keep up with changes to file formats (at least for popular instruments), users will sometimes be forced to use manufacturer software to read the data, or at least convert to a format that can be manually read into R, and then converted to an adp object using the as.adp() function.

Suppose f is a character string naming a data file. In many cases,

- +
d <- read.oce(f)

will read the file, because read.oce will look inside the file to try to guess the format. If this doesnt work, one ought to try e.g.

- +
d <- read.adp(f)

which, again, will try to determine which subtype of adp file has been provided. If that fails, it might help to try a specialized function, such as read.adp.rdi() for an RDI file, read.adp.nortek() for a Nortek file, or read.adp.sontek() for a Sontek file. (Actually, each of these can handle a variety of sub-formats).

If any of these read. functions is called with just a single argument naming the file, then the entire dataset will be read in. This can be slow for large files (e.g. taking of order 0.1s per Mb), and so it can help to use the by argument, e.g. for a 138Mb file on the author’s computer,

- +
f <- "/data/archive/sleiwex/2008/moorings/m09/adp/rdi_2615/raw/adp_rdi_2615.000"
+dall  <- read.oce(f)

takes 15.2s of user time, but

- +
d100 <- read.oce(f, by=100)

takes 0.2s. It is also possible to window the data, using the from and to arguments, which (like by) can be integers that count ensembles (or “profiles”) or a POSIX time (see ?read.adp). For example, the adp dataset provided by oce was read with

- +
read.oce(f, from=as.POSIXct("2008-06-26", tz="UTC")`,
+            to=as.POSIXct("2008-06-27", tz="UTC"),
+            by="60:00"
+            latitude=47.88126, longitude=-69.73433)

with f as defined above. This also illustrates that it is possibly to specify the location of an instrument; the further processing that was done to create data(adp) will be explained in due course.

2 Summarizing adp data

The generic function summary provides a useful overview of an adp data object. For example,

- +
data(adp)
+summary(adp)

(results not shown the goal is to encourage the user to try this!) shows the sort of information normally present in an RDI file. Other manufacturer types will have somewhat different summaries, because both the metadata and the data differ from instrument to instrument. In all cases, there will be a vector holding time (named time), another holding the distance from instrument to the centre of a measurement bin, measured along a line equidistant between the beams (named distance), along with an array holding velocity (named v). (Other arrays may include, as in this RDI case, q for data quality, a for backscatter amplitude, and g for percentage goodness.)

Note that the summary reveals that instrument had been set up to record in beam coordinates, but that processing had been done to convert to enu (east-north-up) coordinates.

Exercise 1. Determine the name of the metadata item that holds the coordinate system of an adp object.

@@ -341,16 +131,16 @@

2 Summarizing adp data

3 Changing coordinate systems

It is critical for analysts to be aware of the coordinate system of an adp dataset, simply because those systems are so different. In what oce calls the beam system, velocities in each of the acoustic beams are measured along the slant angle of the beam. These data may be combined (using a transformation matrix which is part of the metadata displayed with the summary() function) into what oce calls a xyz system that is Cartesian and fixed to the instrument. Since instruments record heading, pitch, and roll, the xyz velocities can be converted easily to enu velocities.

In some applications, the scientific focus is primarily on the enu data, but it is a good idea to calculate and study the data at all stages of processing. Thus, although the toEnu() function will convert data measured in any coordinate system into an enu variant, many analysts will do something along the lines of

- +
beam <- read.oce(f)
+xyx <- beamToXyz(beam)
+enu <- xyzToEnu(xyz, declination=-18.1)

where the declination illustrated is the value used in creating the data(adp) dataset. (Compensating for declination is necessary because compasses are calibrated at a particular place and time, and compasses have no way of “knowing” the local angle between true north and magnetic north.)

4 Plotting

The generic plot() function is specialized to handle adp objects in a variety of simple but helpful ways. The type of plot is determined by the which argument. The default for which is to create a multi-panel image plot that shows time-space dependence for each velocity component. The number of panels depends on the number of beams used by the instrument, and each one has a label indicating the component name. For example,

- -

+
plot(adp)
+

has velocity panels named east, north and up, along with one named error, which is an estimate of the velocity-inference error.

Options to this plot variant can be used to control how distance from instrument is shown (which can be useful in visualizing data from upward- and downward-aligned instruments), what colours are used to represent the signals, etc.

Dozens of other plot types are also provided; see ?"plot,adp-method" for details – but be warned that a complex plot requires many arguments!

@@ -359,8 +149,8 @@

4 Plotting

5 Subsetting

Adp data can be subsetted in a wide variety of ways, e.g.

- -

+
plot(subset(adp, distance<20))
+

plots only data within 20m of the instrument. See ?"subset,adp-method" for more on subsetting.

Exercise 3. Plot the adp data for the first half of the sampling interval.

@@ -368,13 +158,13 @@

5 Subsetting

6 Accessing data within adp objects

All oce objects share a common framework, with so-called “slots” for metadata, data, and processing log. It is possible to examine, and modify, the elements of each of these slots using base R constructs, but it is usually better to use oce accessor functions. The details are provided with ?"[[,adp-method" and ?"[[<-,adp-method", but a few examples should suffice to sketch the general pattern.

For example,

- +
time <- adp[["time"]]
+distance <- adp[["distance"]]

extracts the times and distances of the profiles within data(adp), while the array holding the velocities is extract with

- +
v <- adp[["v"]]

Altering values uses the [[ to the left of the assignment operator, but this is not commonly required in simple work, so it is best for readers to focus first on accessing data. The exercises may help with this.

Some of the fields contained in adp objects, particularly the amplitude, correlation, and percent good fields (a, q, and g) are stored as raw-type. This is done primarly to save memory space, as expanded all of those fields into numeric types can be slow and they get large quickly. However, to actually work with the numbers, users can extract a “numeric” version by supplying the second argument to the [[ function, e.g.

- +
a <- adp[["a", "numeric"]]

This has the advantage of returning an object with the same matrix dimensions as the original (note that functions like as.numeric() typically discard dimension information, thereby converting an array to a vector).

Exercise 4. Plot a time series of mid-depth and depth-averaged eastward velocities.

Exercise 5. Perform an eigen analysis of eastward and northward velocity.

@@ -387,123 +177,123 @@

7.1 Solution 1

Exercise 1. Determine the name of the metadata item that holds the coordinate system of an adp object.

Solution.

With adp as loaded in the text, we may the names in the metadata with:

- +
sort(names(adp[["metadata"]]))

(results not shown), after which a keen eye will notice that two elements relate to coordinates:

- +
adp[["originalCoordinate"]]
+#> [1] "beam"
+adp[["oceCoordinate"]]
+#> [1] "enu"

and comparison with the output of summary(adp) reveals that the former is the coordiante system in which the data were originally measured, while the latter is the coordinate system of the transformed data. Using

- +
processingLogShow(adp)
+#> * Processing Log
+#>     - 2018-04-25 13:43:11.586 UTC: `read.oce("/data/archive/sleiwex/2008/moorings/m09/adp/rdi_2615/raw/adp_rdi_2615.000", ...)`
+#>     - 2018-04-25 13:43:11.589 UTC: `beamToXyzAdp(x = beam)`
+#>     - 2018-04-25 13:43:11.612 UTC: `xyzToEnu(x, declination=-18.1, debug=0)`

reveals the steps in the data transformation.

7.2 Solution 2

Exercise 2. Create a u-v scatterplot.

Solution.

- -

+
plot(adp, which="uv")
+

7.3 Solution 3

Exercise 3. Plot the adp data for the first half of the sampling interval.

Solution.

- -

+
plot(subset(adp, time < median(adp[["time"]])))
+

7.4 Solution 4

Exercise 4. Plot a time series of mid-depth and depth-averaged eastward velocities.

Solution.

- -

+
time <- adp[["time"]]
+v <- adp[["v"]]
+# The second index is for bin number, the third for beam number
+midIndex <- dim(v)[2]/2
+eastMid <- v[, midIndex, 1] # third index is beam
+distance <- adp[["distance"]][midIndex]
+oce.plot.ts(time, eastMid, ylab="Eastward velocity [m/s]")
+## Depth mean; note that na.rm, is passed by apply() to mean()
+eastMean <- apply(v[,,1], 1, mean, na.rm=TRUE)
+lines(time, eastMean, col=2)
+legend("top", lwd=1, col=1:2,
+       legend=c(paste("At", distance, "m"), "Depth Avg"))
+

Here, oce.plot.ts() is used instead of the basic R function plot(), because it provides a margin note on the detailed time interval.

7.5 Solution 5

Exercise 5. Perform an eigen analysis of eastward and northward velocity.

Solution.

- +
u <- adp[["v"]][,,1]
+v <- adp[["v"]][,,2]
+ok <- is.finite(u) & is.finite(v) # remove NA values
+u <- u[ok]
+v <- v[ok]
+eigen(cov(data.frame(u, v)))
+#> eigen() decomposition
+#> $values
+#> [1] 0.48350161 0.01557049
+#> 
+#> $vectors
+#>           [,1]       [,2]
+#> [1,] 0.5441538 -0.8389855
+#> [2,] 0.8389855  0.5441538

Note that this action is equivalent to principal component analysis, and it might be wiser to use

- +
pr <- prcomp(data.frame(u,v))

to do the analysis, which yields the same principal axes

- +
pr
+#> Standard deviations (1, .., p=2):
+#> [1] 0.6953428 0.1247818
+#> 
+#> Rotation (n x k) = (2 x 2):
+#>         PC1        PC2
+#> u 0.5441538 -0.8389855
+#> v 0.8389855  0.5441538

but has the advantage of having associated methods, including for plotting.

7.6 Solution 6

Exercise 6. Plot a time series of pressure variations, to reveal tidal height.

Solution.

- -

+
time <- adp[["time"]]
+pressure <- adp[["pressure"]]
+oce.plot.ts(time, pressure)
+

As a matter of interest, a tidal analysis may be done with

- +
m <- tidem(as.sealevel(pressure, time))
+#> Note: the tidal record is too short to fit for constituents:  SA SSA MSM MM MSF MF ALP1 2Q1 SIG1 Q1 RHO1 O1 TAU1 BET1 NO1 CHI1 PI1 P1 S1 PSI1 PHI1 THE1 J1 SO1 OO1 UPS1 OQ2 EPS2 2N2 MU2 N2 NU2 GAM2 H1 H2 MKS2 LDA2 L2 T2 S2 R2 K2 MSN2 ETA2 MO3 M3 SO3 MK3 SK3 MN4 M4 SN4 MS4 MK4 S4 SK4 2SK5 2MN6 M6 2MS6 2MK6 2SM6 MSK6 M8

where the list of undetermined tidal constituents is long, because this tidal record is so short (see ?tidem). The fitted constituents are as follows

- +
summary(m)
+#> tidem summary
+#> -------------
+#> 
+#> Call:
+#> tidem(t = as.sealevel(pressure, time))
+#> RMS misfit to data:  0.06615227 
+#> 
+#> Fitted Model:
+#>         Freq Amplitude Phase       p    
+#> Z0    0.0000   39.0465   0.0 < 2e-16 ***
+#> K1    0.0418    0.0689 334.9   0.089 .  
+#> M2    0.0805    1.1730 203.8 4.6e-16 ***
+#> 2MK5  0.2028    0.0312 112.8   0.417    
+#> 3MK7  0.2833    0.0265  87.8   0.618    
+#> ---
+#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
+#> * Processing Log
+#>     - 2018-06-07 09:13:55.971 UTC: `create 'tidem' object`
+#>     - 2018-06-07 09:13:55.971 UTC: `tidem(t = as.sealevel(pressure, time))`

(Note that it fitted for M2, but not S2, because the Rayleigh criterion prevents inferring both, and the conventional tidal analysis favours M2 over S2 if the time series is too short to fit for more than one semi-diurnal constituent.)

A useful step might be to draw a tidal fit along with the data. It makes sense to construct a finer time grid, to get a smooth curve, as in the following (see ?predict.tidem).

- -

+
oce.plot.ts(time, pressure, type="p", col="blue")
+timePredict <- seq(min(time), max(time), length.out=200)
+pressurePredict <- predict(m, timePredict)
+lines(timePredict, pressurePredict, col="red")
+

diff --git a/inst/doc/flags.html b/inst/doc/flags.html index a83d2278d2..fc1fad673e 100644 --- a/inst/doc/flags.html +++ b/inst/doc/flags.html @@ -12,7 +12,7 @@ - + Using data-quality flags @@ -20,256 +20,46 @@ - + @@ -280,7 +70,7 @@

Using data-quality flags

Dan Kelley

-

2018-06-03

+

2018-06-07

@@ -307,7 +97,7 @@

2018-06-03

Abstract. This vignette explains the basics of working with data-quality flags in oce. As in the vignette about working with adp data, the material is organized to reflect typical workflow. See the main oce vignette for general information about oce. The sample code provided here relies on the oce package, loaded as follows.

- +
library(oce)

1 Overview

@@ -351,45 +141,45 @@

2 Sample working procedure

2.1 CTD range checks

(This section is an expansion of Example 1 in the setFlags,ctd-method documentation.)

The oce package provides a dataset

- +
library(oce)
+data(ctdRaw)

that contains some clearly anomolous values that are revealed clearly in a summary plot (Figure 1):

- +
plot(ctdRaw)
-Figure 1. Summary plot for raw CTD profile. +Figure 1. Summary plot for raw CTD profile.

Figure 1. Summary plot for raw CTD profile.

As Figure 1 shows, ctdTrim removes most of the anomolous data by examining the variation of the pressure signal over time, it may also be of interest to see how well simple range checks can perform in cleaning up the data. Salinity certainly cannot be negative, but in an oceanographic setting it is common to relax that criterion somewhat, perhaps insisting that Absolute Salinity \(S_A\) exceed \(25\)g/kg. This value might work in other situations as well, and the same could be said of an upper limit of \(40\)g/kg. Similarly, it might make sense to bound temperature between, say, \(-2^\circ\)C and \(40^\circ\)C for application throughout much of the world ocean.

These criteria can be supplied to setFlags in various ways, but the simplest is to create logical vectors, e.g.

- +
badS <- with(ctdRaw[["data"]], salinity < 25 | 40 < salinity)
+badT <- with(ctdRaw[["data"]], temperature < -2 | 40 < temperature)

In the above, with has been used to avoid inserting salinity and temperature in the namespace, but it would also be common to use e.g.

- +
S <- ctdRaw[["salinity"]]
+T <- ctdRaw[["temperature"]]
+bad <- (S < 25 | 40 < S) || (T < -2 | 40 < T)

Since the goal here is to illustrate setting multiple flags, the badS and badT values will be used. The first step is to copy the original data, so that the flag operations will not alter ctdRaw:

- +
qc <- ctdRaw

Work flow is best documented if a flag scheme is established, and the “WHP CTD exchange” scheme is a reasonable choice; using

- +
qc <- initializeFlagScheme(qc, "WHP CTD")

stores a note on the scheme in the metadata slot of qc. Importantly, however, it does not store any flag values. The next step is to initialize flag values. For example, to set up flag storage for salinity and temperature, use e.g.

- +
qc <- initializeFlags(qc, "salinity", 2)
+qc <- initializeFlags(qc, "temperature", 2)

to create storage, and initialize all entries to the "acceptable" value (which, for this flag scheme, is the number 2). Once storage has been initialized for a given variable, further calls to initializeFlags will have no effect, apart from a warning.

Note that it is entirely possible to use initializeFlags with a numerical value for its third argument. One advantage of using initializeFlagScheme is that it clarifies code, but the bigger advantage is that it embeds the scheme within the object, so that a second analyst could examine it later and be clear on the meanings of the numerical codes. This is very important, because there is not much agreement within the oceanographic community on which flag scheme to use, for many data types (Argo being an exception).

At this stage, individidual data can be flagged with setFlags. This function can be called any number of times. Continuing along with our example, we may mark bad salinities with

- +
qc <- setFlags(qc, "salinity", badS, value="bad")

We can see that the flag got inserted by using summary(qc), but for brevity here another method is:

- +
names(qc[["flags"]])
+#> [1] "salinity"    "temperature"

Now, temperature flags may be inserted with

- +
qc <- setFlags(qc, "temperature", badT, value="bad")

Readers ought to use summary(qc) to get more details of how flags were handled, and how many bad salinities and temperatures were flagged.

If qc is plotted with plot(qc), the results will match those of plot(ctdRaw). This is because setting flags has no effect on plots, because it alters flags but not data. One more step is required to test whether this procedure has cleaned up the data significantly: we must “handle” the flags, using

- +
qch <- handleFlags(qc, flags=list(c(1, 3:9)))

Comparing Figure 1 with the summary plot for qch (Figure 2), constructed with

- +
plot(qch)
-Figure 2. Summary plot for range-checked CTD profile. +Figure 2. Summary plot for range-checked CTD profile.

Figure 2. Summary plot for range-checked CTD profile.

shows significant improvements. The downcast and the upcast can be seen quite clearly now, although there appears to be an issue of low salinity at the turnaround point. Setting a flag for pressure increase with time will isolate the downcast somewhat, although some smoothing will be required. Another issue related to the path of the instrument is that it may have been held below the surface for a while to equilibrate. Again, a flag could be set up to remove such data. However, it ought to be noted that ctdTrim can be used to address issues relating to instrument movement.

@@ -399,59 +189,59 @@

2.2 Interactive editing of CTD pr

(This section is an expansion of Example 2 in the setFlags,ctd-method documentation.)

The ctd dataset provided by the oce package is similar to ctdRaw, except that only downcast data are provided. Even so, there are still some points that might be considered suspicious. A common way to find such points is to plot TS diagrams, looking for decreases in density with depth.

Running the following code in an interactive session will demonstrate a simple way to use a TS diagram to identify suspicious data. Pasting this into an R console will show a plot with lines between measurements made at successive depths. Clicking on any point will flag it, and so the point will then disappear on the plot. Clicking to the right of the plot frame will exit the procedure, after which qc is a ctd object with flags as set, and data set to NA where these flags indicate bad data.

- +
options(eos="gsw")
+data(ctd)
+qc <- ctd
+qc <- initializeFlagScheme(qc, "WHP CTD")
+qc <- initializeFlags(qc, "salinity", 2)
+Sspan <- diff(range(qc[["SA"]]))
+Tspan <- diff(range(qc[["CT"]]))
+n <- length(qc[["SA"]])
+par(mfrow=c(1, 1))
+plotTS(qc, type="o")
+message("Click on bad points; quit by clicking to right of plot")
+for (i in seq_len(n)) {
+    xy <- locator(1)
+    if (xy$x > par("usr")[2])
+        break
+    i <- which.min(abs(qc[["SA"]] - xy$x)/Sspan + abs(qc[["CT"]] - xy$y)/Tspan)
+    qc <- setFlags(qc, "salinity", i=i, value="bad")
+    qc <- handleFlags(qc)
+    plotTS(qc, type="o")
+}

It would be a simple matter to extend this simple example to a shiny application that displays other data. For example, there could be panels for profiles, as well as the TS plot. The system could track clicks in any of the panels, taking appropriate actions. It would be sensible to have a staged procedure, in which clicking (or brushing) in one panel would cause a replot of all panels, with the selected data indicated in some way, so that the analyst could then choose whether to go to the next stage, of clicking a button to indicate bad data. Another button might be provided to undo such operations, or to show the original data for comparison. The point is that a wide variety of flag operations are handled very easily in R, with oce.

2.3 Hydrographic sections

The flag-handling functions work for a variety of oce objects. For example, sections, which are built up from a sequence of ctd profiles, are handled easily with functions of the same names as for the ctd case.

As a simple example, the following shows how to clean up the “A03” Atlantic section that is provided with oce.

- -

+
data(section)
+s <- handleFlags(section, flags=list(c(1, 3:9)))
+par(mfrow=c(2, 1))
+plotTS(section)
+plotTS(s)
+

Note, in the above, that the archiving agency had evidently flagged not just wild data (e.g. the salinity near \(26\)g/kg) but also data that were anomalous in more subtle ways (e.g. cleaning up several points that stood out from the data cloud, below \(10^\circ\)C). In fact, the flags in this data set, as in most archived hydrographic data sets, are the result of a multifaceted inspection scheme that is more demanding and useful than simple range checks.

3 Acoustic-Doppler profiler data

By now, the reader should be able to understand the following use of data-quality flags for the adp dataset. Note that there is only a small difference (at 8h and 20h) between the (gray) bad data as identified by the instrument and those identified by the scheme illustrated here; look near 8h and 20h. However, altering the values for G and V4 will reveal that the schemes do differ.

- -

+
data(adp)
+v <- adp[["v"]]
+i2 <- array(FALSE, dim=dim(v)) # construct array to match 'v'
+g <- adp[["g", "numeric"]]
+G <- 25                  # for percent-good field, named 'g'
+V4 <- 0.45               # for error velocity field, in beam 4
+for (k in 1:3)
+    i2[,,k] <- ((g[,,k]+g[,,4]) < G) | (v[,,4] > V4)
+adpQC2 <- initializeFlags(adp, "v", 2)
+adpQC2 <- setFlags(adpQC2, "v", i2, 3)
+adpClean2 <- handleFlags(adpQC2, flags=list(3), actions=list("NA"))
+par(mfrow=c(2,1))
+plot(adp, which="u1")                  # original
+plot(adpClean2, which="u1")            # altered
+

Note, in the above, that initializeFlagScheme has not been used, because the author is unaware of any common notation for such data. The values used for G and V4 were provided by colleagues at the Bedford Institute of Oceanography.

diff --git a/inst/doc/oce.html b/inst/doc/oce.html index af0e1c92c8..9dc1c22941 100644 --- a/inst/doc/oce.html +++ b/inst/doc/oce.html @@ -12,7 +12,7 @@ - + Oce @@ -20,256 +20,46 @@ - + @@ -280,7 +70,7 @@

Oce

Dan Kelley and Clark Richards

-

2018-06-03

+

2018-06-07

@@ -346,11 +136,11 @@

2 Architecture of oce objects

3 Accessing data in oce objects

For detailed analysis, users may want to access data within oce objects. While it is possible to descend through the object using the “slot” and “list” notation (e.g. d@data$salinity yields the salinity within an oce object named d), this approach is not recommended. It is better to use the [[ notation, which derives from the generic “Extract” method for accessing parts of an object. A general introduction is provided by

- +
?`[[,oce-method`

and detaile are provided for individual object classes with e.g.

- +
?`[[,ctd-method`

for the ctd class. The notation is very simple. For example, suppose that d is an object that stores salinity in either its data slot or metadata slot. Then,

- +
S <- d[['salinity']]

will extract the salinity.

The [[ method first looks in the metadata slot, but if the item is not found there, it proceeds to look in the data slot. This two-step scheme is helpful because it frees the user from having to know where data are stored, e.g. in a ctd object it might be stored in the metadata (for a conventional CTD cast) or in the data slot (for the slantwise sampling of a glider).

In addition to accessing data within an object, the [[ notation also permits the extraction of derived information. There are two reasons for this.

@@ -359,45 +149,45 @@

3 Accessing data in oce objectsIt provides a uniformity of notation that can be helpful to users. In oce, objects derived from data files tend to hold just the information in those files, not derived information. For example, For example, CTD datasets often provide in-situ temperature but not potential temperature (since the latter is not measured). The in-situ temperature is found with e.g. d[["temperature"]], and so it seems natural to write d[["theta"]] to get the potential temperature. If this is done, then the ctd version of [[`` first looks in the data slot, and returnsthetaif it is found there, as may sometimes be the case (depending on the choice of the person who created the dataset). However, if it is not found, then[[callsswTheta()` to calculate the value, and returns the result.

Finally, it is also possible to extract the entirety of either the metadata or data slot, e.g.

- +
data <- d[['data']]

yields the full data slot, which is a list with elements that can be accessed in the conventional way, e.g. for a ctd object,

- +
data$temperature

retrieves the temperature. For obvious reasons, the method of derived-quantity access (e.g. the theta of the example above) will not work.

4 Modifying data in oce objects

There are several schemes for modifying the data within oce objects. By analogy with the [[ notation of the previous section, e.g. the following

- +
data(ctd)
+ctd[["temperatureAboveFreezing"]] <- ctd[["temperature"]] - swTFreeze(ctd)

will store the excess over freezing temperature into the ctd object.

Further information on this notation is provided by

- +
?"[[<-,oce-method"

The above works only within the data slot. To store within the metadata slot, consider using e.g.

- +
ctd[["metadata"]]$scientist <- "Dalhousie Oceanography 4120/5120 Class of 2003"

sets the “scientist”.

For archival work, it is important to store reasons for changing things in an object. Two functions are provided for this purpose: oceSetData and oceSetMetadata. For example, a better way to change the scientist might be to write

- +
ctd <- oceSetMetadata(ctd, name="scientist",
+                      value="Dalhousie Oceanography 4120/5120 Class of 2003",
+                      note="give credit where it's due")

and a better way to store temperatureAboveFreezing would be

- +
ctd <- oceSetData(ctd, name="temperatureAboveFreezing",
+                      value=ctd[["temperature"]] - swTFreeze(ctd), 
+                      unit=list(unit=expression(degree*C), scale="ITS-90"),
+                      originalName="-",
+                      note="add temperatureAboveFreezing, for ice-related calculations")

which illustrates that this notation, as opposed to the [[<- notation, permits the specification of a unit and an originalName, both of which, together with the note, are displayed by summary(ctd).

5 Oce functions

The uniformity of the various oce objects helps users build skill in examining and modifying objects. Fine-scale control is provided throughout oce, but the best way to learn is to start with simplest tools and their defaults. For example, the following will read a CTD file named "station1.cnv", summarize the contents, and plot an overview of the data, with profiles, a TS diagram, and a map (Figure 2).

- +
library(oce)
+d <- read.oce("station1.cnv")
+summary(d)
+plot(d)

The reader should stop now and try this on a file of their own. The pattern will work with a fairly wide variety of file types, because read.oce() examines the file name and contents to try to discover what it is. For an example, if read.oce() is given the name of a file created by an Acoustic Doppler Profiler, it will return an object inheriting from class "adp", so the summary() and plot() calls will be tailored to that type, e.g. the graph will show images of time-distance variation in each of the measured velocity components.

Notes on oce function names.

    -

  1. As just illustrated, the general function to read a dataset ends in .oce, and the name is a signal that the returned object is of class oce. Depending on the file contents, d will also have an additional class, e.g. if the file contains CTD data, then the object would inherit from two classes, oce and ctd, with the second being used to tailor the graphics by passing control to the ctd variant of the generic plot function (use help("plot,ctd-method") to learn more).

    2. +

  2. As just illustrated, the general function to read a dataset ends in .oce, and the name is a signal that the returned object is of class oce. Depending on the file contents, d will also have an additional class, e.g. if the file contains CTD data, then the object would inherit from two classes, oce and ctd, with the second being used to tailor the graphics by passing control to the ctd variant of the generic plot function (use help("plot,ctd-method") to learn more).

  3. Generally, oce functions employ a “camel case” naming convention, in which a function that is described by several words is named by stringing the words together, capitalizing the second and subsequent words. For example, ctdAddColumn() takes a ctd object, and returns a new ctd object that matches the original except for the added column (and an added entry in the processingLog to indicate the addition).

  4. Function names begin with oce in cases where a more natural name would be in conflict with a function in the base R system or a package commonly used by Oceanographers. For example, oceContour() is a function that provides an alternative to contour() in the graphics package.

@@ -439,40 +229,40 @@

7 CTD data

7.1 Example with pre-trimmed data

Many of the object types supported by oce come with built-in data. For an example, data(ctd) yields a CTD profile that has been trimmed to just the downcast portion of the sampling. (See the next section to learn how to do this trimming.) A summary and plot (Figure 2) are created as follows.

- +
library(oce)
+data(ctd)
+summary(ctd)
+#> CTD Summary
+#> -----------
+#> 
+#> * Instrument:          SBE 25 
+#> * Temp. serial no.:    1140
+#> * Cond. serial no.:    832
+#> * File:                /Users/kelley/git/oce/create_data/ctd/ctd.cnv
+#> * Original file:       c:\seasoft3\basin\bed0302.hex
+#> * Start time:          2003-10-15 11:38:38
+#> * Sample interval:     1 s
+#> * Cruise:              Halifax Harbour
+#> * Vessel:              Divcom3
+#> * Station:             Stn 2
+#> * Location:            44.684N  63.644W 
+#> * Data
+#> 
+#>                               Min.   Mean   Max.   Dim. OriginalName
+#>     scan                      130    220    310    181  scan        
+#>     pressure [dbar]           1.48   22.885 44.141 181  pr          
+#>     depth [m]                 1.468  22.698 43.778 181  depS        
+#>     temperature [°C, IPTS-68] 2.919  7.5063 14.237 181  t068        
+#>     salinity [PSS-78]         29.916 31.219 31.498 181  sal00       
+#>     flag                      0      0      0      181  flag        
+#> 
+#> * Processing Log
+#>     - 2017-09-09 12:00:27.237 UTC: `create 'ctd' object`
+#>     - 2017-09-09 12:00:27.648 UTC: `read.ctd.sbe(file = file, processingLog = processingLog)`
+#>     - 2017-09-09 12:00:27.716 UTC: `oce.edit(x = ctd, item = "startTime", value = ctd[["systemUploadTime"]],     reason = "startTime in file has Y2K error", person = "Dan Kelley")`
+plot(ctd)
-Figure 2. An overview of a ctd dataset. +Figure 2. An overview of a ctd dataset.

Figure 2. An overview of a ctd dataset.

Accessing the data within this ctd object can be done directly, e.g. ctd@data$pressure holds the pressure record, but it is usually better to use an accessor function that is provided with oce. This function is named [[, and it takes a character string as an argument, e.g. ctd[["pressure"]] yields the pressure column. The accessor notation is preferable to direct access because it is simpler for the user. For example, several oce objects store the data in single-byte or two-byte chunks, to match the raw format used by the instruments, and the accessor function takes care of translating these values to what are sometimes called “science” units.

@@ -481,41 +271,41 @@

7.1 Example with pre-trimmed data

7.2 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. Spurious data from these phases must be trimmed as a first step in processing. For example, consider the following code.

- +
data(ctdRaw)
+plotScan(ctdRaw)
-Figure 3. Scanwise plot of the ctdRaw sample data set. Note the spike at the start, the equilibration phase before the downcast, and the spurious freshening signal near the start of the upcast. +Figure 3. Scanwise plot of the ctdRaw sample data set. Note the spike at the start, the equilibration phase before the downcast, and the spurious freshening signal near the start of the upcast.

Figure 3. Scanwise plot of the ctdRaw sample data set. Note the spike at the start, the equilibration phase before the downcast, and the spurious freshening signal near the start of the upcast.

This produces a two-panel plot (Figure 3) of the data as a time-series, revealing not just the desired downcast, but also an earlier equilibration phase and a later upcast. The x-axis in Figure 3 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

- +
plotScan(ctdTrim(ctdRaw, "range",
+                 parameters=list(item="scan", from=140, to=250)))
+plotScan(ctdTrim(ctdRaw, "range",
+                 parameters=list(item="scan", from=150, to=250)))

This method of making decisions based on plotted information is probably the most robust method of trimming data. However, for quick work, users may be satisfied with the results of automatic downcast detection, e.g.

- +
ctdTrimmed <- ctdTrim(ctdRaw)

It should be noted that ctdTrim() inserts entries into the object’s log file, so that the details of how the trimming was done are recorded together with the data.

Once the profile has been trimmed, one may wish to use ctd.decimate() to smooth the data and interpolate the smoothed results to uniformly-spaced pressure values.

Taking these things together, a quick visual examination of a CTD file takes just one line of code:

- +
plot(ctdDecimate(ctdTrim(read.ctd("stn123.cnv"))))

7.3 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 CTD,20060609WHPOSIODAM (or BOTTLE,...). Scanning the item to the left of the comma is not difficult (although there are variants to the two shown, e.g. CTDO sometimes occurs). The part to the right of the comma is more difficult. The first part is a date (yyyymmdd) so that is 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 WHOI occurs in the header line), the approach taken with oce is to ignore such issues relating to quirky headers, on the assumption that users can scan human-written headers with high skill.

Quite commonly, CTD profiles taken during a cruise are collected together in a sequence of files in a given directory. For a real-world example, one might visit the website mentioned in the code provided below, download and expand the zip file, enter the directory thus formed, and run the code to get an overall TS plot for all the CTD stations of this cruise. (Caution: the link seems to change from time to time.)

- +
library(oce)
+# http://cchdo.ucsd.edu/data/7971/ar18_58JH19941029_ct1.zip
+# setwd("~/Downloads/ar18_58JH19941029_ct1")
+files <- system("ls *.csv", intern=TRUE)
+for (i in 1:length(files)) {
+    x <- read.ctd(files[i])
+    if (i == 1) {
+        plotTS(x, Slim=c(31, 35.5), Tlim=c(-2, 10), type='o')
+    } else {
+        points(x[["salinity"]], x[["potential temperature"]])
+        lines(x[["salinity"]], x[["potential temperature"]])
+    }
+}

The [[ notation is explained in Section 3, but this example conveys the gist, that it permits accessing data, or derived data, from within an object.

In the above, lines connect the points within a given profile. This can be a useful method for a quick scan looking for outliers. Another is to colour-code the profiles, although this gets confusing with large datasets, in which case the method of the following exercise might be useful.

Exercise 3. (advanced) Make a multi-file plot summarizing the TS relationship, with each file showing the overall relationship in gray and the individual profile in black.

@@ -524,17 +314,17 @@

7.3 Example with WOCE archive dat

8 Section plots

The commands

- +
data(section)
+plot(section, which=c(1, 2, 3, 99))

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, the section variant of the generic plot function can also create individual plots of individual properties (use help("plot,section-method") to learn more).

Some section datasets are supplied in a pre-gridded format, but it is also common to have different pressure levels at each station. For such cases, sectionGrid() may be used, e.g. the following produces Figure 4. The ship was travelling westward from the Mediterranean towards North America, taking 124 stations in total; the stationId value selects the last few stations of the section, during which the ship headed toward the northwest, crossing isobaths (and perhaps, the Gulf Stream) at nearly right angles.

- +
library(oce)
+data(section)
+GS <- subset(section, 102 <= stationId & stationId <= 124)
+GSg <- sectionGrid(GS, p=seq(0, 1600, 25))
+plot(GSg, which=c(1,99), map.xlim=c(-85,-(64+13/60)))
-Figure 4. Portion of the CTD section designated A03, showing the Gulf Sream region. The square on the cruise track corresponds to zero distance on the section. +Figure 4. Portion of the CTD section designated A03, showing the Gulf Sream region. The square on the cruise track corresponds to zero distance on the section.

Figure 4. Portion of the CTD section designated A03, showing the Gulf Sream region. The square on the cruise track corresponds to zero distance on the section.

Exercise 4. Draw a \(T\)-\(S\) diagram for the section data, using black symbols to the east of 30W and gray symbols to the west, thus highlighting Mediterranean-derived waters. Use handleFlags() (see using data-quality flags) to discard questionable data, and use the accessor function [[.

@@ -545,16 +335,16 @@

9 Maps

9.1 About projections

There are several ways to handle map projections in R and other systems. Oce uses the rgdal package to calculate the details of map projections. A wide group of functions is provided for plotting maps. The first step is to call mapPlot() to construct the map, after which points can be added with mapPoints(), lines with mapLines(), etc. For example, the following plots a world coastline in the Winkel Tripel projection, with the cruise track for the section dataset in red.

- +
library(oce)
+data(coastlineWorld)
+par(mar=rep(0, 4))
+mapPlot(coastlineWorld, projection="+proj=robin", col="lightgray")
+data(section)
+lon <- section[["longitude", "byStation"]]
+lat <- section[["latitude", "byStation"]]
+mapLines(lon, lat, col='red', cex=0.5)
-Figure 5. World in Robinson projection, with the section profile sites indicated. +Figure 5. World in Robinson projection, with the section profile sites indicated.

Figure 5. World in Robinson projection, with the section profile sites indicated.

There are far too many projections to illustrate here. See http://dankelley.github.io/r/2015/04/03/oce-proj.html for a blog item that provides examples of the available projections. Note that some have a problem with spurious horizontal lines. This can result from coastline segments that cross the edge of the plotting area, and getting rid of this is tricky enough to be the heart of the longest-lived bug in the oce issue list, i.e. https://github.com/dankelley/oce/issues/388. In some instances the function coastlineCut() can help, but it is provisional and subject to change.

@@ -563,47 +353,47 @@

9.1 About projections

9.2 Popular world-view projections

A few common projections are worth highlighting. In addition to the Robinson projection shown in Figure 5, it is also common to employ Mollweide (projection="+proj=moll"), Eckert IV (projection="+eck4") or Winkel Tripel (projection="+proj=wintri") for world-wide views. An attractive feature of the Winkel Tripel projection is familiarity, since it is used by the National Geographical Society, but its lack of a formal inverse can cause problems for plotting in oce unless the operating system has a recent version of the gdal library that is used for calculating projections.

Before moving on to other projections, it may be helpful to show how to plot gridded data with projection, illustrated in Figure 6 with bathymetry near eastern Canada.

- +
par(mar=c(1.5, 1, 1.5, 1))
+data(topoWorld)
+topo <- decimate(topoWorld, 2) # coarsen grid: 4X faster plot
+lon <- topo[["longitude"]]
+lat <- topo[["latitude"]]
+z <- topo[["z"]]
+cm <- colormap(name="gmt_globe")
+drawPalette(colormap=cm)
+mapPlot(coastlineWorld, projection="+proj=moll", grid=FALSE, col="lightgray")
+mapImage(lon, lat, z, colormap=cm)
-Figure 6. World bathymetry in Molleweide projection. +Figure 6. World bathymetry in Molleweide projection.

Figure 6. World bathymetry in Molleweide projection.

@@ -613,11 +403,11 @@

10 Sea-level data

10.1 Time-domain analysis

The following code graphs a built-in dataset of sea-level time series (Figure 9). The top panel 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 log spectrum, with a few tidal constituents indicated. The section variant of the generic plot function provides other possibilities for plotting, including a cumulative spectrum that can be quite informative (use help("plot,sealevel-method") to learn more).

- +
library(oce)
+data(sealevel)
+plot(sealevel)
-Figure 9. Sea-level timeseries measured in 2003 in Halifax Harbour. +Figure 9. Sea-level timeseries measured in 2003 in Halifax Harbour.

Figure 9. Sea-level timeseries measured in 2003 in Halifax Harbour.

Exercise 6. Illustrate Halifax sea-level variations during Hurricane Juan, near the end of September, 2003.

@@ -631,12 +421,12 @@

10.2 Tidal-harmonic analysis

11 Acoustic Doppler Profiler data

The following commands produce Figure 10, showing one velocity component of currents measured in the St Lawrence Estuary Internal Wave Experiment. This plot type is just one of many provided by the adp variant of the generic plot function (see ?"plot,adp-method"). See the adp vignette for much more on acoustic Doppler profiler data.

- +
library(oce)
+data(adp)
+plot(adp, which=1)
+lines(adp[['time']], adp[['pressure']], lwd=2)
-Figure 10. Measurements made with a bottom-mounted ADP in the St Lawrence Estuary. The line near the surface indicates pressure measured by the ADP. +Figure 10. Measurements made with a bottom-mounted ADP in the St Lawrence Estuary. The line near the surface indicates pressure measured by the ADP.

Figure 10. Measurements made with a bottom-mounted ADP in the St Lawrence Estuary. The line near the surface indicates pressure measured by the ADP.

@@ -647,10 +437,10 @@

12 Handling data-quality flags

13 Working in non-English languages

Many of the oce plotting functions produce axis labels that can be displayed in languages other than English. At the time of writing, French, German, Spanish, and Mandarin are supported in at least a rudimentary way. Setting the language can be done at the general system level, or within R, as indicated below (results not shown).

- +
library(oce)
+Sys.setenv(LANGUAGE="fr")
+data(ctd)
+plot(ctd)

Most of the translated items were found by online dictionaries, and so they may be incorrect for oceanographic usage. Readers can help out in the translation effort, if they have knowledge of how nautical words such as Pitch and Roll and technical terms such as Absolute Salinity and Potential Temperature should be written in various languages.

@@ -659,62 +449,62 @@

14 Extending oce

14.1 Initializing oce objects

Oce handles a wide variety of data types, creating sub-classes the inherit from the base class that is named oce. However, even if a dataset does not fit within these sub-classes, it is still possible to use some of the features of oce, such as (FILL IN … accessors, plotting, summaries ??).

For example, suppose that there is a dataset that contains two variables, mm and ss, that have been measured at times t. Suppose that the data are recorded at a fixed location. To illustrate how to deal with such data, we first create some fake data:

- +
lon <- -60
+lat <- -30
+t <- seq(as.POSIXct("2017-08-30"), as.POSIXct("2017-09-01"), by="15 min")
+tsec <- as.numeric(t)
+u <- sin(tsec * 2 * pi / (3600 * 12.4206))
+v <- 0.5 * cos(tsec * 2 * pi / (3600 * 12.4206))

These may be inserted into an oce object as follows. First, create the object.

- +
o <- new("oce")

This contains slots without data:

- +
str(o)
+#> Formal class 'oce' [package "oce"] with 3 slots
+#>   ..@ metadata     :List of 2
+#>   .. ..$ units: list()
+#>   .. ..$ flags: list()
+#>   ..@ data         : list()
+#>   ..@ processingLog:List of 2
+#>   .. ..$ time : POSIXct[1:1], format: "2018-06-07 08:49:26.787"
+#>   .. ..$ value: chr "Create oce object"

As noted previously, it’s possible to insert t, u, and v directly into the data slot, but it’s better to use oceSetData to do that, because it can set up units, and it leaves a trail in the processing log.

- +
o <- oceSetData(o, "t", t, list(unit=expression(s), scale=""))
+o <- oceSetData(o, "u", u, list(unit=expression(m/s), scale=""))
+o <- oceSetData(o, "v", v, list(unit=expression(m/s), scale=""))

Astute readers will realize that these data are akin to tidal velocities, and so might agree with the author that lon and lat belong in the metadata slot, not the data slot.

- +
o <- oceSetMetadata(o, "longitude", lon)
+o <- oceSetMetadata(o, "latitude", lat)

At this stage, o is a full-fledged oce object. It may be summarized with

- +
summary(o)
+#> * Data
+#> 
+#>             Min.       Mean       Max.       Dim. OriginalName
+#>     t [s]   2017-08-30 2017-08-31 2017-09-01 193  -           
+#>     u [m/s] -0.99982   -0.027415  1          193  -           
+#>     v [m/s] -0.49996   -0.0048211 0.49998    193  -           
+#> 
+#> * Processing Log
+#>     - 2018-06-07 08:49:26.787 UTC: `Create oce object`
+#>     - 2018-06-07 09:14:23.914 UTC: `oceSetData(object = o, name = "t", value = t, unit = list(unit = expression(s),     scale = ""))`
+#>     - 2018-06-07 09:14:23.917 UTC: `oceSetData(object = o, name = "u", value = u, unit = list(unit = expression(m/s),     scale = ""))`
+#>     - 2018-06-07 09:14:23.919 UTC: `oceSetData(object = o, name = "v", value = v, unit = list(unit = expression(m/s),     scale = ""))`
+#>     - 2018-06-07 09:14:23.938 UTC: `oceSetMetadata(object = o, name = "longitude", value = lon)`
+#>     - 2018-06-07 09:14:23.941 UTC: `oceSetMetadata(object = o, name = "latitude", value = lat)`

and individual data components may be extracted with the [[ notation

- +
o[["latitude"]]                        # from metadata
+#> [1] -30
+head(o[["t"]])                         # from data
+#> [1] "2017-08-30 00:00:00 ADT" "2017-08-30 00:15:00 ADT"
+#> [3] "2017-08-30 00:30:00 ADT" "2017-08-30 00:45:00 ADT"
+#> [5] "2017-08-30 01:00:00 ADT" "2017-08-30 01:15:00 ADT"
+head(o[["u"]])
+#> [1] 0.9961884 0.9772306 0.9426639 0.8930403 0.8291526 0.7520211
+head(o[["v"]])
+#> [1] -0.04361357 -0.10608997 -0.16687184 -0.22498833 -0.27951117 -0.32956949

A crude plot can be done with

- +
plot(o)
-Figure 11. Crude plot of constructed oce object. +Figure 11. Crude plot of constructed oce object.

Figure 11. Crude plot of constructed oce object.

Note that the plot uses the pairs function to show how all elements of the data slot vary with each other. To customize the plot, a new object class should be created, as in the next section.

@@ -722,77 +512,77 @@

14.1 Initializing oce objects

14.2 Creating new oce object classes

Let’s call this new object class uv. With this, inventing a new object sub-class is as simple as invoking

- +
setClass("uv", contains="oce")

For most objects, it makes sense to define three specialized methods: one to create new objects (called "initialize"), one to plot them ("plot"), and one to summarize them ("summary").

The initialization can be handled with

- +
setMethod(f="initialize",
+          signature="uv",
+          definition=function(.Object, time, u, v, longitude, latitude) {
+              if (missing(time)) stop("must provide 'time'")
+              ## repeat the above also for u, v, longitude, and latitude
+              .Object@metadata$units$u <- list(unit=expression(m/s), scale="")
+              .Object@metadata$units$v <- list(unit=expression(m/s), scale="")
+              .Object@metadata$longitude <- longitude
+              .Object@metadata$latitude <- latitude
+              .Object@data$time <- time
+              .Object@data$u <- u
+              .Object@data$v <- v
+              .Object@processingLog$time <- as.POSIXct(Sys.time())
+              .Object@processingLog$value <- "create 'uv' object"
+              return(.Object)
+          })

This is fairly arcane, but simply copying and pasting this will be enough to get the reader started. Let’s try it:

- +
oo <- new("uv", t, u, v, lon, lat)

Next, customize summary for this new class:

- +
setMethod(f="summary",
+          signature="uv",
+          definition=function(object, ...) {
+              cat("uv Summary\n-----------\n\n", ...)
+              cat(paste("* Location:           ",
+                        sprintf("%.5f N", object@metadata$latitude), ", ",
+                        sprintf("%.5f E", object@metadata$longitude), "\n", sep=''))
+              callNextMethod()
+          })

The callNextMethod will summarize the entries in the data slot adequately for this case. Let’s try it:

- +
summary(oo)
+#> uv Summary
+#> -----------
+#> 
+#> * Location:           -30.00000 N, -60.00000 E
+#> * Time ranges from 2017-08-30 to 2017-09-01 with 193 samples and mean increment 15 min
+#> * Data
+#> 
+#>             Min.     Mean       Max.    Dim. OriginalName
+#>     u [m/s] -0.99982 -0.027415  1       193  -           
+#>     v [m/s] -0.49996 -0.0048211 0.49998 193  -           
+#> 
+#> * Processing Log
+#>     - 2018-06-07 09:14:24.289 UTC: `create 'uv' object`
+#> NULL

Finally, let’s make a custom method for plot. This is usually the most complicated part of handling any new object class, because there are usually many ways to look at data. (Some oce objects have over a dozen specialized plot styles.) Let’s suppose there are only two desired plot styles: one showing timeseries of u and v in a two panel plot, and the other showing how u and v covary.

- +
setMethod(f="plot",
+          signature=signature("uv"),
+          definition=function(x, which=1, ...)
+          {
+              if (which == 1) {
+                  oce.plot.ts(x[["time"]], x[["u"]], ylab="u [m/s]", ...)
+              } else if (which == 2) {
+                  oce.plot.ts(x[["time"]], x[["v"]], ylab="v [m/s]", ...)
+              } else if (which == 3) {
+                  plot(x[["u"]], x[["v"]], xlab="u [m/s]", ylab="v [m/s]", ...)
+              }
+          })

An now, we can try it out

- +
par(mfrow=c(2,1))
+plot(oo, which=1)
+plot(oo, which=2)
-Figure 12. Time-series u(t) plot of newly-constructed uv object. +Figure 12. Time-series u(t) plot of newly-constructed uv object.

Figure 12. Time-series u(t) plot of newly-constructed uv object.

- +
plot(oo, which=3, asp=1)
-Figure 13. Hodograph plot of newly-constructed uv object. +Figure 13. Hodograph plot of newly-constructed uv object.

Figure 13. Hodograph plot of newly-constructed uv object.

@@ -809,110 +599,110 @@

16 Development website

17 Answers to exercises

Exercise 1. Seawater properties. In the UNESCO system we may write

- +
library(oce)
+swRho(34, 10, 100, eos="unesco")
+#> [1] 1026.624
+swTheta(34, 10, 100, eos="unesco")
+#> [1] 9.988599
+swRho(34, swTheta(34, 10, 100, eos="unesco"), 0, eos="unesco")
+#> [1] 1026.173
+swRho(34, swTheta(34, 10, 100, 200, eos="unesco"), 200, eos="unesco")
+#> [1] 1027.074
+plotTS(as.ctd(c(30,40),c(-2,20),rep(0,2)), eos="unesco", grid=TRUE, col="white")

and in the Gibbs SeaWater system, we use eos="gws" and must supply longitude and latitude arguments to the sw*() calls and also to the as.ctd() call.

Exercise 2. Plot a profile of \(\sigma_\theta\) and \(N^2\), for the data in the pycnocline. (Hint: use subset().)

Although one may argue as to the limits of the pycnocline, for illustration let us say it is in 5 dbar to 12 dbar range. One way to do this is

- +
library(oce)
+data(ctd)
+pycnocline <- ctdTrim(ctd, "range",
+                      parameters=list(item="pressure", from=5, to=12))
+plotProfile(pycnocline, which="density+N2")

Another is

- +
library(oce)
+data(ctd)
+pycnocline <- subset(ctd, 5 <= pressure & pressure <= 12)
+plotProfile(pycnocline, which="density+N2")

Exercise 3. Make a multi-file plot summarizing the TS relationship, with each file showing the overall relationship in gray and the individual profile in black. (This is an advanced exercise requiring some R skills.)

The code provided below creates 91 PNG files, with names ar18_01.png, ar18_02.png, etc. Loading these in a view that permits quick paging through this file list is an easy way to spot suspicious data, since each plot has the station number at the top. (Users trying this example should bear in mind that this is a fairly large dataset, so the processing will take up to a minute.)

- +
library(oce)
+# http://cchdo.ucsd.edu/data/7971/ar18_58JH19941029_ct1.zip
+# setwd("~/Downloads/ar18_58JH19941029_ct1")
+files <- system("ls *.csv", intern=TRUE)
+n <- length(files)
+ctds <- vector("list", n) # to hold the CTD objects
+station <- vector("list", n)
+for (i in 1:n) {
+    ctds[[i]] <- read.ctd(files[i])
+    station[[i]] <- ctds[[i]][["station"]]
+}
+S <- unlist(lapply(1:n, function(i) ctds[[i]][["salinity"]]))
+T <- unlist(lapply(1:n, function(i) ctds[[i]][["temperature"]]))
+p <- unlist(lapply(1:n, function(i) ctds[[i]][["pressure"]]))
+overall <- as.ctd(S, T, p)
+png("ar18_%02d.png")
+for (i in 1:n) {
+    plotTS(overall, col='gray')
+    lines(ctds[[i]][["salinity"]], ctds[[i]][["potential temperature"]])
+    mtext(station[i], side=3, line=0)
+}
+dev.off()

Exercise 4. Draw a \(T\)-\(S\) diagram for the section data, using black symbols to the east of 30W and gray symbols to the west, thus highlighting Mediterranean-derived waters. Use handleFlags() (see using data-quality flags) to discard questionable data, and use the accessor function [[.

We will use the Gibbs SeaWater system, so as.ctd() needs location information.

- +
library(oce)
+data(section)
+s <- handleFlags(section, flags=list(c(1, 3:9)))
+ctd <- as.ctd(s[["salinity"]], s[["temperature"]], s[["pressure"]],
+              longitude=s[["longitude"]], latitude=s[["latitude"]])
+col <- ifelse(s[["longitude"]] > -30, "black", "gray")
+plotTS(ctd, col=col, eos="gsw")

Exercise 5. Plot dynamic height and geostrophic velocity across the Gulf Stream. (Hint: use the swDynamicHeight function.)

(Try ?swDynamicHeight for hints on smoothing.)

- +
library(oce)
+data(section)
+GS <- subset(section, 102 <= stationId & stationId <= 124)
+dh <- swDynamicHeight(GS)
+par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0))
+plot(dh$distance, dh$height, type="l", xlab="", ylab="Dyn. Height [m]")
+grid()
+# 1e3 metres per kilometre
+latMean <- mean(GS[["latitude"]])
+f <- coriolis(latMean)
+g <- gravity(latMean)
+v <- diff(dh$height)/diff(dh$distance) * g / f / 1e3
+plot(dh$distance[-1], v, type="l", xlab="Distance [km]", ylab="Velocity [m/s]")
+grid()
+abline(h=0, col='red')

Exercise 6. Halifax sea-level during Hurricane Juan, near the end of September, 2003.

A web search will tell you that Hurricane Juan hit about midnight, 2003-sep-28. The first 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)
+# Focus on 2003-Sep-28 to 29th, the time when Hurricane Juan caused flooding
+plot(sealevel,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("Juan", at=as.POSIXct("2003-09-29 04:00:00", tz="UTC"), col="red")

Exercise 7. Plot the de-tided Halifax sea level during Autumn 2003, to see whether Hurricane Juan is visible. (Hint: use predict with the results from tidem.)

- +
library(oce)
+data(sealevel)
+m <- tidem(sealevel)
+oce.plot.ts(sealevel[['time']], sealevel[['elevation']] - predict(m),
+            ylab="Detided sealevel [m]", 
+            xlim=c(as.POSIXct("2003-09-20"), as.POSIXct("2003-10-08")))

The spike reveals a surge of about 1.5m, on the 29th of September, 2003.

Exercise 8. Manipulating flags.

- +
library(oce)
+d <- read.oce(system.file("extdata", "d201211_0011.cnv", package="oce"))
+
+## The following is in the .cnv file:
+# Processing Notes:  Flag word has 3 columns: Temperature, Conductivity, and Oxygen
+# Processing Notes:  Flag_code:  0 = no QC; 2 = good; 6 = interpolated, or replaced by dual sensor or upcast value;
+
+flag <- d[['flag']]
+flag1 <- floor(flag/100)
+flag2 <- floor((flag-100*flag1)/10)
+flag3 <- floor((flag-100*flag1-10*flag2))
+d[["flags"]]$temperature <- flag1
+d[["flags"]]$conductivity <- flag2
+d[["flags"]]$oxygen <- flag3