setwd("~/Western")
nino<-readRDS("nino.rds")
ts.plot(nino)
tsp(nino)
## [1] 1950.000 2016.917 12.000
summary(nino)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.24300 -0.59420 0.05350 0.07715 0.69620 3.04400
require(dplR)
## Loading required package: dplR
## Warning: package 'dplR' was built under R version 3.3.3
##Multivariate El Nino Index anomalies from the shiny app
#I decided on 10 year, 25 year, and 40 year moving averages
ma10 <- filter(x=nino, filter=rep(x=1/10,times=10), sides=2)
ma25 <- filter(x=nino, filter=rep(x=1/25,times=25), sides=2)
ma40 <- filter(x=nino, 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(nino,type="l",xlab="Year",ylab="RWI",col="grey")
abline(h=1)
lines(ma10, col = "red", lwd = 2)
lines(ma25, col = "green", lwd = 2)
lines(ma40, col = "blue", lwd = 2)
axis(3);axis(4)
##The moving averages revealed interesting cycles of 10 years, which just so happens to correspond to the Pacfic Decadal Oscillation. This may help explain factors determining El Nino years.
###Lowess smoother using a fraction of the data points to calculate polynomial fits
f0.25.lo<-lowess(nino, f=1/4)
f0.15.lo<-lowess(nino,f=1/6)
f0.10.lo<-lowess(nino,f=1/10)
plot(nino)
lines(f0.25.lo,col="red",lwd=2)
lines(f0.15.lo,col="green",lwd=2)
lines(f0.10.lo,col="blue",lwd=2)
## The lowess smoother using 10% of the data points showed the most clear decadal trend
##First used the least smooth lowess filter to remove low frequency cycles
rwi<-nino/f0.10.lo$y
par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.25,0.25,0),tcl=0.5)
plot(nino,type="l",ylab="mei",xlab="",col="grey",axes=F)
lines(f0.10.lo$y,lwd=2)
axis(2);axis(3);axis(4);box()
par(mar=c(3,3,0.1,3))
plot(rwi,type="l",ylab="RWI",xlab="time",col="black",axes=F)
abline(h=1)
axis(2);axis(1);axis(4);box()
###
##Next used the 25% smoother with the flatter line to remove low frequency cycles
rwi2<-nino/f0.25.lo$y
par(mfcol=c(2,1),mar=c(0.1,3,3,3),mgp=c(1.25,0.25,0),tcl=0.5)
plot(nino,type="l",ylab="mei",xlab="",col="grey",axes=F)
lines(f0.25.lo$y,lwd=2)
axis(2);axis(3);axis(4);box()
par(mar=c(3,3,0.1,3))
plot(rwi2,type="l",ylab="RWI2",xlab="time",col="black",axes=F)
abline(h=1)
axis(2);axis(1);axis(4);box()
## This data was taken monthly but it is useful to look at annual averages so I will use filtering to take an annual average
nino.filt <- filter(x = nino,filter = rep(x=1/12,times=12),sides=2)
head(cbind(nino,nino.filt))
## nino nino.filt
## [1,] -1.043 NA
## [2,] -1.162 NA
## [3,] -1.312 NA
## [4,] -1.104 NA
## [5,] -1.433 NA
## [6,] -1.391 -1.10725
annual.index<-seq(from=6,to=length(nino.filt),by=12)
nino.annual<-nino.filt[annual.index]
length(nino.annual)
## [1] 67
par(mar=rep(4,4))
plot(nino.annual,type="b",xlab="Year")
##another way to aggregate data from monthly into annual. As you can see there are 10 points for 10 years, one point every year! The graph looks identical to the one created using the filter function, so now I actually know how the black magic of the aggregate function works!
ag.nino<-aggregate(nino)
plot(ag.nino,type="b")
##Next we will switch gears and look at some data in the frequency domain using spectral analysis to find which frequencies dominate the summer insolation set, in the hopes of finding the periocidy of ice ages
jul65N <- readRDS("jul65N.rds")
plot(jul65N$kya,jul65N$W.per.m2,type="l",
xlab="Thousands of years ago",
ylab=expression(W/m^2),main="July insolation at 65N")
##spectral analysis of summer insolation, without a log axis
iceage.spec <- spectrum(jul65N, span=5, plot= FALSE)
plot(iceage.spec,log="no", type="h", xlab="Frequency (cycles/thousands of years ago)", ylab="Summer insolation W/m-squared",main="Annual Insolation")
##It appears that the strongest frequencies are at 0.05-0.06 which correspond with a period of 16,000-20,000-years cycles. This is consistent with findings of interglacial periods lasting 10,000-30,000 years (Kukla, 2004)
##There is also a dominant periocidy at 40,000 years so there may be two trends happening here at two different levels, where the 40,000 cycle correspond to the Earth's tilt and the Milankovitch cycle, which switches every 41,000 years.
##Kukla, George. 2004. Saalian supercycle, Mindel/Riss interglacial, and Milankovitch's dating. Quaternary Science Reviews. 24(15): 1573-1583.