Aligning Tide Gauge Stations

I’ve spoken before about my new way of aligning tide gauge stations. Maybe it’s time for me to outline some of the details, and share my program (written in R) for doing the computation.

Warning: there’s math ahead. If you want to use the program but you really don’t want to read a lot of equations, skip ahead to “How to Run the Program” and have some fun.

Velocity-adjustment Weighted alignment

Tide gauge data certainly give us information about the height of the sea surface (SSH). But their data are relative to a completely arbitrary “zero point” so to combine data from different stations, we need to estimate how far they’re offset from each other (we’ll call it the “offset”). If we have $N$ stations, there will be a separate offset for each. We’ll call them $\mu_A$ , where the index $A$ runs from 1 to $N$ (the number of stations we’re combining).

The idea is that the data at time $t_\alpha$ (where the index $\alpha$ runs from 1 to the number of observation times) isn’t just some oceanic signal, it also includes the offsets, and there’s noise in the measurement to boot, so

$x_{\alpha A} = \gamma_\alpha + \mu_A + \varepsilon_{\alpha A}$ .

Where $x_{\alpha A}$ is the observation from station $A$ at time $t_\alpha$ , $\gamma_\alpha$ is the oceanic signal at time $t_\alpha$ , $\mu_A$ is the offset for station A, and $\varepsilon_{\alpha A}$ is the noise value at time $t_\alpha$ from statino $A$ .

We can estimate the offsets by defining a quantity, the sum of squares of differences of each observation from the momentary “global” average

$Q = \Sigma_{stations~A} ~ \Sigma_{times~t_\alpha} ~ ( x_{\alpha A} - \gamma_\alpha - \mu_A )^2$ .

We then find the quantities which minimize this, which means finding not only the offsets $\mu_A$ , but the “global average” $\gamma_\alpha$ as well.

This is my old way of doing things, which I’ve sometimes called the “Berkeley method.” In my opinion it’s the best way to align many things, including temperature data from different stations.

But tide gauges are different. There’s another complication called vertical land movement (VLM). If the land rises or falls, the sea appears to fall or rise because tide gauges don’t measure absolute sea level, they measure relative sea level.

A saving grace about VLM is that at most stations, it goes at a steady pace. If at station $A$ the land is falling at a rate $\beta_A$ , then relative sea level will show an artificial rise at a rate of $\beta_A$ . Our new model is that the data at a given time $t_\alpha$ , at a given station $A$ , is given by

$x_{\alpha A} = \gamma_{\alpha} + \mu_A + \beta_A t_\alpha + \varepsilon_{\alpha A}$ .

Where $x_{\alpha A}$ is the observation from station $A$ at time $t_\alpha$ , $\gamma_\alpha$ is the oceanic signal at time $t_\alpha$ , $\mu_A$ is the offset for station $A$ , $\beta_A$ is the negative of the rate of vertical land movement at station $A$ , and $\varepsilon_{\alpha A}$ is the noise value at time $t_\alpha$ from statino $A$ .

Now we minimize a quantity similar to the one before, but including the effect of VLM:

$Q = \sum_A \sum_\alpha ( x_{\alpha A} - \gamma_\alpha - \mu_A - \beta_A t_\alpha )^2$ .

In addition to finding not only the offsets $\mu_A$ and the “global average” $\gamma_\alpha$ , we also find the VLM-induced rates $\beta_A$.

I have most recently introduced one additional complication. Some stations show a larger noise level than others, and I wanted to use a weighted mean for the global average, to take this into account. The new model requires the quantity to be minimized is

$Q = \sum_A \sum_\alpha \frac{1}{\sigma_A^2} ( x_{\alpha A} - \gamma_\alpha - \mu_A - \beta_A t_\alpha )^2$ ,

where $\sigma_A$ is the noise level associated with station $A$ .

How to Run the Program

vw.align is an R program to align records from different stations which do not all have observations at the same set of times. The basic inputs are:

tall — All the times for all the observations. There will be identical values
when multiple stations report at the same time.
xall — All the observations.
part — Station identifier for each observation. This can be a number or a name, as long as it’s unique for each individual station.

Optional paramters:

sloped — Should the alignment include a separate slope for each station?
Defaults to TRUE.
weighted — Should the computation use weighting (to downweight stations with
large variance)? Defaults to FALSE.
offres –Resolution in station offsets required for convergence. Defaults to 0.001.
Lower means more resolution and takes longer to compute; higher means less
resolution and computes faster.
full — Should it report input data and “adjusted” data for individual stations?
Defaults to FALSE.
tround — Rounding for input times, before matching times to identify when
stations are reporting the “same” time. Defaults to NULL.
xround — Rounding for output values. Defaults to NULL

The output is a list containing:

