NH Sea Level Reconstruction

The following is my attempt to reconstruct the NH Mean Sea Level for the period 1930-2009.  I will describe the method used and compare my results with Church and White (2011) Global values.  I see that C&W applied a GIA adjustment to their data whereas I made no such attempt.

Method

Setup
  • Download PSMSL tide gauge data.
  • Download C&W data.
  • Create a file of sea surface area relative weighting by latitude.  This was done by taking a 1x1 gridded SST dataset, calculating the fraction of grid cells reporting for each latitude, and then multiplying by the latitudinal cosine.
Processing

  • First create a data frame of PSMSL data by concatenating the individual files and attaining the station id from the filenames.
  • Define the latitudinal bands to be calculated and a percentage of reporting periods threshold that each station must meet in order to be selected.  In my case I defined five bands as follows:

South LatitudeNorth LatitudeReporting Threshold
02530%
254080%
405080%
505590%
556090%

So for the 0-25 band, a station must report on 30% of the 1930-2009 period in order to make it into my sample.  This was an arbitrary choice.  Wider bands and lower thresholds were selected in the more southerly latitudes due to a dearth of station data.  The more northerly bands were narrower and a higher threshold required due to more data and GI effects.

  • For each band, select the sample and create a matrix of MSL by station and year.  For illustration purposes, here's an example:

     yr
stn   1930 1931 1932 1933 1934 1935 1936 1937 1938 1939 1940
  155 6984 6996 6991 6978 6943 6998 7003 6972 6980 6981 7057
  163 7047 7000 7047 7008 7009 7021 7036 7061 7018 7032 7054
  269 7001 7025 6991 7006 6985 6976 7016 7019 7014 7007 6969
  331 6965 6943 6968 7022 6997 7014 6998 7020 7012 7009 7030
  43  7006 7046 7043 7063 7055 7020 7017 7066 7075 7045 7076
  44  6909 6930 6921 6933   NA   NA   NA 6958   NA   NA 6983

  • Then take the year on year difference for each station, appending a zero for the first year:

      1931 1932 1933 1934 1935 1936 1937 1938 1939 1940
