# RSS MSU Satellite Data - TLT - Temp Anomaly Trend - Lower Troposhpere - Land & Ocean C
## R script by D Kelly O'Day, 2/5/09
## Get Monthly anomaly data
#Set plot parameters
     par(mfrow=c(1,1))
     par(oma=c(3,0,4,1)); par(mar=c(2,5,0,0))
     par(las=1); par(ps = c(9)) ; par(bty="l")
link1<- "http://www.remss.com/data/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_TLT_Anomalies_Land_and_Ocean_v03_2.txt"
mo_RSS_in <- read.table(link1,
             skip = 3, sep = "", dec=".",
             row.names = NULL, header = FALSE,
             as.is = T,
             colClasses = c(rep("numeric",3), rep("NULL", 8)),
             comment.char = "#", na.strings = c("*", "-",-99.9, -999.9),
             col.names = c("yr", "mo", "RSS_lo", rep("",8)))
    yr_frac <- mo_RSS_in$yr + (mo_RSS_in$mo-1)/12     # yr_frac simplifies calcs
# Create new data frame   & attach
  mo_RSS_lo <- data.frame(mo_RSS_in, yr_frac)
  attach(mo_RSS_lo)
# Get last mo, yr and value
  c <- nrow(mo_RSS_lo)                    # Find number of data rows
    lmo <- mo_RSS_lo$mo[c]                           # Fimd last month of data
    lyr <- mo_RSS_lo$yr[c]                           # Find last year of data
    lRSS <- RSS_lo[c]                         # Find last reading
# Create RSS plot
   plot(yr_frac, RSS_lo,  type="l",  col ="grey60",  xlab = "",
     ylab = expression(paste("RSS Temperature Anomaly - ",degree,"C")),
     xlim=c(1979, 2010), ylim=c(-0.6, 0.80),
     xaxs="i", yaxs="i", axes = TRUE,
     cex.axis = 0.95, cex.lab = 0.95)
     abline(h=0, col = "grey")
#     axis(1,  col = "grey", at = NULL, labels = TRUE, cex=0.75)
#     axis(2, col = "grey", at = NULL, labels = TRUE)
# Calculate RSS Linear Model for rate of change
    lm_fit <- lm(RSS_lo ~ yr_frac)
        a <- coef(lm_fit)[1]
        b <- coef(lm_fit)[2]
   ## Determine min & max yr_frac
   yr1 <- min(yr_frac)
   yr2 <- max(yr_frac)
        y1 <- a + b* yr1
        y2 <- a + b *yr2
        x_val <- c(yr1, yr2)
        y_val <- c(y1, y2)
      lines(x_val,y_val, type = "l",  col = "red")
      b <- signif(b, 3)
# Plot annotation
  n_text <- 0.9                 # use to control note font size
  # Establish plot titles
  title_3 <- expression(paste("RSS Land & Ocean Temperature Anomalies - ",degree, "C"))
  title_2 <- paste("Period: 1/1979  to ",lmo,"/",lyr)
   mtext(title_3, side = 3, line = 3, adj = 0.5)
   mtext(title_2, side = 3, line = 2, adj = 0.5, cex = n_text)
   mtext( "RSS Version 3.2 - MSU/AMSU TLT", side = 3,line = 1, adj=.5, cex = n_text)
  # Establish plot notes
   text(1995,-0.4, "Overall slope = ", pos = 1, col = "red", cex = n_text)
   text(1998.7, -0.4, b, pos = 1, col ="red", cex = n_text)
   text(2002, -0.4, expression(paste(degree, "C per year")), pos = 1,cex = n_text, col = "red")
   text(1980, 0.7, " Area: 70S to 82.5N", col = "black", pos =4, cex = n_text)
  # Highlight last reading
     points( yr2, lRSS, pch=19, col = "black")
     points(1996, -0.5, pch=19)
     note <- paste(lmo, "/", lyr, " @ ", lRSS," C")
      text(1996, -0.5, note, pos = 4, col = "black", cex = n_text)
  # Add margin notes for name and plot date
      mtext("D Kelly O'Day, http://ChartsGraphs.Wordpress.Com", side = 1, line = 3, adj = 0, cex = 0.9)
      mtext(format(Sys.time(), "%m/%d/ %Y"), side = 1, line = 3, adj = 1, cex = 0.9)
detach(mo_RSS_lo)