data — Results, including time, averaged value, its standard error,
and the number of stations reporting for each time.
offset — Offset for each individual station
slopes — Slope for each individual station (if used)
weights — Weights for each individual station (if used)

In addition, there are two further entries if the option “full” is set to “TRUE”

raw — Raw input data
adj — Station data corrected for offset and individual slope

Here is a tiny R script to run an example for you:

source(“vw_align.r”)

t1 = 0:999; x1 = rnorm(1000); p1 = rep(“A”,1000)
t2 = 100:899; x2 = rnorm(800)+.01*t2; p2 = rep(“B”,800)
t3 = 100:899; x3 = rnorm(800,,10)-.01*t3; p3 = rep(“C”,800)

plot(t1,x1,ylim=range(x1,x2,x3))
lines(t2,x2,col=”red”)
lines(t3,x3,col=”blue”)

tall = c(t1,t2,t3)
xall = c(x1,x2,x3)
part = c(p1,p2,p3)

q1 = vw.align(tall,xall,part,sloped=FALSE)
q2 = vw.align(tall,xall,part)
q3 = vw.align(tall,xall,part,weighted=TRUE,full=TRUE)

t = q1$d$t; x = q1$d$x; y = q2$d$x; w = q3$d$x

plot(t,x,type=”l”,col=”cornsilk4″)
lines(t,y,col=”red”)
lines(t,w,col=”blue”)

The Program

Wordpress kept telling me “this file type is not allowed for security reasons.” So here it is as inline text; copy and paste into your own file (which I call “vw_align.r”)

##################################
##################################
# vw.align
# combine station data
# t is time
# x is data
# part is station ID
# offres is resolution for offsets
# tround is rounding digits for t
# xround is rounding digits for x
#
# with options to use
# using variable slope by station
# and variable weight by station
##################################
##################################
vw.align = function(t,x,part,
sloped=TRUE,weighted=FALSE,full=FALSE,
tround=NULL,xround=4,offres=.001){
##########################
# make “part” a factor
##########################
part = as.factor(part)
nparts = length(levels(part))
if (nparts < 2){
zout = data.frame(t,x)
zout = zout[order(t),]
zout = list(data=zout,offsets=0)
} else{
#######################
# form giant data frame
#######################
for (jpart in 1:nparts){
tt = t[part == levels(part)[jpart]]
xx = x[part == levels(part)[jpart]]
ndx = order(tt)
xx = xx[ndx]
tt = tt[ndx]
if (length(tround) > 0){
tt = round(tt,tround)
}
z1 = data.frame(tt,xx)
names(z1) = c(“t”,levels(part)[jpart])
if (jpart == 1){
zz = z1
} else{
zz = merge(zz,z1,by=1,all=T)
}
}
###############################
# iteratively determine offsets
# and merged values
###############################
zdat = zz[,2:(nparts+1)]
zdat = as.matrix(zdat)
ntimes = dim(zz)[1]
tau = zz[,1]
tau = tau-mean(tau)
################
# initial values
################
mu = rep(0,nparts) # station A: zero point offset
beta = rep(0,nparts) # station A: slope offset
wate = rep(1,nparts) # station A: statistical weight
gam = rep(NA,ntimes) # time alpha: global mean
gam2 = gam # UNweighted mean when gam is weighted
dmu = 9999; dbeta = 9999
######################
# do until convergence
######################
while (dmu > offres){
for (jj in 1:ntimes){
xx = zdat[jj,]
xx = xx-mu-beta*tau[jj]
if (weighted){
gam[jj] = weighted.mean(xx,wate,na.rm=TRUE)
gam2[jj] = mean(xx,na.rm=TRUE)
} else{
gam[jj] = mean(xx,na.rm=TRUE)
}
}
oldmu = mu
if (sloped){
for (jj in 1:nparts){
xx = zdat[,jj]-gam
xxfit = lm(xx~tau)
mu[jj] = xxfit$coef[1]
beta[jj] = xxfit$coef[2]
}
} else{
for (jj in 1:nparts){
xx = zdat[,jj]-gam
mu[jj] = mean(xx,na.rm=TRUE)
}
}
mu = mu-median(mu)
beta = beta-median(beta)
#####################
# compute new weights
#####################
if (weighted){
for (jj in 1:nparts){
xx = zdat[,jj]-gam2-mu[jj]-beta[jj]*tau
wate[jj] = 1/var(xx,na.rm=TRUE)
}
}
#########################
# difference in mu values
#########################
dmu = mean(abs(mu-oldmu))
}
############
# final mean
############
num = rep(NA,ntimes)
se = rep(NA,ntimes)
std.dev = rep(NA,ntimes)
for (jj in 1:ntimes){
xx = zdat[jj,]
xx = xx-mu-beta*tau[jj]
if (weighted){
gam[jj] = weighted.mean(xx,wate,na.rm=TRUE)
} else{
gam[jj] = mean(xx,na.rm=TRUE)
}
xx = xx[is.finite(xx)]
num[jj] = length(xx)
std.dev[jj] = sd(xx)
}
se = std.dev/sqrt(num)
zout = data.frame(x=gam,num,se,std.dev)
if (length(xround) > 0){
zout = round(zout,xround)
}
zout = data.frame(t=zz[,1],zout)
if (full){
zd = data.frame(zdat)
names(zd) = levels(part)
zd2 = zd
for (jj in 1:nparts){
zd2[,jj] = zd2[,jj]-mu[jj]-beta[jj]*tau
}
zout = list(data=zout,offsets=mu,slopes=beta,weights=wate,
raw=zd,adj=zd2)
} else{
zout = list(data=zout,offsets=mu,slopes=beta,weights=wate)
}
}
zout
}

