############## RClimate Script: GISS Trends by PDO, AMO, ENSO Phases   #####################
## Script stored on http://chartsgraphs.wordpress.com account for Users to source()       ##                                                                                       ## 
## Download and process     Consolidated monthly file                                     ##
## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ##
##                   http:chartsgraphs.wordpress.com    4/27/10                           ##
##                   Updated 5/28/2010 to correct lm of AMO                               ##
############################################################################################
## Read consolidated data file
  link_in <- "http://processtrends.com/Files/pdo_enso_amo_comb_index_data.csv" 
  df <- read.csv(link_in)
  par(ask=T)

## Figure 1: Boxplot of PAE vs GISS LOTA-C
 par(mfrow=c(1,1))
 par(oma=c(2,1,2,1)); par(mar=c(4,3,2,1))
 pae_lab <- c("1: PDO:-\nAMO:-\nNino34:-", 
              "2: PDO:-\nAMO:-\nNino34:+",
              "3: PDO:-\nAMO:+\nNino34:-",
              "4: PDO:-\nAMO:+\nNino34:+",
              "5: PDO:+\nAMO:-\nNino34:-", 
              "6: PDO:+\nAMO:-\nNino34:+",
              "7: PDO:+\nAMO:+\nNino34:-",
              "8: PDO:+\nAMO:+\nNino34:+")
 boxplot(df$giss ~ df$pae_score, axes=F, xaxs="i", ylim=c(-1,1),
   yaxs="i", pch=16, ylab="GISS LOT Anomaly - C",
   main = "Boxplots of GISS LOTA-C by PAE Codes\n1900-2010")
 axis(1, at=c(1:8), tick=F,labels=pae_lab, cex.axis=0.65)
 axis(2)
 abline(h=-1)
 mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,1, adj = 0, cex = 0.7, outer=TRUE)
 mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.7, outer=TRUE)

####################################################################################
## Figure 2: Trend charts of PAE vs time
   pae<-matrix(data = NA, nrow = 8, ncol = 6, byrow = FALSE, dimnames = NULL) 
   score_lab <- c("---", "--+", "-+-", "-++", "+--", "+-+", "++-", "+++")
###########################################################################
  par(mfrow=c(4,2)); par(las=1)
  par(mar=c(1.5,4,1,1)); par(oma=c(2.5,2,2,2))
  for(score in 1:8) { 
    nnn <- subset(df, df$pae_score == score)
    pae[score,1] <- nrow(nnn)
    pae[score,2] <- min(nnn$giss, na.rm=T)
    pae[score, 3] <- signif(mean(nnn$giss, na.rm=T),3)
    pae[score,4] <- signif(median(nnn$giss, na.rm=T),3)
    pae[score,5] <- signif(max(nnn$giss, na.rm=T),3)
    plot(nnn$yr_frac, nnn$giss, type = "p", pch=20, col = "grey", ylim=c(-0.5,1),
      xlab="", ylab = "",xlim=c(1900, 2011)  )
    reg <- lm(nnn$giss~nnn$yr_frac)
    pae[score,6] <- signif(coef(reg)[2],3)
    abline(reg, col="red")
    score_des <- paste(score, ". ", score_lab[score], sep="")
    text(1905,0.45, paste("no obs: ", pae[score,1], sep=""), adj=0,cex=0.85)
    text(1905, 0.6, paste("mean: ",pae[score,3], sep=""),adj=0,cex=0.85)
    text(1905, 0.75, paste("slope: ",signif(coef(reg)[2],2), sep=""),adj=0,cex=0.85)  
    text(1905,0.95, score_des, adj=0,cex=1.1)  
  }
  mtext("GISS LOTA - C", side = 2, line =0, cex = 1,las=0, outer = T, adj = 0.5)## Add script name to each plot
  mtext("GISS LOTA Trend by Combined PDO, AMO, ENSO Phase",side=3, line=0.5, font=2,cex=1.1, outer=T, adj=0.5,las=1)
  mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,1, adj = 0, cex = 0.7, outer=TRUE)
  mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.7, outer=TRUE)
