##  RClimate_YTD_G_5_LOTA_for_latest_month.R             ###########################################################
##  Calc YTD temp anomalies for LOTA 5 global series: based on R script from N Urban                               #                                                                                              #                  
##  ClimateSight post comment by Nathan Urban http://climatesight.org/2010/10/18/odds-and-ends/#comments           # 
##  Nathan Urban    http://www.princeton.edu/~nurban/                                                              #  
##  Nathan Used my consolodated LOTA data to calc ytd trends by series                                             #
## Reads most recent LOTA file, calcs YTD for latest mo in file                                                    #
##  Orig 11/5/10                                                                                                   #
####################################################################################################################
  source("http://processtrends.com/Files/RClimate.txt")
  library(plotrix)                        # needed for addtable2plot() function
  series =  c("GISS", "HAD", "NOAA", "RSS", "UAH")                  
# read temperature anomaly data
  temps = read.csv("http://processtrends.com/Files/RClimate_consol_temp_anom_latest.csv")
# calculate years and months
  temps$yr = as.integer(temps$yr_frac) # round down fractional year
  temps$mo = temps$yr_mon - temps$yr*100 # extract month
  years = sort(unique(temps$yr)) # years of data
  last_year = years[length(years)] # last year in data
  which_mo <- temps$mo[length(temps$mo)]
  months = 1:which_mo 
  min_a <- min(temps[,seq(3,7)], na.rm=T)
  max_a <- max(temps[,seq(3,7)], na.rm=T)

## create matrix for ytd temps
  ytd_mat<-matrix(data = NA, nrow = length(years), ncol = 1, byrow = FALSE, dimnames = NULL)
  ytd_mat[,1] <- years

 par(mfcol=c(3,2)); par(mar=c(0.5,1,0,0)); par(oma=c(4.5,3.5,4,1))
 for (s in 1:5) {
  temp = temps[,s+2] # select temperature record to analyze
 
# calculate ytd means for selected months
  temp.ytd = sapply(years, function(y) mean(temp[temps$yr==y & temps$mo %in% months]))
  ytd_mat <- cbind(ytd_mat, temp.ytd)
  
################################################################################################################
  plot(years, temp.ytd, type="o", pch=16, col="grey", cex=0.3,las=1,
   xlab="", ylab="YTD Temperature Anomaly - C",  cex.main=0.85,axes=F,
   ylim=c(min_a, max_a))
  abline(h=0, col="grey")

# loess/30 years,  polynomial  degree 1
  regloess=loess(temp.ytd ~ years, span = 30/((max(years)-min(years))+1), degree = 1)
  v=predict(regloess, years)
  points(years, v, type = "l", lwd = 1.5, col= "blue")

## To highlight last point
 points(years[length(years)], temp.ytd[length(temp.ytd)], type = "p",pch=16, col = "red")
 last_temp.ytd <- signif(temp.ytd[length(temp.ytd)],2)
 last_val_note <- paste(last_year, " @ ", last_temp.ytd, sep="")
 ytd_note <- paste("YTD Average (Months 1 to ", which_mo,")", sep="")
 ## 3x2 panel axes: isees 3 & 4
  axis(3, at=NULL, labels=F, tcl=0.5)
  axis(4, at=NULL, labels=F, tcl=0.5)
 ## side 1 labels: yes if bottom panel, otherwidw no
  ifelse (s != 3 &  s!= 5,
    axis(1, at=NULL, labels=F, tcl=0.5),
    axis(1, at=NULL, labels=T))
 ## side 2 labels: yes if panel 1, 2 or 3, no if panel 4,5  
  ifelse (s <= 3, 
    axis(2, at=NULL, labels=T, las=1, cex=0.75),
    axis(2, at=NULL, tcl=0.5,labels=F))

  text(1885, 0.5, paste(s, " - ", series[s], sep=""),adj=0, font=2)
  text(1888, 0.325,  last_val_note, adj=0)
  points(1885, 0.325, pch=16, col="red")
  box()
 }
 m_title <- paste("Year To Date Land Ocean  Temperature Anomaly - C\nJanuary - ",mon_name[which_mo]," & Loess Regression (30 yr span)", sep="") 
 mtext(m_title, 3,1,outer=T, font=2, adj=0.5)
 mtext("YTD Temperature Anomaly - C", 2,2,adj=0.5, outer=T,cex=0.7,, las=0)
 mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com",1 ,2, adj = 0, cex = 0.7, outer=T)
 mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 2, adj = 1, cex = 0.7, outer=T)

# loop to calc ytd ranks by series
  r_df <- data.frame(years)
  colnames(ytd_mat) <- c("year", series)
  for (s in 1:5)  {
    s_rank <-rank(ytd_mat[,s+1], na.last="keep", ties.method="max")
    s_rank_max <- max(s_rank, na.rm=T)
    s_pos <- (s_rank_max +1)-s_rank
   r_df <- cbind(r_df, s_pos)
 }
  ytd_rank <- numeric(5)
  ytd_anom <- numeric(5)
  names(r_df) <- c("Year", series)
  num <- nrow(ytd_mat)
  ## loop to develop table of  latest ytd anomaly & ranking for 
   for(s in 1:5) {
    ytd_anom[s] <-   signif(ytd_mat[num,s+1],2)
    ytd_rank[s] <-  r_df[r_df$Year ==2010,s+1]
   }

#### Make table showing rankings 
  baseline <- c("1951-1980", "1961-1990", "1961-1990", "1979-1998", "1979-1998")
 st_yr <- c(rep(1880,3), rep(1979,2))
 ranking <- data.frame(series,  baseline, ytd_anom, ytd_rank)
 names(ranking) <- c("Series", "Baseline", "Anom-C","Ranking") 

 plot_tab <- function() {
   table_title <- paste("YTD though ", mon_name[which_mo], ", ",last_year, " Summary",sep="")
   plot(x=c(0,1,2,3), y= c(0,1,2,3), type="n", axes=F,ann=F, xlim=c(0,4), ylim=c(0,4))
   addtable2plot(2,1.5,ranking,bty="o",display.rownames=F,hlines=T,cex=1,
   title=table_title, xjust=0.5, yjust=0.5)
   }
  plot_tab()