week4
Sara Cendejas-Zarelli
April 20, 2016
Part 1: Filtering:
co2 data set filtering
plot(co2)
The data is composed of monthly co2 data (ppm) and extends from 1959 to 1998
Creating a time (annual) component for the data shows the following:
tm <- time(co2)Using a basic filter at 9 and 12 months shows the following patterns:
ma9 <- filter(x=co2, filter=rep(x=1/9,times=9), sides=2)
ma12 <- filter(x=co2,filter=rep(x=1/12,times=12), sides=2)
plot(co2,col="grey2")
lines(ma9,col="red",lwd=2)
lines(ma12,col="green",lwd=2)
The better filter is the 12 month…this makes sense because it is smoothing the data in to a more annual mean, and thus becomes a more or less even line throught the data.
Another way to look at this is to aggregate the monthly data into annual data:
co2.ag <- aggregate(co2,FUN=mean)
plot(co2.ag,ylab=expression(CO[2]~(ppm)))
We can see the line is identical to the green m12 line in the previous figure.
Sea ice data set
library("repmis")githubURL <- "https://github.com/AndyBunn/TeachingData/raw/master/climate_Time_Series_Extravaganza.Rdata"
source_data(githubURL)## Downloading data from: https://github.com/AndyBunn/TeachingData/raw/master/climate_Time_Series_Extravaganza.Rdata## SHA-1 hash of the downloaded data file is:
## 5f5b467b42520fdfbb85cfc8a1fbcc45197e0e56## [1] "loti.ts" "loti.zoo" "ice.ts" "ice.zoo"
## [5] "sealevel.ts" "sealevel.zoo" "sunspots.ts" "sunspots.zoo"
## [9] "enso.ts" "enso.zoo" "amo.ts" "amo.zoo"
## [13] "pdo.ts" "pdo.zoo" "ohc.ts" "ohc.zoo"
## [17] "co2.ts" "co2.zoo"plot(ice.ts)
Aggregating by month
ice.ag <- aggregate(ice.ts,FUN=mean)
plot(ice.ag,ylab=expression(ice.ts))
tm.ice<-time(ice.ag)Moving average filter
ma3 <- filter(x=ice.ag, filter=rep(x=1/3,times=3), sides=2)
ma8 <- filter(x=ice.ag, filter=rep(x=1/8,times=8), sides=2)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(ice.ag,type="l",xlab="Year",ylab="Sea ICe",col="grey2")
abline(h=1)
lines(ma3, col = "red", lwd = 2)
lines(ma8, col = "green", lwd = 2)
axis(3);axis(4)
The ends of the filters lack critical data!
Lowess filtering
f310 <- 310/tm.ice
f310.lo <- lowess(x = tm.ice, y = ice.ag, f = f310)
f410 <- 410/tm.ice
f410.lo <- lowess(x = tm.ice, y = ice.ag, f = f410)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(ice.ag,xlab="Year",ylab="Sea Ice",col="grey2")
lines(f310.lo, col = "red", lwd = 2)
lines(f410.lo, col = "green", lwd = 2)
It took some very large numbers to smooth the data! I’m not entirely sure why…?
Part 2
Ocean temperature variablility
library(repmis)githubURL <- "https://github.com/AndyBunn/TeachingData/raw/master/climate_Time_Series_Extravaganza.Rdata"
source_data(githubURL) # will load all the objects.## Downloading data from: https://github.com/AndyBunn/TeachingData/raw/master/climate_Time_Series_Extravaganza.Rdata## SHA-1 hash of the downloaded data file is:
## 5f5b467b42520fdfbb85cfc8a1fbcc45197e0e56## [1] "loti.ts" "loti.zoo" "ice.ts" "ice.zoo"
## [5] "sealevel.ts" "sealevel.zoo" "sunspots.ts" "sunspots.zoo"
## [9] "enso.ts" "enso.zoo" "amo.ts" "amo.zoo"
## [13] "pdo.ts" "pdo.zoo" "ohc.ts" "ohc.zoo"
## [17] "co2.ts" "co2.zoo"oceans3 <- cbind(ENSO=enso.ts,PDO=pdo.ts,AMO=amo.ts)
par(tcl=0.5,mar=rep(2.5,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(oceans3,main="Three Ocean Patterns")
PDO
plot(pdo.ts)
Aggregating PDO by month:
pdo.ag <- aggregate(pdo.ts,FUN=mean)
plot(pdo.ag,ylab=expression(pdo.ts))
tm.pdo <- time(pdo.ag)Moving average filter
ma5 <- filter(x=pdo.ag, filter=rep(x=1/5,times=5), sides=2)
ma25 <- filter(x=pdo.ag, filter=rep(x=1/25,times=25), sides=2)
ma40 <- filter(x=pdo.ag, filter=rep(x=1/40,times=40), sides=2)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(pdo.ag,type="l",xlab="Year",ylab="PDO",col="grey2")
abline(h=1)
lines(ma5, col = "red", lwd = 2)
lines(ma25, col = "green", lwd = 2)
lines(ma40, col = "blue", lwd = 2)
axis(3);axis(4)
The ends of the filters lack critical data!
Lowess filtering
f200 <- 200/tm.pdo
f200.lo <- lowess(x = tm.pdo, y = pdo.ag, f = f200)
f300 <- 300/tm.pdo
f300.lo <- lowess(x = tm.pdo, y = pdo.ag, f = f300)
f128 <- 128/tm.pdo
f128.lo <- lowess(x = tm.pdo, y = pdo.ag, f = f128)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(pdo.ag,xlab="Year",ylab="PDO",col="grey2")
abline(h=1)
lines(f200.lo, col = "red", lwd = 2)
lines(f300.lo, col = "green", lwd = 2)
lines(f128.lo, col = "blue", lwd = 2)
axis(3);axis(4)
A the Lowess method, we can see that the 128 year filter does a better job of capturing the upward trend that occurs at the right end of the data than either the 200 or 300 filters. The 300 year filter seems to do the worst in capturing that upward trend, but overall it smoothes the data the best.
Spectral Analysis
n.vec <- pdo.ag
pdo.wav5 <- 0.5 * sin(2 * pi / 5 * n.vec)
pdo.wav50 <- 0.75 * sin(2 * pi / 50 * n.vec)
pdo.wav100 <- sin(2 * pi / 100 * n.vec)
pdo.wav250 <- sin(2 * pi / 250 * n.vec)
plot(pdo.wav250)
lines(pdo.wav5,col="blue")
lines(pdo.wav50,col="green")
lines(pdo.wav100,col="red")
lines(pdo.wav250,col="purple")
I don’t think I plotted this right, but it looks interesting!
noise <- rnorm(pdo.ag)
wav <- pdo.wav250 + pdo.wav100 + pdo.wav50 + pdo.wav5 + noise
plot(wav,type="l")
spectrum(wav,log="no")
I need to play around with the size of each wav, but I can see a some pretty large peaks at at ~0.39 and again at ~0.5. In order to get the frequency (1/0.39=2.6) So there is a larger occurrence of data at 2.6 years and again at 2 years. This suggests that there may be a two year cycle in sea PDO
AMO
plot(amo.ts)
Aggregating AMO by month:
amo.ag <- aggregate(amo.ts,FUN=mean)
plot(amo.ag,ylab=expression(amo.ts))
tm.amo<-time(amo.ag)Moving average filter
ma5 <- filter(x=amo.ag, filter=rep(x=1/5,times=5), sides=2)
ma25 <- filter(x=amo.ag, filter=rep(x=1/25,times=25), sides=2)
ma40 <- filter(x=amo.ag, filter=rep(x=1/40,times=40), sides=2)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(amo.ag,type="l",xlab="Year",ylab="AMO",col="grey2")
abline(h=1)
lines(ma5, col = "red", lwd = 2)
lines(ma25, col = "green", lwd = 2)
lines(ma40, col = "blue", lwd = 2)
axis(3);axis(4)
The shortest filter (ma5) extends the furthest outward in the data set than either the 25 or 40 year filters. This makes sense because it only needs 5 years on either side to draw its line. You can see the 25 year filter (green) is a little longer the the 40 year filter (blue) for this same reason.
Lowess filtering
f100 <- 100/tm.amo
f100.lo <- lowess(x = tm.amo, y = amo.ag, f = f100)
f200 <- 200/tm.amo
f200.lo <- lowess(x = tm.amo, y = amo.ag, f = f200)
f300 <- 300/tm.amo
f300.lo <- lowess(x = tm.amo, y = amo.ag, f = f300)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(amo.ag,xlab="Year",ylab="AMO",col="grey2")
abline(h=1)
lines(f100.lo, col = "red", lwd = 2)
lines(f200.lo, col = "green", lwd = 2)
lines(f300.lo, col = "blue", lwd = 2)
axis(3);axis(4)
A filter of 300 smoothes the line the best, but does not do the greatest job of finding the mean of the line at the ends.
Spectral Analysis
n.vec1 <- amo.ag
amo.wav5 <- 0.5 * sin(2 * pi / 5 * n.vec1)
amo.wav50 <- 0.75 * sin(2 * pi / 50 * n.vec1)
amo.wav100 <- sin(2 * pi / 100 * n.vec1)
amo.wav250 <- sin(2 * pi / 250 * n.vec1)
plot(amo.wav250)
lines(amo.wav5,col="blue")
lines(amo.wav50,col="green")
lines(amo.wav100,col="red")
lines(amo.wav250,col="purple")
I don’t think I plotted this right, but it looks interesting!
noise1 <- rnorm(amo.ag)
wav1 <- amo.wav250 + amo.wav100 + amo.wav50 + amo.wav5 + noise1
plot(wav1,type="l")
spectrum(wav1,log="no")
I can see a large peak at at ~0.01 and again at ~3.8 In order to get the frequency (1/0.01=100) So there is a larger occurrence of data at 100 and 2.6 years. If I were confident in my plot, I would interpret this to mean there is a 100 year cycle in the AMO.
ENSO
plot(enso.ts)
Aggregating by month:
enso.ag <- aggregate(enso.ts,FUN=mean)
plot(enso.ag,ylab=expression(enso.ts))
tm.enso<-time(enso.ag)Moving average filter
ma5 <- filter(x=enso.ag, filter=rep(x=1/5,times=5), sides=2)
ma25 <- filter(x=enso.ag, filter=rep(x=1/25,times=25), sides=2)
ma40 <- filter(x=enso.ag, filter=rep(x=1/40,times=40), sides=2)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.0,0.1,0),xaxs="i")
plot(enso.ag,type="l",xlab="Year",ylab="ENSO",col="grey2")
abline(h=1)
lines(ma5, col = "red", lwd = 2)
lines(ma25, col = "green", lwd = 2)
lines(ma40, col = "blue", lwd = 2)
axis(3);axis(4)
Again, the same trend occurs in the lack of data at the ends
Lowess filtering
f150 <- 150/tm.enso
f150.lo <- lowess(x = tm.enso, y = enso.ag, f = f150)
f200 <- 200/tm.enso
f200.lo <- lowess(x = tm.enso, y = enso.ag, f = f200)
f300 <- 300/tm.enso
f300.lo <- lowess(x = tm.enso, y = enso.ag, f = f300)
par(tcl=0.5,mar=rep(3,4),mgp=c(1.1,0.1,0),xaxs="i")
plot(enso.ag,xlab="Year",ylab="ENSO",col="grey2")
abline(h=1)
lines(f150.lo, col = "red", lwd = 2)
lines(f200.lo, col = "green", lwd = 2)
lines(f300.lo, col = "blue", lwd = 2)
axis(3);axis(4)
Once again, the largest (300) spline smoothes the data the best…
Spectral Analysis
n.vec2 <- enso.ag
enso.wav5 <- 0.5 * sin(2 * pi / 5 * n.vec2)
enso.wav50 <- 0.75 * sin(2 * pi / 50 * n.vec2)
enso.wav100 <- sin(2 * pi / 100 * n.vec2)
enso.wav250 <- sin(2 * pi / 250 * n.vec2)
plot(enso.wav250)
lines(enso.wav5,col="blue")
lines(enso.wav50,col="green")
lines(enso.wav100,col="red")
lines(enso.wav250,col="purple")
I don’t think I plotted this right, but it looks interesting!
noise2<- rnorm(enso.ag)
wav2 <- enso.wav250 + enso.wav100 + enso.wav50 + enso.wav5 + noise2
plot(wav2,type="l")
spectrum(wav2,log="no")
I can see some serious peaks at ~0.23 and ~3.5… In order to get the frequency (1/0.23=4.3) So there is a larger occurrence of data at 4.3 and 2.9 years. I can see a large peak at at ~0.01 and again at ~3.8
If I were confident in my plot, I would interpret this to mean there is a 3 and 4 year cycle in ENSO.