pae
################################################

# Figure 3: GISS Regressions by Year, PDO, AMO, ENSO Phases
 ## Updated AMO regression based on reader comment ######
  r_pdo <- lm(giss~pdo, data=df)
  r_amo <- lm(giss~amo, data=df)
  r_nino <- lm(giss~nino34, data=df)
  r_pan <- lm(giss~ pdo + amo+ nino34, data = df)
  r_yr <- lm(giss~yr_frac , data=df)
  r_yr_pan <- lm(giss~yr_frac + pdo + amo+ nino34, data = df)

 sum_r_pdo <- summary(r_pdo)
 sum_r_amo <- summary(r_amo)
 sum_r_nino <- summary(r_nino)
 sum_r_pan <- summary(r_pan)
 sum_r_yr <- summary(r_yr)
 sum_r_yr_pan <- summary(r_yr_pan)

# make prediction series for plotting
 pred_pdo <- coef(r_pdo)[1] + df$pdo * coef(r_pdo)[2]
 pred_amo <- coef(r_amo)[1] + df$amo * coef(r_amo)[2]
 pred_nino <- coef(r_nino)[1] +df$nino34 * coef(r_nino)[2]
 pred_pan <- coef(r_pan)[1] + df$pdo*coef(r_pan)[2] +df$amo *coef(r_pan)[3] + df$nino34 * coef(r_pan)[4]
 pred_yr <- coef(r_yr)[1] + df$yr_frac * coef(r_yr)[2]
 pred_yr_pan <- coef(r_yr_pan)[1] + df$yr_frac * coef(r_yr_pan)[2]+ df$pdo*coef(r_yr_pan)[3] +df$amo *coef(r_yr_pan)[4]+ df$nino34 * coef(r_yr_pan)[5]

 df <- data.frame(df, pred_pdo, pred_amo, pred_nino, pred_pan,pred_yr, pred_yr_pan)
  
# Plot Function ############################################################
  p_func <- function() { plot(df$yr_frac, df$giss, type = "l", col = "grey",
    xlim = c(1900, 2010), xlab="", ylab="GISS LOTA- C", las=1)
    }
############################################################################
## Produce 6 plots 
  par(oma=c(2,2.5,2.5,1)); par(mar=c(2,3,1,0))
  par(mfrow=c(3,2))
  p_func()
   text(1905, 0.8, "PDO", adj=0, cex=1, font=2)
   points(df$yr_frac, df$pred_pdo, type = "l", col = "green")
  p_func()
   text(1905, 0.8, "AMO", adj=0, cex=1, font=2)
   points(df$yr_frac, df$pred_amo, type = "l", col = "green")
  p_func()
   text(1905, 0.8, "NINO34", adj=0, cex=1, font=2)
   points(df$yr_frac, df$pred_nino, type = "l", col = "green")
  p_func()
   text(1905, 0.8, "PDO + AMO + NINO34", adj=0, cex=1, font=2)
   points(df$yr_frac, df$pred_amo, type = "l", col = "green")
  p_func()
   text(1905, 0.8, "Year", adj=0, cex=1, font=2)
   points(df$yr_frac, df$pred_yr, type = "l", col = "green")
  p_func()
   text(1905, 0.8, "Year + PDO + ADO + ENSO", adj=0, cex=1, font=2)
   points(df$yr_frac, df$pred_yr_pan, type = "l", col = "orange")
 ## Plot annotation
   mtext("GISS LOTA-C", side =2, line =1, adj=0.5, outer=T, las=0)
   mtext("GISS LOTA-C Regressions Using Year, PDO, ADO, ENSO Phases", side=3, line=1, font=2,outer=T)
   mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,0.5, adj = 0, cex = 0.7, outer=TRUE)
   mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 0.5, adj = 1, cex = 0.7, outer=TRUE)