Final note: If you’d like to make a donation to this blog, this would be a particularly good time to follow the link below.

This blog is made possible by readers like you; join others by donating at My Wee Dragon.

21 responses to “Aligning Tide Gauge Stations”

ecoquant | October 21, 2019 at 9:04 pm |

Question, and I haven’t gone through the code to see …

Do you do anything in the resulting expansion to keep, for instance, the $\gamma_{\alpha}^{2}$ and $\beta_{A}^{2}t_{\alpha}^2$ terms identifiable? In other words, it seems when minimizing $Q$ there could be a path along the surface, given constant $A$ , where $\gamma_{\alpha}^{2}$ and $\beta_{A}^{2}t_{\alpha}^2$ trade off against one another keeping $Q$ constant. There are other cross terms where the same might be true.

In other applications, identifiability is often handled by a change of variables where, using the notation here, reference $\gamma_{\alpha_{0}}$ is chosen and then others are re-expressed as $\gamma_{\alpha} = \Delta \gamma_{\alpha} + \gamma_{\alpha_{0}}$ , but I’m not sure how that would work here?

Are there separate constraints on these? Priors perhaps?

I probably should just sit down and study the code, but if this kind of thing is in there it would be good to have some explanation in the text.

Just my two cents.

Thanks.

[Response: I should have provided more details. There are ambiguities, because adding a constant to all the offsets while subtracting it from the global signal, gives the Q value unchanged. Also, adding a constant to all the slopes leaves Q unchanged.

Fortunately one gets around that by imposing simple constraints. I originally constrained the constants (offset and slope) to be zero for the first station, but now I constrain the median offset and median slope to be zero.]
Matthias U | October 23, 2019 at 8:12 am |

Owch. Is it possible to do this in a way that doesn’t toast all the indents?

[Response: Not easily on wordpress. Let me see what I can come up with.]
bindidon | October 23, 2019 at 2:18 pm |

Tamino

Isn’t it possible to load txt files into e.g. Libre Office Writer, from which you may obtain a pdf file, and finally to upload the latter via Google Drive or similar?

I know it’s ugly but it works, like here:
https://drive.google.com/file/d/1kLZZQH-zQjkWMvwT72JxAM34cLPOUjMr/view
AJ | October 23, 2019 at 3:50 pm |

