tufteTicks <- function(data, side=1) {
	# Function to calculate tick locations for 
	# a "range frame" graph, a la Edward Tufte's
	# "Visual Display of Quantitative Information."
	# Works analagously to axTicks()
	#
	# Arguments:
	#   data: Values to calculate ticks and range from.
	#   side: Side of the graph to place the axis on.
	#         Defaults to 1 (bottom).
	ax <- axTicks(side)
	n <- length(ax)
	lo <- min(data)
	hi <- max(data)
	if (ax[1] < lo)
		ax[1] <- lo
	else if (ax[1] > lo)
		ax <- c(lo, ax)
	if (ax[n] < hi)
		ax <- c(ax, hi)
	else if (ax[n] > hi)
		ax[n] <- hi
	return(round(ax, 2))
}

#################################################
# Script
#################################################

# Read in data
rain <- read.csv('qlnd_rain.csv')
soi <- read.csv('soi.csv')
# Clip SOI data to after 1900
soi <- soi[soi$year >= min(rain$year), ]

# Add decimal years
rain$year.dec <- rain$year + rain$month * 1/12
soi$year.dec <- soi$year + soi$month * 1/12 

# Find monthly means, remove them from rainfall
mod.month <- lm(rain ~ as.factor(month), data=rain)
rain$rain.avg <- mod.month$fit
rain$rain.ds <- rain$rain - rain$rain.avg
monthly.avg  <- predict(mod.month, newdata=list(month=1:12))
# calculate monthly standard deviations
monthly.sd <- rep(0, 12)
for (i in 1:12) {
	monthly.sd[i] <- sd(rain$rain[rain$month == i])
}

mod.1 <- lm(rain$rain.ds ~ soi$soi)
summary(mod.1)
anova(mod.1)

soi.rain.ccf <- ccf(soi$soi, rain$rain.ds, lag.max=36, plot=FALSE)

#################################################
# Plots
#################################################

# rainfall, 1970-
plot(rain ~ year.dec, data=rain[rain$year > 1970, ],
	ty='l', bty='n', yaxt='n', xlab='', ylab='Rainfall (mm)',
	main='Monthly Queensland Rainfall, 1970-2010')
axis(1, at=1970:2011, labels=FALSE, tcl=-0.3)
axis(2, axis(2, at=round(tufteTicks(rain$rain, 2), 0)))

# Climatology (1900-)
# (this one isn't in the blog post)
plot(monthly.avg, bty='n', ylim=c(0, 200), xaxt='n', xlab='Month', ylab='Rainfall (mm)')
abline(h=seq(0, 200, 25), lty=1, lwd=0.5, col=rgb(0.9, 0.9, 0.9))
axis(1, at=1:12)
upper <- monthly.avg + monthly.sd
lower <- monthly.avg - monthly.sd
arrows(1:12, monthly.avg, 1:12, upper, angle=90, length=0.05, col='black')
arrows(1:12, monthly.avg, 1:12, lower, angle=90, length=0.05, col='black')
lines(monthly.avg, ty='b', pch=16)
rain.10 <- rain$rain[rain$year == 2010]
lines(rain.10, ty='b', pch=15, col='red')

# Climatology (1900-)
# Plot each year as an individual line
plot(rain ~ month, data=rain[rain$year == 1900, ],
	ty='l', bty='n', col=rgb(0,0,0,0.15), ylim=c(0,475),
	xaxt='n', yaxt='n', xlab='Month', ylab='Rainfall (mm)')
axis(1, at=1:12)
axis(2, at=round(tufteTicks(rain$rain, 2), 0))
for (y in 1901:2010) {
	lines(rain ~ month, data=rain[rain$year == y, ], ty='l', col=rgb(0,0,0,0.15))
}
lines(rain ~ month, data=rain[rain$year == 1973, ][7:12, ],
	ty='l', col='red')
lines(rain ~ month, data=rain[rain$year == 1974, ][1:7, ],
	ty='l', col='red')
lines(rain ~ month, data=rain[rain$year == 2010, ], ty='l', col='blue')
lines(monthly.avg ~ month, data=rain[rain$year == 2010, ],
	ty='l', col='black')
points(monthly.avg ~ month, data=rain[rain$year == 2010, ],
	col='black', pch=15)
legend('topright', c('1973-74', '2010', 'Average'), lty=1, col=c('red', 'blue', 'black'), bty='n')

# rainfall - climatology (1970-)
# deseasonalized rainfall
plot(rain.ds ~ year.dec, data=rain[rain$year > 1970, ], ty='l')
axis(1, at=1970:2011, labels=FALSE, tcl=-0.3)

# rainfall vs. soi
plot(rain$rain.ds ~ soi$soi, cex=0.5, xlab='SOI (mPa)',
	ylab='Rainfall (mm)', col=rgb(0,0,0,0.5), bty='n',
	xlim=c(-45, 40), ylim=c(-100, 400))
abline(mod.1, col='red')
text(-40, 300, pos=4, expression(R^2==0.1292))
text(-40, 260, pos=4, 'p << 0.001')

# Cross-correlation
plot(soi.rain.ccf, xaxt='n', yaxt='n', bty='n', 
	xlab='Lag (months)', ylab='Cross-Correlation', 
	main='Cross-Correlation of SOI with Queensland Rainfall Anomaly')
axis(side=1, at=seq(-36, 36, by=12))
axis(side=2, at=tufteTicks(soi.rain.ccf$acf, 2))