##                         RClimate_Arctic_osc_monthly.R                                   ###################
##  Download & Plot Monthly AO Index: 1950 to Present                                      ###################
##  http://www.cpc.noaa.gov/products/precip/CWlink/daily_ao_index/monthly.ao.index.b50.current.ascii        ##
##  D Kelly O'DAy http://chartsgraphs.wordpress.com                                                         ##
##  Orig: 9/23/10                                                                                           ##
##  Updated: 1/14/11 Legend note on filter corrected to 60 month moving avg; reader skrafner comment        ##
##############################################################################################################
  link <- "http://www.cpc.noaa.gov/products/precip/CWlink/daily_ao_index/monthly.ao.index.b50.current.ascii "
  AO_df <- read.table(link, skip = 0, na.strings = c("-0.99900E+34"), header=F)
  names(AO_df)<- c("yr", "mo", "AO")
  yr_frac <- as.numeric(AO_df$yr) + ((AO_df$mo-0.5)/12)
  AO_df <- data.frame(yr_frac, AO_df)
  AO_df <- subset(AO_df, AO_df$AO != "NA")
# Create pos & neg series for color coding AOs
  AO_p <- subset(AO_df, AO >= 0)
  AO_n <- subset(AO_df, AO <0)
# Calc 60 month moving avg 
 AO_ma <- filter(AO_df$AO, rep(1/60,60),sides=1)
# Specify min & max years to provide zoom capability
 x_min = 1950
 x_max = 2012

## Information on Last data point in series for highlighting & legend 
 mon_name <- c("January", "February", "March", "April", "May", "June", "July", "August", "September", "October", "November", "December")
 c <- nrow(AO_df)
 last_AO <- signif(AO_df[c,4],3)
 last_yr_mon <- paste(mon_name[AO_df[c,3]], ", ", AO_df[c,2], sep = "") 
 last_val_note <- paste(last_yr_mon, " AOI @ ", last_AO, sep="")
 
##    Plot 
 par(xaxs="i"); par(yaxs="i")
 par(oma=c(3,1,1,1)); par(mar=c(2,4,3,1))
 title <- paste("Monthly Arctic Oscillation Index\n January, 1950 to ", last_yr_mon, sep="")
 plot(AO_df$yr_frac, AO_df$AO, type = "n", xlab="", ylab="Arctic Oscillation Index",
   main = title,  col = "red", ylim = c(-5, 4),
   xlim=c(x_min,x_max),las=1)
 #grid(ny=NA,col = "lightgrey", lty=1)
 points(AO_p$yr_frac, AO_p$AO, col = "red", type="h")
 points(AO_n$yr_frac, AO_n$AO, col = "blue", type = "h")
 abline(h=0, col="grey")
 points(AO_df$yr_frac, AO_ma, col="black", type="l",lwd=2.5)
 points(AO_df$yr_frac[c], AO_df$AO[c], type = "p",pch=16, col = "yellow")

 legend(1950, -3.3,  c("Positve AOI", "Negative AOI", "60-Month Moving Average AOI", last_val_note ), col = c("red", "blue","black", "yellow"),
       text.col = "black", lty = c(1,1,1,0),pch=c(0,0,0,16),pt.cex=c(0,0,0,1),
       merge = F, bg = "white", bty="n", cex = .7)
# Generate and add bottom footer KOD , source, System date notes
       source_note <- paste("Data Source: ", link)
       mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,1, adj = 0, cex = 0.8, outer=TRUE)
       mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.8, outer=TRUE)
       mtext(source_note, 1,0,adj=0.5, cex=0.8, outer=T)
#############################################################################################