# 3 Panel plot of RSS Temp Anomaly Trend, SATO Index, NINO3.4 Index
## Multiple trend chart version - Monthly anomaly data
## D Kelly O'Day; http://processtrends.com; http://chartsgraphs.wordpress.com
## Feb 6, 2009
#########################################################################
# Setup
#########################################################################
options(digits=6)
par(ps=10); par(cex.axis=c(1.1)); par(cex.lab=c(1));par(las="1")
 #Set outer and figure margins to reduce gap between plots
     par(mfrow=c(3,1)); par(oma=c(2.5,5,5,5));  par(mar=c(0,5,0,0))
link<- "C:\\R_Home\\Charts & Graphs Blog\\GISS_Temp_Anom\\Consol_GISS_SATO_Nino34A_RSS.csv"
## Start & end years for plot
    yr_st <- 1980
    yr_end<- 2010
    x_grid_gap <-5
##################################################################
# Read Data
##################################################################
ts_data <- read.csv(link,
             skip = 1, sep = ",", dec=".",
             row.names = NULL, header = FALSE,
             colClasses = rep("numeric",7),
             comment.char = "#", na.strings = c("*", "-",-99.9, -999.9),
             col.names = c("Yr", "Mo","Yr_frac", "GISS", "SATO", "NINO34", "RSS") )
attach(ts_data)
# Find last mo & year for title
  c <- nrow(ts_data)
 lst_mo <- ts_data$Mo[c]
 lst_yr <- ts_data$Yr[c]
 mo_names <- c("Jan", "Feb", "Mar", "Apr", "May", "June", "July", "Aug", "Sept",  "Oct", "Nov", "Dec")
 lst_RSS <- RSS[c]
 lst_NINO34 <- NINO34[c]
#################################################################
# Create Plot 1: RSS Temperature Anomaly
#################################################################
  plot(Yr_frac,RSS,  type = "S", col = "darkgrey",
     xlim = c(yr_st, yr_end),  ylim=c(-0.4, 1),
     xaxs="i", yaxs="i", axes= FALSE,
     ylab = expression(paste("RSS Anomaly ",degree, "C")))
     axis(1,  col = "lightgrey", at = NULL, labels = FALSE)
     axis(2, col = "lightgrey", labels = TRUE )
     abline(v=seq(yr_st,yr_end, x_grid_gap),col = "lightgrey")
     axis(3, col = "lightgrey", at = NULL, labels = TRUE)
     rect(1992,0.8, 2000, .95,col = "white", border = NA)
     text(1990, 0.87, "RSS Land & Ocean Temp Anomaly V3.2",pos = 4)
     rect(1982,0.55, 1986, .65,col = "white", border = NA)
     text(1983, 0.6, "Area: 70S to 82.5N", col = "grey")
     points( Yr_frac[c], lst_RSS, pch=19, col = "darkgrey")
     points(2001.7, -0.2, pch=19, col = "darkgrey")
     note <- paste(mo_names[lst_mo], ",", lst_yr, " @ ", lst_RSS,"C")
     rect(2001.9, -0.3, 2007, -0.1, col = "white",border = NA)
     text(2002, -0.2, note, pos = 4, col = "black", cex = 1)
########################################################################
# Create plot 2: Nino34 index
########################################################################
  n_red <- subset(ts_data, NINO34 >=0, select=c(Yr_frac, NINO34))
  n_blue <- subset(ts_data, NINO34<0, select=c(Yr_frac, NINO34))
  plot(Yr_frac, NINO34, type = "n", col = "grey",
     xlim = c(yr_st, yr_end), ylim=c( -2.5, 2.9),
     xaxs="i", yaxs = "i", axes = FALSE,  xlab ="",
     ylab=expression(paste("NINO3.4 Index - SST ", degree,"C")))
  points(n_red$Yr_frac, n_red$NINO34, type = "h", col = "red")
  points(n_blue$Yr_frac, n_blue$NINO34, type = "h", col = "blue")
  axis(2, col = "lightgrey", at = NULL, labels = TRUE)
  abline(h=0, col = "lightgrey",lwd = 1)
  abline(v=seq(yr_st, yr_end, x_grid_gap),col = "lightgrey")
  rect(1995,1.8, 1997,2.5,col = "white", border = NA)
  text(1995, 1.9, "NINO3.4 Index", pos = 3)
  clr="yellow"
  if(lst_NINO34 >0) clr="red" else clr="blue"
  points( Yr_frac[c], lst_NINO34, pch=19, col = clr)
     points(2001.7, -2, pch=19,col=clr)
     note_1 <- paste(mo_names[lst_mo], ",", lst_yr, " @ ", lst_NINO34,"C")
     rect(2001.9, -1.8, 2007, -2.3, col = "white",border = NA)
     text(2002, -2, note_1, pos = 4, col = "black", cex = 1)
#####################################################################
# Create plot 3: SATO index
#####################################################################
  par(mar=c(2,5,0,0))
  plot(Yr_frac, SATO, type = "h", col = "grey",
     xlim = c(yr_st, yr_end), ylim=c(0, 0.15),
     xaxs="i", yaxs = "i", axes = FALSE,  xlab ="",
     ylab=expression(paste("SATO Index ")))
  axis(1, col = "lightgrey", at = NULL, labels = TRUE )
  axis(2, col = "lightgrey", at = NULL, labels = TRUE)
  abline(v=seq(yr_st, yr_end, x_grid_gap),col = "lightgrey")
  rect(1996,0.04, 2009, 0.08,col = "white", border = NA)
  text(2002, 0.06,"Volcanoes increase SATO Index; \n Tend to decrease temperatures.\nNote Pinatobo effect in 1992")
  axis(3, col = "lightgrey", at = NULL, labels = FALSE)
  text(1988, 0.09,  "Pinatubo & \n Hudson \n Eruptions")
  arrows(1988.3,0.073, 1992, 0.05,  col = "black", length = 0.0)
  rect(1994.5,0.12, 2000, 0.145,col = "white", border = NA)
  text(1995.5, 0.125,  "SATO Index", pos = 3)
  text(1982, 0.12,  "El Chichon\n Eruption")
  arrows(1982,0.1, 1982.5, 0.08,  col = "black", length = 0.0)
############################################################################
# Overall Chart Title   & Footer Notes
############################################################################
  mtext("RSS and NINO3.4 Temperature Anomalies and SATO Index Trends Since 1980", 3,3, outer = TRUE)
  mtext(paste("Data Updated Through ", mo_names[lst_mo], ",", lst_yr), 3,2, outer = TRUE, cex = 0.7)
# Add KOD & print date to chart
   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)
 detach(ts_data)