# R package: oce

Dan Kelley is one of those heroes that creates and maintain software libraries for oceanographers. His language of choice is R, I particularly do not like R syntax. I use it only when I have really have to.

However, R does have some amazing libraries, like Dan's oce package.

This post is just a copy-and-paste of few examples from Dan's oce blog, to show how to use the oce library with the new IPython/Jupyter R kernel.

Every example has a link to the original post and I encourage you to go check it out.

# Installing oce

library(devtools)

options(unzip = "unzip")

install_github('ocedata', 'dankelley', 'master')
install_github('oce', 'dankelley', 'develop')


# Day length calculation

In [1]:
library(oce)

daylength <- function(t, lon=-38.5, lat=-13)
{
t <- as.numeric(t)
alt <- function(t)
sunAngle(t, longitude=lon, latitude=lat)$altitude rise <- uniroot(alt, lower=t-86400/2, upper=t)$root
set <- uniroot(alt, lower=t, upper=t+86400/2)$root set - rise } t0 <- as.POSIXct("2014-01-01 12:00:00", tz="UTC") t <- seq.POSIXt(t0, by="1 day", length.out=1*356) dayLength <- unlist(lapply(t, daylength)) par(mfrow=c(2,1), mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0)) plot(t, dayLength/3600, type='o', pch=20, xlab="", ylab="Day length (hours)") grid() solstice <- as.POSIXct("2013-12-21", tz="UTC") plot(t[-1], diff(dayLength), type='o', pch=20, xlab="Day in 2013", ylab="Seconds gained per day")  packageStartupMessage in packageStartupMessage(gettextf("Loading required package: %s", : Loading required package: mapproj packageStartupMessage in packageStartupMessage(gettextf("Loading required package: %s", : Loading required package: maps packageStartupMessage in packageStartupMessage(gettextf("Loading required package: %s", : Loading required package: proj4  # Box models in R In [2]: library(deSolve) IC <- 0 parms <- list(A=1e6, gamma=10) sperday <- 86400 # seconds per day times <- seq(0, 10*sperday, length.out=200) F <- function(t) { ifelse (t > sperday & t < 2*sperday, 1, 0) } DE <- function(t, y, parms) { h <- y[1] dhdt <- (F(t) - parms$gamma * h) / parms$A list(dhdt) } sol <- lsoda(IC, times, DE, parms) par(mfrow=c(2,1), mar=c(3,3,1,1), mgp=c(2,0.7,0)) h <- sol[,2] Day <- times / 86400 plot(Day, F(times), type='l', ylab="River input [m^3/s]") plot(Day, h, type='l', ylab="Lake level [m]")  # oce runlm function In [3]: library(oce) x <- 1:100 y <- 1 + x/100 + sin(x/5) yn <- y + rnorm(100, sd=0.1) L <- 4 calc <- runlm(x, y, L=L, deriv=0) plot(x, y, type='l', lwd=7, col='gray') points(x, yn, pch=20, col='blue') lines(x, calc, lwd=2, col='red')  # N2 with runlm() In [4]: library(oce) data(ctd) rho <- swRho(ctd) z <- swZ(ctd) drhodz <- runlm(z, rho, deriv = 1) g <- 9.81 rho0 <- mean(rho, na.rm = TRUE) N2 <- -g * drhodz/rho0 plot(ctd, which = "N2") lines(N2, -z, col = "blue") legend("bottomright", lwd = 2, col = c("brown", "blue"), legend = c("spline", "runlm"), bg = "white")  In [5]: # Alter next three lines as desired; a and b are watermasses. Sa <- 30 Ta <- 10 Sb <- 40 library(oce) # Should not need to edit below this line rho0 <- swRho(Sa, Ta, 0) Tb <- uniroot(function(T) rho0-swRho(Sb,T,0), lower=0, upper=100)$root
Sc <- (Sa + Sb) /2
Tc <- (Ta + Tb) /2
## density change, and equiv temp change
drho <- swRho(Sc, Tc, 0) - rho0
dT <- drho / rho0 / swAlpha(Sc, Tc, 0)

plotTS(as.ctd(c(Sa, Sb, Sc), c(Ta, Tb, Tc), 0), pch=20, cex=2)
drawIsopycnals(levels=rho0, col="red", cex=0)
segments(Sa, Ta, Sb, Tb, col="blue")
text(Sb, Tb, "b", pos=4)
text(Sa, Ta, "a", pos=4)
text(Sc, Tc, "c", pos=4)
legend("topleft",
legend=sprintf("Sa=%.1f, Ta=%.1f, Sb=%.1f  ->  Tb=%.1f, drho=%.2f, dT=%.2f",
Sa, Ta, Sb, Tb, drho, dT),
bg="white")


# Inferring halocline depth

In [6]:
findHalocline <- function(ctd, deltap=5, plot=TRUE)
{
S <- ctd[['salinity']]
p <- ctd[['pressure']]
n <- length(p)
## trim df to be no larger than n/2 and no smaller than 3.
N <- deltap / median(diff(p))
df <- min(n/2, max(3, n / N))
spline <- smooth.spline(S~p, df=df)
SS <- predict(spline, p)
dSSdp <- predict(spline, p, deriv=1)
H <- p[which.max(dSSdp$y)] if (plot) { par(mar=c(3, 3, 1, 1), mgp=c(2, 0.7, 0)) plotProfile(ctd, xtype="salinity") lines(SS$y, SS$x, col='red') abline(h=H, col='blue') mtext(sprintf("%.2f m", H), side=4, at=H, cex=3/4, col='blue') mtext(sprintf(" deltap: %.0f, N: %.0f, df: %.0f", deltap, N, df), side=1, line=-1, adj=0, cex=3/4) } return(H) } # Plot two panels to see influence of deltap. par(mfrow=c(1, 2)) data(ctd) findHalocline(ctd) findHalocline(ctd, 1)  Out[6]: [1] 4.292  Out[6]: [1] 4.721  # Smoothing CTD profiles In [7]: library(oce) library(signal) data(ctd) S <- ctd[['salinity']] p <- ctd[['pressure']] dp <- median(diff(p)) pp <- seq(min(p), max(p), dp) S0 <- approx(p, S, pp)$y

W <- dp / 2
f1 <- butter(1, W)
f2 <- butter(2, W)

par(mfrow=c(1, 3))

plotProfile(ctd, xtype="salinity", type='l')
S0f1 <- filtfilt(f1, S0)
S0f2 <- filtfilt(f2, S0)
lines(S0f1, pp, col='red')
lines(S0f2, pp, col='blue')
mtext("(a) ", side=3, adj=1, line=-5/4, cex=3/4)

plotProfile(ctd, xtype="salinity", type='l')
Sd <- detrend(pp, S0)
S1f1 <- filtfilt(f1, Sd$Y) + Sd$a + Sd$b * pp S1f2 <- filtfilt(f2, Sd$Y) + Sd$a + Sd$b * pp
lines(S1f1, pp, col='red')
lines(S1f2, pp, col='blue')
mtext("(b) ", side=3, adj=1, line=-5/4, cex=3/4)

spline <- smooth.spline(pp, S0, df=3/W)
S2 <- predict(spline)\$y
plotProfile(ctd, xtype="salinity", type='l')
lines(S2, pp, col='red')
mtext("(c) ", side=3, adj=1, line=-5/4, cex=3/4)

packageStartupMessage in packageStartupMessage(gettextf("\nAttaching package: %s\n", sQuote(package)), :
Attaching package: â€˜signalâ€™

packageStartupMessage in packageStartupMessage(msg): The following object is masked from â€˜package:oceâ€™:

decimate

packageStartupMessage in packageStartupMessage(msg): The following objects are masked from â€˜package:statsâ€™:

filter, poly



There are many more interesting posts, like the Gulf Stream center and the solstice. I was unable to run them due to data and/or difficulty to install certain R packages :-(

PS: Note that, due to a bug in the R-ipython kernel, a new figure is created every time R adds something new to the plot. There is a work around here.

In [2]:
HTML(html)

Out[2]:

This post was written as an IPython notebook. It is available for download or as a static html.