Global mean sea level

dataset - GMSL time series (NASA)

In [1]:
library(zoo)
library(waveslim)

options(repr.plot.width=8, repr.plot.height=4)

Attaching package: ‘zoo’

The following objects are masked from ‘package:base’:

    as.Date, as.Date.numeric


waveslim: Wavelet Method for 1/2/3D Signals (version = 1.7.5)
In [11]:
dat=read.table(file="GMSL_TPJAOS_V4.2_199209_201711.txt",header=FALSE)

# year+fraction of year (mid-cycle) 
time.msl=dat$V3

# smoothed (60-day Gaussian type filter) GMSL (GIA applied) variation (mm); annual and semi-annual signal removed
msl.mean=dat$V12
plot(zoo(msl.mean,time.msl), xlab="",ylab="MSL (mm)")

# linear detrending
msl.detrended=residuals(lm(msl.mean~time.msl))
plot(zoo(msl.detrended,time.msl), xlab="",ylab="MSL - detrended (mm)")

ENSO events - MEI index

In [94]:
mei=t(read.table(file="mei.data",header=FALSE)[,2:13])

library(tidyr)
mei.ts=ts(gather(as.data.frame(mei),V2:V13)$value,start=1993, frequency=12)

plot(time(mei.ts)[which(mei.ts>0)],mei.ts[which(mei.ts>0)],type="h",xlab="",ylab="MEI enso index",col=2,ylim=c(-2,3))
lines(time(mei.ts)[which(mei.ts<0)],mei.ts[which(mei.ts<0)],type="h",xlab="",ylab="MEI enso index",col=4,ylim=c(-2,3))
In [12]:
plot(zoo(msl.detrended,time.msl),type="h",xlab="",ylab="MSL (mm)")
In [14]:
acf(msl.detrended, main="GMSL - detrended",lag.max=100)
pacf(msl.detrended, main="GMSL - detrended",lag.max=100)
In [31]:
library(forecast)

auto.arima(msl.detrended)

Series: msl.detrended 
ARIMA(4,1,2) 

Coefficients:
         ar1      ar2      ar3     ar4      ma1     ma2
      1.2175  -0.5316  -0.0313  0.0158  -0.0193  0.4298
s.e.  0.0999   0.1522   0.1368  0.0707   0.0945  0.0903

sigma^2 estimated as 0.02844:  log likelihood=331.71
AIC=-649.42   AICc=-649.3   BIC=-615.7
In [18]:
arima(msl.detrended,order=c(1,0,0))

Call:
arima(x = msl.detrended, order = c(1, 0, 0))

Coefficients:
         ar1  intercept
      0.9945     0.6819
s.e.  0.0032     2.0298

sigma^2 estimated as 0.159:  log likelihood = -459.17,  aic = 924.33

Normality assessment

In [15]:
par(mfrow=c(1,2))

hist(msl.detrended,main="", xlab="mm")

qqnorm(msl.detrended)
qqline(msl.detrended)

Wavelet spectrum

In [40]:
msl.modwt <- modwt(msl.detrended, "haar", 6)

msl.modwt.bw <- brick.wall(msl.modwt, "haar")

msl.modwt.var <- wave.variance(msl.modwt.bw)

matplot(2^(0:5), msl.modwt.var[-7,], type="b", log="xy",
        xaxt="n", ylim=c(.001, 6), pch="*LU", lty=1, col=c(1,4,4),
        xlab="wavelet Scale", ylab="", main="Wavelet spectrum (Haar wavelet filter)")
axis(side=1, at=2^(0:5))
In [46]:
msl.modwt <- modwt(msl.detrended, "la8", 6)

msl.modwt.bw <- brick.wall(msl.modwt, "la8")

msl.modwt.var <- wave.variance(msl.modwt.bw)

matplot(2^(0:5), msl.modwt.var[-7,], type="b", log="xy",
        xaxt="n", ylim=c(.001, 6), pch="*LU", lty=1, col=c(1,4,4),
        xlab="wavelet Scale", ylab="",main="Wavelet spectrum (LA8 filter)")
axis(side=1, at=2^(0:5))