155 0   12   -5  -13  -35   55    5  -31    8    1   76
163 0  -47   47  -39    1   12   15   25  -43   14   22
269 0   24  -34   15  -21   -9   40    3   -5   -7  -38
331 0  -22   25   54  -25   17  -16   22   -8   -3   21
43  0   40   -3   20   -8  -35   -3   49    9  -30   31
44  0   21   -9   12   NA   NA   NA   NA   NA   NA   NA

  • Then calculate the average change for each year (excluding NA's) and calculate the cumulative sum:

       1931  1932  1933  1934  1935  1936  1937  1938  1939  1940 
 0.00  4.67  8.17 16.33 -1.27  6.73 14.93 28.53 20.73 15.73 38.13 

  • Then take the weighted average of the sampled bands.  The weighting was by the relative sea surface area of each band.

And here's the resulting plot:

And here are the residuals:

Discussion

My results were largely "consistent with" C&W's results without consideration for GIA and NH vs Global allowances.  My 1930-2009 slope is little higher but not significantly so (1.90mm/yr vs 1.82mm/yr).  My 1993-2009 slope is also a little higher (2.93mm/yr vs 2.80mm/yr).  The R^2 value is 0.95.  I note that my time series is noisier.  If I specify my threshold as 10% for all bands the results reverse with my 1930-2009 slope being a little lower.

I probably should have labelled my results as a Method or Reconstruction instead of Model.  I'd change that but it's early Saturday afternoon here and the sun just came out and also I'm just too damn lazy.  I tried to come up with a better name for my method [Sea Height I Tried makes for an interesting acronym:)], but just settled on using my initials.  

Source Code and Instructions

Data is available from the following sources:

Church White (2011)
http://www.cmar.csiro.au/sealevel/GMSL_SG_2011.html

PSMSL RLR Annual
http://www.psmsl.org/data/obtaining/rlr.annual.data/rlr_annual.zip

These files and my sea-area-weight-by-lat.txt file are attached to this page for your convenience.  Save/extract them in your working directory.

Source Code:

#
# This program is an attempt to reconstruct NH sea levels
# for the period 1930-2009.  Also this will be compared
# against Church and White (2011) Global Mean Sea Levels
# over the same period.
#

# read in self constructed table of sea surface area
# relative weights by latitude.

wgt_df = read.table("sea-area-weight-by-lat.txt",sep=",")
names(wgt_df) = c("lat","wgt")

# read in Church and White (2011) GMSL data

recon_df = read.table("CSIRO_Recons_gmsl_yr_2011.txt")
names(recon_df) = c("t","anom")
recon_df$t = floor(recon_df$t)

# read in the PSMSL station data

stn_df = read.table("rlr_annual/filelist.txt",sep=";",as.is=T,strip.white=T,quote="")
names(stn_df) = c("stn","lat","lon","name","cst_id","old_id","flag_ind")

# read in and concatenate the PSMSL station data including filename

files = list.files(path="rlr_annual/data",full.names=T)
rlr_df = do.call("rbind", lapply(files, function(fn) data.frame(stn=fn, read.table(fn,sep=";"))))

# attain the station id from the filename

stnstr = strsplit(as.vector(rlr_df$stn),".",fixed=T)
stnstr = sapply(stnstr, "[", 1)
stnstr = strsplit(stnstr,"/",fixed=T)
rlr_df$stn = sapply(stnstr, "[", 3)

names(rlr_df)[2:3] = c("t","anom")
rlr_df$anom[rlr_df$anom==-99999] = NA

##############

# Function: rlranom
# Parameters:
#
#    rlr_df - dataframe of PSMSL sea level data
#    stn_df - dataframe of PSMSL station data
#    wgt_df - dataframe of latitudinal sea surface area relative weights
#    yr_s   - starting year of period
#    yr_e   - ending year of period
#    lat_s  - southern boundary
#    lat_n  - northern boundary
#    pcut   - fraction of period that a station must report on to make the cut
#
# Returns the following in a list:
#
#    rlr_a - timeseries of yearly mean sea level anomaly (is zero at year zero)
#    mlat  - the mean latitude of the sampled stations
#    wgt   - calculated weight for sampled region
#    rlr_n - timeseries of the sample size for each year
#

rlranom = function(rlr_df, stn_df, wgt_df, yr_s, yr_e, lat_s, lat_n, pcut){

 yrs = yr_s:yr_e

 # attain stations within sampled region
 samp_df = subset(stn_df, lat > lat_s & lat <= lat_n & flag_ind == "N")

 # subset PSMSL sea level data of sample period, exclude NA's
 rlrx_df = subset(rlr_df,t %in% yrs & !is.na(rlr_df$anom))
 
 # find number of reporting period for each station
 agg_df = aggregate(rlrx_df$stn, by=list(stn=rlrx_df$stn), FUN=length)
 
 # calculate number of years in sample period (used for cutoff test)
 testcnt = length(yrs)
 
 # get sample stations, data for stations making cut and in sample region
 samp = agg_df$stn[agg_df$x/testcnt >= pcut & agg_df$stn %in% samp_df$stn]
 
 # subset PSMSL sea level data for sampled stations
 rlrx_df = subset(rlrx_df,stn %in% samp)
 
 # create a station x year sample matrix of NA's
 rlr_m = matrix(NA,length(samp),length(yrs))
 dimnames(rlr_m) = list(stn=samp,yr=yrs)

 # populate sample matrix
 rlr_m[cbind(rlrx_df$stn,rlrx_df$t)] = rlrx_df$anom

 # apply diff function to station data and bind a zero at year 0
 #   (note: probably retarded to bind the zero, but I like
 #          seeing the sample size in year zero in histogram
 #          plots one can produce - see below)
 
 rlr_d = cbind(rep(0,dim(rlr_m)[1]),t(apply(rlr_m,1,diff)))
 
 # calculate timeseries of number of sampled stations for each year
 rlr_n = apply(rlr_d,2,function(x)length(which(!is.na(x))))
 
 # calculate anomaly by first calculation the mean of the change
 # for each year and the calculating the cumulative sum.
 rlr_a = cumsum(apply(rlr_d,2,mean,na.rm=T))
 
 # calculate the mean latitude of the sampled stations
 mlat = mean(samp_df$lat)

 # calculate the sea surface area relative weight for sampled region
 wgt = sum(subset(wgt_df,lat > lat_s & lat <= lat_n)$wgt)
 
 # return list of calculated values
 list(rlr_a=rlr_a,mlat=mlat,wgt=wgt,rlr_n=rlr_n)
 
}

#####################
# MAIN PROCESS
#####################

# define sample period
yr_s = 1930
yr_e = 2009
yrs = yr_s:yr_e

# calculate the anomalies for the 5 chosen latitudinal bands
ra1 = rlranom(rlr_df, stn_df, wgt_df, yr_s, yr_e, 00, 25, .3)
ra2 = rlranom(rlr_df, stn_df, wgt_df, yr_s, yr_e, 25, 40, .8)
ra3 = rlranom(rlr_df, stn_df, wgt_df, yr_s, yr_e, 40, 50, .8)
ra4 = rlranom(rlr_df, stn_df, wgt_df, yr_s, yr_e, 50, 55, .9)
ra5 = rlranom(rlr_df, stn_df, wgt_df, yr_s, yr_e, 55, 60, .9)

# Plot the sample sizes for each band, if your wish
# plot(t,ra1$rlr_n,type="h")
# plot(t,ra2$rlr_n,type="h")
# plot(t,ra3$rlr_n,type="h")
# plot(t,ra4$rlr_n,type="h")
# plot(t,ra5$rlr_n,type="h")

# calculate the weighted average
wsum=sum(ra1$wgt,ra2$wgt,ra3$wgt,ra4$wgt,ra5$wgt)

wc1 = ra1$rlr_a * ra1$wgt/wsum
wc2 = ra2$rlr_a * ra2$wgt/wsum
wc3 = ra3$rlr_a * ra3$wgt/wsum
wc4 = ra4$rlr_a * ra4$wgt/wsum
wc5 = ra5$rlr_a * ra5$wgt/wsum

rlr_a = apply(rbind(wc1,wc2,wc3,wc4,wc5),2,sum)

# do the stats
yrsq= yrs^2

# quadratic of my model
rlr_mod=lm(rlr_a ~ yrsq + yrs)
summary(rlr_mod)

# setup the comparison, baseline is the entire period
rec_df = subset(recon_df,t %in% yrs)
rec_df$anom = rec_df$anom + mean(rlr_a) - mean(rec_df$anom) 

# quadratic of C&W model
rec_mod=lm(rec_df$anom ~ yrsq + yrs)
summary(rec_mod)

# plot results
windows()
# png("aj-gmsl.png", bg="transparent", width=450, height=450)

plot(yrs,rlr_a,type='n',main="AJ's Model vs C&W",xlab="year",ylab="MSL(mm)")
lines(yrs,rlr_a)
lines(yrs,predict(rlr_mod))
lines(yrs,rec_df$anom,col=2)
lines(yrs,predict(rec_mod),col=2)
grid()
legend("topleft",inset=0.01,legend=c("AJ's Model","C&W"),bg="white",
         text.col=1,col=c(1:2),lwd=1)
# dev.off()

# plot residuals (aj's model - C&W)
windows()
# png("aj-cw-gmsl-resid.png", bg="transparent", width=450, height=450)
plot(yrs,(rlr_a-rec_df$anom),type='l',
                main="Residuals AJ - C&W", xlab="year",ylab="MSL(mm)")
grid()
# dev.off()

# calculate slope 1930-2009
ty = 1930:2009

tc = rlr_a[yrs %in% ty]
summary(lm(tc~ty))

tc = rec_df$anom[yrs %in% ty]
summary(lm(tc~ty))

# calculate slope 1993-2009 (for comparison against altimetry)
ty = 1993:2009

tc = rlr_a[yrs %in% ty]
summary(lm(tc~ty))

tc = rec_df$anom[yrs %in% ty]
summary(lm(tc~ty))

# compare my model against C&W
summary(lm(rlr_a ~ rec_df$anom))

Year Anomaly(mm)
1930 0.0
1931 12.0
1932 16.1
1933 25.8
1934 3.1
1935 20.4
1936 25.3
1937 38.2
1938 37.7
1939 29.2
1940 31.2
1941 32.9
1942 42.3
1943 32.4
1944 31.6
1945 40.2
1946 50.5
1947 45.9
1948 66.3
1949 53.2
1950 61.1
1951 68.3
1952 75.8
1953 82.3
1954 74.3
1955 82.9
1956 76.8
1957 69.5
1958 80.3
1959 77.9
1960 86.1
1961 87.9
1962 82.9
1963 55.0
1964 67.0
1965 70.3
1966 73.5
1967 75.4
1968 65.0
1969 91.0
1970 96.4
1971 101.4
1972 91.8
1973 114.6
1974 116.1
1975 106.7
1976 84.9
1977 91.7
1978 111.6
1979 98.6
1980 93.7
1981 110.0
1982 83.2
1983 118.0
1984 122.1
1985 101.7
1986 96.7
1987 100.0
1988 122.1
1989 121.5
1990 110.1
1991 114.3
1992 127.5
1993 121.4
1994 126.1
1995 140.6
1996 153.9
1997 134.1
1998 152.1
1999 165.8
2000 159.2
2001 162.0
2002 147.5
2003 165.9
2004 164.7
2005 168.1
2006 162.7
2007 159.2
2008 186.0
2009 178.8

ċ
Aj J,
Sep 21, 2013, 8:39 AM
ċ
rlr_annual.zip
(1235k)
Aj J,
Sep 21, 2013, 8:51 AM
ċ
Aj J,
Sep 21, 2013, 8:39 AM
Comments