Yet another reconstruction. This one employing kriging:
https://climateadj.wordpress.com/2019/05/14/kriged-sea-level-rise/
- bindidon | October 23, 2019 at 4:12 pm |
  
  AJ
  
  Could you please upload the list of those PMSL stations which survived all your successive filtering operations (in a form allowing me to quickly extract them out of PMSL’s station file (their id is enough for that)? I’m highly interested in it.
  
  I start from the assumption that you too use only RLR records.
  - AJ | October 24, 2019 at 1:45 pm |
    
    With regard to filtering, yes it’s always a question if you’re filtering noise or information. More high quality data would be nice. Most stations were filtered due to lack of GPS. I also tried including stations with a low seismic hazard, without VLM adjustment, but this led to a significant variogram nugget.
- bindidon | October 23, 2019 at 7:48 pm |
  
  AJ [2]
  
  My reason to ask for your final PMSL station subset is the amazing difference between your gauge selection:
  
  and the available set of all stations giving valuable data even if for short periods:
  https://drive.google.com/file/d/1iCIoZqp0ImvktVLUkJ0yNet1DVVafhzG/view
  
  It is clearly visible that your successive flitering steps eliminate a great part of the Globe’s coastal corners, what might create a far bigger flaw than the sum of those you tried to get rid of.
  - AJ | October 24, 2019 at 5:13 am |
    
    There is a link to source code and intermediary files provided on the page
    The file rlr-obs-df.rds is a data frame which contains psmsl station id.
- bindidon | October 24, 2019 at 9:14 am |
  
  AJ
  
  Thanks, that I had imagined already, but… not everybody on Earth uses R :-(.
  I’m a good old C-imple-minded programmer.
- bindidon | October 24, 2019 at 7:44 pm |
  
  AJ [4]
  
  I forgot above all to mention that the file contents of the directory
  https://drive.google.com/file/d/18KXQoi8qLXYMAolqb5gYlMDgsLKKQJ1c/view
  
  and those of its subdirectory ‘source’ are not accessible.
  - AJ | October 24, 2019 at 11:37 pm |
    
    Crap. I’ll fix that next week when I have access to my laptop.
  - AJ | October 29, 2019 at 5:41 pm |
    
    Bindidon… I updated the link and tried in incognito mode. I got a network error message, but was still able to download. Not sure the reason for the difficulties. Thought for sure I had tested it previously.
  - bindidon | October 29, 2019 at 6:28 pm |
    
    AJ [5]
    
    Sorry, still no change on my side. Only the ‘source’ directory is accessible within ‘slr-kriged.zip’, but no file contents are either.
    
    Please upload the gauge list as a single pdf file.
  - AJ | October 29, 2019 at 7:32 pm |
    
    Here’s the updated link:
    https://www.dropbox.com/s/mn1et2lv7srb7l8/slr-kriged.zip?dl=0
  - bindidon | October 30, 2019 at 8:49 pm |
    
    AJ [end]
    
    “.rds files can’t be previewed.”
    I give it up: it’s too boring to help you understanding even simplest things.
  - ecoquant | October 31, 2019 at 2:12 am |
    
    @bindidon,
    
    You can be as intransigent about SLR and climate as you want, but I draw the line at insulting R‘s .rds files!
  - AJ | October 30, 2019 at 11:40 pm |
    
    bindidon [out]
  - bindidon | October 31, 2019 at 7:40 pm |
    
    ecoquant
    
    “… but I draw the line at insulting R‘s .rds files!”
    
    Wow. Sir, I propose that you read the sequence a bit more carefully before writing such unnecessary kauderwelsch.
    
    Why the heck should I little layman insult a worldwide known programming language or one of its technical components?
  - ecoquant | October 31, 2019 at 11:07 pm |
    
    @Bindidon,
    
    Wow. Sir, I propose that you read the sequence a bit more carefully before writing such unnecessary kauderwelsch.
    
    Because I was making a joke.
    
    $\textttt{.rds}$ files are actually pretty handy.
ecoquant | October 24, 2019 at 1:25 am |

The code could be tossed into a Github as it popular these days. I just put it into a Google Drive directory … Google even gives you some version control.
bindidon | October 24, 2019 at 12:28 pm |

It would be very interesting to obtain a sorted trend chart for PMSL stations like this one

https://drive.google.com/file/d/1MKN4zBzgKIwM9BPBOS6M7Vwoqzm7ahSi/view

but then including all corrections for vertical land movement.

The chart above compares, for the 679 (of 1513) PMSL stations having presented sufficient data to compute anomalies for the reference period 1993-2013, the trends for (1) their entire lifetime and for (2) the reference period.

Trends below -10 mm/yr and above +10 mm/yr were excluded in this chart (this was of course an arbitrary decision).

The life time average trend was 2.1 mm/yr regardless of trend exclusion; the average trend for the reference period was 3.1 mm/yr without, and 2.7 mm/yr with trend exclusion.

	kinimod on A High Schooler’s Take o…
	geoffbeacon on A High Schooler’s Take o…
	Bob Loblaw on Not Even Wrong?
	Bob Loblaw on Not Even Wrong?
	Tom Passin on Not Even Wrong?
	convictstreak on A High Schooler’s Take o…
	Bob Loblaw on Not Even Wrong?
	Cliff Mass on Not Even Wrong?
	jgnfld on Not Even Wrong?
	Bob Loblaw on Not Even Wrong?
	kinimod on A High Schooler’s Take o…
	jgnfld on Not Even Wrong?
	geoffbeacon on A High Schooler’s Take o…
	climatefreak on A High Schooler’s Take o…
	Lowlander on Not Even Wrong?

Aligning Tide Gauge Stations

21 responses to “Aligning Tide Gauge Stations”

Support Your Global Climate Blog

Recent Comments

Recent Posts

Buy the book

astronomy

Blogroll

Global Warming

mathematics

Aligning Tide Gauge Stations

Share this:

Like this:

Related

21 responses to “Aligning Tide Gauge Stations”

Support Your Global Climate Blog

Recent Comments

Recent Posts

Buy the book

astronomy

Blogroll

Global Warming

mathematics