#################################################################################################
## RClimate script to calculate common baseline LOTA for 5 series;                              # 
## Plot trendslines and scatter plots                                                           # 
## Original: 1/5/11                                                                             # 
## D Kelly O'Day      http://chartsgraphs.wordpress.com                                          # 
#################################################################################################
  graphics.off()
  save_dir<- "C:\\R_Home\\Charts & Graphs Blog\\RClimateTools\\Global_LOTA_consolidated_data\\"
## Establish baseline period - user can set baseline to any desired period
  baseline <- c(1979, 2008)
#################################################################################################
## Retreive LOTA source data from http://Processtrends.com
 series_link  <- "http://processtrends.com/files/RClimate_consol_temp_anom_latest.csv"
 series <- read.table(series_link,header=T, as.is=T,sep=",")
# n <- nrow(series)
 series_names<- c("GISS", "HadleyCRU", "NOAA", "RSS", "UAH")
 series_col <- c("red", "blue", "green", "yellow", "grey")

## Set data.frame to just 5 LOTA series 
 lota <- series[,c(1, 3:7)]

## Calculate series offsets 
  lota_s <- subset(lota, lota$yr_frac > baseline[1] & lota$yr_frac <baseline[2]+1)
  s_m <- numeric(5)
 ## Calculate series mean for baseline period
 for(i in 1:5) {
        s_m[i] <- mean(lota_s[,i+1],na.rm=T)
     }
 s_m
###############################################################################################
## Develop 5 LOTA Series with Common Baseline
   offset_mat <- as.matrix(lota)
 ## Build Matrix Columns
     for (j in 2:6) {
       offset_mat[,j] <- lota[,j] - s_m[j-1]
       }
   ## set column names
   m_names <- c("yr_frac", series_names)
   colnames(offset_mat) <- m_names
############################################################################################### 
## Generate 5 LOTA series trend chart - use x1 & x2 to specify plot xlim
  x1 <- 1979
  x2 <- 2011 
  y_lab <- paste("Anomaly ", "\u00B0", "C - Baseline (",baseline[1], " to ", baseline[2],")", sep="")
  plot_5_lota <- paste("Global LOTA Series Trends: ", x1, " to ", x2,"\nCommon Baseline (",baseline[1], " to ", baseline[2],")", sep="")
  par(oma=c(2,2,1,1)); par(mar=c(3,4,3,1)); par(ps=11)
  f_cb_1 <- paste(save_dir, "cb_1.png", sep="")
#  png(file=f_cb_1,width=1200,height=1000, res=200) 
  plot(offset_mat[,1], offset_mat[,2], col="red", type="n", xlab="", ylab= y_lab,
      xlim=c(x1,x2), ylim=c(-1, 1), main=plot_5_lota)
   abline(h=0, col="grey")
   legend("bottomright", series_names, col = series_col,
       text.col = "black", lty = c(1),pch=c(0),pt.cex=0,lwd=2,
       merge = F, bg = "white", bty="o", box.col = "black", cex = .7)
  points(offset_mat[,1], offset_mat[,2], type="l", col=series_col[1])
  points(offset_mat[,1], offset_mat[,3],type="l", col = series_col[2])
  points(offset_mat[,1], offset_mat[,4], type="l", col=series_col[3])
  points(offset_mat[,1], offset_mat[,5], type="l", col = series_col[4])
  points(offset_mat[,1], offset_mat[,6],type="l",  col=series_col[5])
  mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,0, adj = 0, cex = 0.8, outer=T)
  mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 0, adj = 1, cex = 0.8, outer=T)
 # dev.off()                                 # turn off PNG dev
  x11()
########################################################################################################
## Moving Avg Plot
 ma_per <- 133
 ma_yrs <- c(1880, 2011)
## Create time series of vector
  g_ts <- ts(offset_mat[,2], start = c(1880,1), freq=12)
  h_ts <- ts(offset_mat[,3], start = c(1880,1), freq=12)
  n_ts <- ts(offset_mat[,4], start = c(1880,1), freq=12)
  r_ts <- ts(offset_mat[,5], start = c(1880,1), freq=12)
  u_ts <- ts(offset_mat[,6], start = c(1880,1), freq=12)
## Create moving avg
  ma_g <- filter(g_ts, rep(1/ma_per,ma_per), sides=1)
  ma_h <- filter(h_ts, rep(1/ma_per,ma_per), sides=1)
  ma_n <- filter(n_ts, rep(1/ma_per,ma_per), sides=1)
  ma_r <- filter(r_ts, rep(1/ma_per,ma_per), sides=1)
  ma_u <- filter(u_ts, rep(1/ma_per,ma_per), sides=1)

