Introduction
Smoothing hydrographic profiles with conventional timeseries methods is problematic for two reasons: (a) the data are commonly not equispaced in depth and (b) the data seldom lack trends in depth. The issues and their solutions are illustrated without much discussion here.
The first step in making the graph shown above is to load the oce
library and create a function that measures daylength by finding sunrise and sunset times. Note that uniroot()
, which is used to find times of zero solar altitude, needs lower and upper limits on t
, and these are calculated by looking back and then forward a halfday. This works well for application to Halifax, but in other timezones other offsets would be needed. Interested readers might want to devised a method based on the longitude, which can be transformed into a timezone.
Methods
The first step is to load the data and extract dependent and independent variables, here S and p.

library(oce)
library(signal)
data(ctd)
S < ctd[['salinity']]
p < ctd[['pressure']]

Next, one must create equispaced data, for filtering to make any sense at all.

dp < median(diff(p))
pp < seq(min(p), max(p), dp)
S0 < approx(p, S, pp)$y

Now comes something that must be established by the actual task at hand: setting the critical frequency for the filter (as a ratio to Nyquist). In this case, a somewhat arbitrary frequency is selected, and then both firstorder and secondorder filters are created.

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

Now, set up for a threepanel plot.

par(mfrow=c(1, 3))

For the lefthand panel, show the raw data in black, and the two filters in red and blue.

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)

For the middle panel, detrended the profile, and then filter.

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)

For the righthand panel, use a smoothing spline instead of a filter.

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)

Results
Filtering a nondetrended profile is a generally a bad idea. There is almost always a zerooffset problem, and also most properties vary dramatically with depth, so detrending is required as well as zero offsetting. The advantage of detrending is illustrated in the lefthand and middle panels.
Smoothing splines provide an attractive alternative to filtering, especially in the notuncommon cases in which derivative are required. However, a disadvantage of splines is that there is no simple way to describe the method to those who are not familiar with them. In some branches of the literature, splines will be understood by readers, and in others, they will be a mystery that will waste a review cycle for the education of referees. There is also a problem in describing the spline simply in terms that have physical meaning. For example, spline smoothness is here controlled by setting the df
argument, but this has no simple analogy to physics, as perhaps the halfpower frequency of a filter might (in terms of spectral models of finestructure, say).
Rbloggers.com offers daily email updates about R news and tutorials on topics such as: Data science, Big Data, R jobs, visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...