I am a physical oceanographer interested in how ocean water is mixed and transformed. I am currently working as Research Scientist at the Bedford Institute of Oceanography in Halifax, Nova Scotia.

Recent Posts:

Calculating buoyancy frequency for argo/section objects using the apply() family

The most recent CRAN release of oce includes some nice new functionality for reading and converting argo objects (see http://www.argo.ucsd.edu/ for more information about the fantastic Argo float program). One question that arose out of this increased functionality was how to calculate (also known as the buoyancy or Brunt-Väisälä frequency) for such objects.

Buoyancy frequency

The definition of is:

where is the acceleration due to gravity, is the fluid density, and is the vertical coordinate. Essentially describes the vertical variation of fluid density (also known as “stratification”).

Calculating for regular ctd objects is easily accomplished with the function oce::swN2(). A caution: readers are encouraged to read the documentation carefully, as the details of the actual calculation can have important consequences when applied to real ocean data.

for station objects

For the case of a station object (which is essentially a collection of ctd stations), the most straightforward way to calculate is to use the lapply() function to “apply” the swN2() function to each of the stations in the object. An example:

section[['station']] <- lapply(section[['station']],
                               function(x) oceSetData(x, 'bvf', swN2(x)))

The line with the lapply() command takes the list of stations from the section object, and evaluates each of the resulting ctd objects using the oceSetData() function to add the result of swN2() back into the station @data slot.

If we wanted to make a nice plot of the result, we could do:

col <- colorRampPalette(c('white', rev(oceColorsViridis(32))))
plot(section, which='bvf', ztype='image', zbreaks=seq(0, 1e-4, 0.5e-5), zcol=col)

plot of chunk stationplot where I’ve defined a custom colormap just for the fun of it.

for argo objects

In an argo object, the default storage for the profiles is a matrix, rather than a list of ctd objects. To calculate and make a plot, the simplest approach would be to use as.section() to convert the argo object to a section class object and then do as above. However, having the field as a matrix allows for greater flexibility in plotting, e.g. using the imagep() function, so one might want to calculated in a manner consistent with the default argo storage format.

Let’s load some example data from the argo dataset included in oce:

argo <- argoGrid(argo)

Note that I’ve gridded the argo fields so the matrices are at consistent pressure levels. Now we create a function that can be applied to each of the matrix columns, to calculate from a single column of the density matrix:

N2 <- function(x) {
    swN2(argo[['pressure']][,1], x)

Now we use the above function N2 to calculate buoyancy frequency and add it back to the original object, like:

argo <- oceSetData(argo, 'N2', apply(swSigmaTheta(argo), 2, N2) )

Note that because of the difference between the “list” and “matrix” approach, the oceSetData() occurs outside of the apply(). Also note the second argument in the apply() call, which specifies to apply the N2() function along the 2nd dimension of the density matrix, i.e. along columns.

Now, lets make a sweet plot of the N2 field using imagep()!

p <- argo[['pressure']][,1]
t <- argo[['time']]
imagep(t, p, t(argo[['N2']]), flipy=TRUE, col=col, breaks=seq(0, 1e-4, 1e-6),
       ylab='pressure [dbar]',
       zlab=expression(N^2~group('[', 1/s^2, ']')), zlabPosition = 'side',
       mar=c(3, 3, 2, 0.5))

plot of chunk argoplot

A thing of stratified beauty, if I do say so myself.

R tutorial at the 2016 AGU Ocean Sciences Meeting

Today I had the pleasure of presenting a talk about R during one of the tutorial sessions at the AGU Ocean Sciences Meeting in New Orleans. I made a deliberate point of saying that my main message was more like: “R is really cool, and here’s why”, rather than: “You should all stop using Matlab”. Being divisive doesn’t help anybody work better.

I was quite surprised at the number of people who raised their hand when I asked if they use R as their main analysis tool – perhaps about a third of the (surprisingly full) room!

Anyway, see below for the pdf of my talk and the Rnw source file. Note that the Rnw has links to some external images so you won’t be able to knitr() it for yourself, but feel free to use the slide content.



A Valentine's day map

In the spirit of Valentine’s day, I made a map that uses the beautifully weird bonne projection. The code uses the proj4 facilities in the oce package.


msgtop <- "You mean the WORLD to me"
msglon <- 130
msglat <- -30

mapPlot(coastlineWorld, col='pink', proj='+proj=bonne +lat_1=85',
        xlim=c(-13700000, 13700000),
        ylim=c(-18079625, 8008557))
mapGrid(longitude = c(-180, 180), lwd=4, col=2, latitude=NULL)
text(0, 8005400, msgtop, cex=2, font=3)
nc <- nchar(msg)
lon <- seq(-msglon, msglon, length.out = nc)
rot <- seq(-45, 45, length.out = nc)
for (i in 1:nc) {
    mapText(lon[i], msglat, substr(msg, i, i), cex=2.2, font=2, srt=rot[i])
    if (substr(msg, i, i+1) == 'ES')  mapText(1.02*lon[i], 0.8*msglat, "`", cex=2.2, font=2, srt=rot[i])

plot of chunk heartmap

oce package version 0.9-18 is released!

Today we released a new version of oce, and it has been uploaded to CRAN. Only the source version is available as of the time of writing, but binary versions for all platforms should become available in the next few days. As always, the best way to install the package is to do:


at an R prompt.

Then you can do stuff like:

## Loading required package: methods
## Loading required package: gsw
## Loading required package: testthat
cm <- colormap(levitus$SSS, col=oceColorsSalinity, breaks=seq(30, 37, 0.5))
drawPalette(colormap=cm, drawTriangles = TRUE)
mapPlot(coastlineWorld, projection='+proj=stere +lat_0=90 +lon_0=-60',
        longitudelim=c(-180, 180), latitudelim=c(50, 90), col='grey')
with(levitus, mapImage(longitude, latitude, SSS, colormap=cm, filledContour = TRUE))
mapPolygon(coastlineWorld, col='grey')
title('Sea Surface Salinity')

plot of chunk unnamed-chunk-2

New features

The previous version of oce was uploaded to CRAN about 9 months ago. In the meantime, we’ve fixed lots of bugs and added even more improvements. A quick look at the NEWS file gives a summary of the enhancements:

- improve plot.coastline() and mapPlot()
- add support for G1SST satellite
- all objects now have metadata items for units and flags
- ctdTrim() method renamed: old A and B are new A; old C is new B
- support more channels and features of rsk files
- as.adp() added
- convert argo objects to sections
- makeSection() deprecated; use as.section() instead
- read.adp.rdi() handles Teledyne/RDI vsn 23.19 bottom-track data
- geodXyInverse() added; geod functions now spell out longitude etc
- read.odf() speeded up by a factor of about 30
- add colour palettes from Kristen Thyng's cmocean Python package
- as.oce() added
- rename 'drifter' class as 'argo' to recognize what it actually handles
- add oceColorsViridis()
- interpBarnes() has new argument 'pregrid'
- binMean2D() has new argument 'flatten'
- data(topoWorld) now has longitude from -179.5 to 180
- ODF2oce() added
- read.odf() handles more data types
- read.adp.rdi() reads more VmDas (navigational) data
- ITS-90 is now the default temperature unit
- ctd objects can have vector longitude and latitude
- logger class renamed to rsk
- bremen class added
- coastlineCut() added
- rgdal package used instead of local PROJ.4 source code
- mapproj-style map projections eliminated

Some of the best additions (in my opinion) are:

  • addition of a units metadata field for objects
  • more tools for working with argo objects (which have been renamed from drifter)
  • new color palettes, particularly the cmocean palettes from python
  • renaming of the logger class to rsk for data from RBR instruments
  • lots of new features for rsk objects, including ability to convert to ctd objects with as.ctd()
  • more robust projection handling using the rgdal package
  • coastlineCut() for helping with UHL (ugly horizontal lines) in map projections


The vignette has also been updated somewhat, using an Rmarkdown source and with some new examples.

Help us out!

You can help! If you find any bugs or have any requests, please open an Issue on the Github development page:


oce has benefited immensely from some great requests recently, and we’d like to keep the momentum going!

R/oce at AGU Ocean Sciences 2016

If you’re going to be at the AGU Ocean Sciences meeting in New Orleans in a couple weeks, make sure you come by to check out the R/oce tutorial I’ll be doing. It’s on Wednesday February 24th at 3:00pm in room R03. See you there!

An R function to shift vectors by a specified lag

Quite of bit of my work involves looking at “shifts” between two time series. There are lots of reasons why shifts are interesting, including such things as:

  • phase differences in the tides at two different locations,

  • physical separation between sensors on a profiling instrument, and

  • clock drifts between two logging sensors.

To accomplish the actual shifting of the vectors (I’m not going to discuss here how to determine the amount by which the series should be shifted, since that depends on the parameters of the problem), I created the following function:

shift <- function(x, lag) {
    n <- length(x)
    xnew <- rep(NA, n)
    if (lag < 0) {
        xnew[1:(n-abs(lag))] <- x[(abs(lag)+1):n]
    } else if (lag > 0) {
        xnew[(lag+1):n] <- x[1:(n-lag)]
    } else {
        xnew <- x

The function takes as input the vector x and a specified integer lag, and shifts the input series by that amount. NA’s are added to either the beginning or the end (depending on the sign of lag) to pad the shifted vector to be the same length as the input. Note that lag is defined so that a positive lag shifts x “to the right”, i.e. moves values to a higher index.