## Plot moving avg series
 ma_min= min(c(ma_g,ma_h, ma_n, ma_r, ma_u),na.rm=T)*1.1
 ma_max <- max(c(ma_g,ma_h, ma_n, ma_r, ma_u),na.rm=T)*1.1
 ma_title <- paste("Global LOTA Series Trends\n",ma_per," Month Moving Avg,  Baseline period (",baseline[1], " to ", baseline[2],")",  sep="")
 par(lwd=2);par(oma=c(2,2,1,1)); par(mar=c(3,4,3,1)); par(ps=11)
 f_cb_4 <- paste(save_dir, "cb_4.png", sep="")
 png(file=f_cb_4,width=1200,height=1000, res=200) 

 plot(ma_g, col="red", type="l", xlab="", ylab = paste("Anomaly - ","\u00B0","C"),
   main = ma_title, ylim= c(ma_min, ma_max), xlim=ma_yrs)
 abline(h=0, col="grey",lwd=1)
 points(ma_g, col = series_col[1], type="l")
 points(ma_h, col = series_col[2], type="l")
 points(ma_n, col = series_col[3], type = "l")
 points(ma_r, col = series_col[4], type = "l")
 points(ma_u, col = series_col[5], type = "l")
 legend("bottomright", series_names, col = series_col,
       text.col = "black", lty = c(1),pch=c(0),pt.cex=0,lwd=2,
       merge = F, bg = "white", bty="o", box.col = "black", cex = .9)
 mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,0, adj = 0, cex = 0.8, outer=T)
 mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 0, adj = 1, cex = 0.8, outer=T)
 dev.off()    
 x11()
################################################################################################
## Generate trend plot of series ? and series ?
   p2_series <- c(3,4)
   p2 <- c(1998, 2011)                    # set plot yr range
   p2_mat <- subset(offset_mat, offset_mat[,1] >p2[1] & offset_mat[,1] <p2[2])
   p2_min= min(p2_mat[,p2_series[1]+1], p2_mat[,p2_series[2]+1],na.rm=T)*1.1
   p2_max <- max(p2_mat[,p2_series[1]+1], p2_mat[,p2_series[2]+1],na.rm=T)*1.1
   
   par(oma=c(2,2,1,1)); par(mar=c(3,4,3,1)); par(ps=11)
  f_cb_2 <- paste(save_dir, "cb_2.png", sep="")
  #png(file=f_cb_2,width=1200,height=1000, res=200) 

  plot_title <- paste("LOTA Series Comparison: ",series_names[p2_series[1]], " and ", series_names[p2_series[2]],
      "\nBaseline period (",baseline[1], " to ", baseline[2],")",  sep="")
  plot(p2_mat[,1], p2_mat[,p2_series[1]+1], type="n",
     col="grey", xlab="", ylab=paste("Anomaly - ","\u00B0","C"), 
      main=plot_title, xlim = p2, ylim <- c(p2_min, p2_max))
  points(p2_mat[,1], p2_mat[,p2_series[1]+1], type="l", col=series_col[p2_series[1]])
  points(p2_mat[,1], p2_mat[,p2_series[2]+1], type="l", col=series_col[p2_series[2]])
  legend("bottomright", c(series_names[p2_series[1]], series_names[p2_series[2]]), col = c(series_col[p2_series[1]], series_col[p2_series[2]]),
       text.col = "black", lty = c(1),pch=c(0),pt.cex=0,
       merge = F, bg = "white", bty="o", lwd=2,box.col = "black", cex = .8)
  baseline_note <- paste("Each series converted to common baseline:", baseline[1], " - ", baseline[2], sep="")
  text(p2[1], p2_max*0.95, baseline_note,adj=0, cex=0.7)
  mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,0, adj = 0, cex = 0.8, outer=T)
  mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 0, adj = 1, cex = 0.8, outer=T)
 # dev.off()    
x11()
########################################################################################################
## Generate scatter plot of series ? versus series ?
  xy_series <- c(3,4)
  par(oma=c(2,2,1,1)); par(mar=c(4,4,3,1)); par(ps=11)

  g_u_lm <- lm(offset_mat[,xy_series[1]+1]~ offset_mat[,xy_series[2]+1])
  reg_sum <- summary(g_u_lm)
  r2 <-  signif(reg_sum$adj.r.squared, digits=3)
  r2_note <- as.expression(substitute(r^2==r2, list(r2=r2)))
  xy_min= min(offset_mat[,xy_series[1]+1], offset_mat[,xy_series[2]+1],na.rm=T)
  xy_max <- max(offset_mat[,xy_series[1]+1], offset_mat[,xy_series[2]+1],na.rm=T)

  xy_title <- paste(series_names[xy_series[1]], " versus ", series_names[xy_series[2]],": ", p2[1], " - ", p2[2], 
      "\nBaseline period (",baseline[1], " to ", baseline[2],")",  sep="")
    f_cb_3 <- paste(save_dir, "cb_3.png", sep="")
#  png(file=f_cb_3,width=1200,height=1000, res=200) 
  plot(offset_mat[,xy_series[1]+1], offset_mat[,xy_series[2]+1], type="p",pch=16, cex=0.35, 
     col="grey", xlab=paste(series_names[xy_series[1]]," Anomaly \u00B0C"), 
     ylab=paste(series_names[xy_series[2]], " Anomaly \u00B0C"), 
      main=xy_title, ylim=c(xy_min,xy_max), xlim=c(xy_min,xy_max))
  abline(h=0, col="grey")
  abline(v=0 , col="grey")
  abline(g_u_lm, col="black")
  text(0.1,-0.4, r2_note, adj=0, cex=0.85)
  mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,0, adj = 0, cex = 0.8, outer=T)
  mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 0, adj = 1, cex = 0.8, outer=T)
 # dev.off()
##################################################################################################