############## RClimate Script: GISS Monthly Temperature Anomaly ###########################
## Script stored on http://chasrtsgraphs.woedpress.com account for Users to source()      ##                                                                                       ## 
## Download and process GISS Monthly Anomaly Data File                                    ##
## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ##
##                   http:chartsgraphs.wordpress.com    2/5/10                            ##
############################################################################################
   library(plyr); library(reshape); par(las=1); par(ps=10)
  par(oma=c(2.5,1,1,1)); par(mar=c(2.5,4,2,1))
   
## File Download and File
## GISS monthly data import script develodped by http://LearnR.wordpress.com   

  url <- c("http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt")
  file <- c("GLB.Ts+dSST.txt")
  download.file(url, file)
  
## 1st 8 rows and the last 12 rows contain instructions 
## Find out the number of rows in a file, and exclude the last 12
    rows <- length(readLines(file)) - 12
## Read file as  char vector, one line per row, Exclude first 8 rows 
    lines <- readLines(file, n=rows)[9:rows]
## Data Manipulation, R vector 
## Use regexp to replace all the occurences of **** with NA
    lines2 <- gsub("\\*{3,5}", " NA", lines, perl=TRUE)
## Convert the character vector to a dataframe
    df <- read.table(
      textConnection(lines2), header=TRUE, colClasses = "character")
    closeAllConnections()
## Select monthly data in first 13 columns
    df <- df[,1:13]
## Convert all variables (columns) to numeric format
    df <- colwise(as.numeric) (df)
## Remove rows where Year=NA from the dataframe
    df <- df [!is.na(df$Year),] 
## Convert from wide format to long format
    dfm <- melt(df, id.var="Year", variable_name="Month")
    GISS_mo_num <- unclass(dfm$Month)
    GISS_mo_frac <- as.numeric((  unclass(dfm$Month)-0.5)/12)
    GISS_yr_frac <- dfm$Year + GISS_mo_frac
    GISS_anom <- dfm$value/100
    dfm <- data.frame(dfm, GISS_mo_num, GISS_yr_frac, GISS_anom)
    dfm <- dfm[order(dfm$GISS_yr_frac), ]
    dfm <- dfm[!is.na(dfm$GISS_anom),]
## Find last report month and last value
    GISS_last <- nrow(dfm)
    GISS_last_yr <- dfm$Year[GISS_last]
    GISS_last_mo <- dfm$Month[GISS_last]
    GISS_last_mo_num <- dfm$GISS_mo_num[GISS_last]
    GISS_last_yr_frac <- GISS_last_yr + (GISS_last_mo_num-0.5)/12
    GISS_last_temp <- dfm$GISS_anom[GISS_last]

## Calc decade averages
  dec_mean<- as.numeric(14)
  dec_st <- as.numeric(14)
  dec_end <- as.numeric(14)
  yr_n <- as.integer(dfm$GISS_yr_frac)
  base_yr <- 1870
  dec_n <-  (as.numeric((yr_n - base_yr) %/% 10) * 10) + base_yr
  dfm <- data.frame(dfm, dec_n)
  for (i in 1:13) {dec_st[i] = base_yr+ i*10
     dec_sub <- subset(dfm, dec_n == dec_st[i], na.rm=T)
     dec_mean[i] <- mean(dec_sub$GISS_anom)
     }
 dec_st[14] <- 2010              # Need to have for last step line across decade
 dec_mean[14] <- dec_mean[13]
 dec<- data.frame(dec_st, dec_mean)

# Trend chart with trend line 
# specify plot yr min & max
  p_xmin <- 1880;   p_xmax <- GISS_last_yr_frac
  title <- paste("GISS Land and Sea Temperature Anomaly Trend \n", p_xmin, " to ", GISS_last_mo, ", ", GISS_last_yr)
  plot(dfm$GISS_yr_frac, dfm$GISS_anom, type = "l", col = "grey",
     xlim = c(p_xmin, p_xmax), ylab = "Temperature Anomaly - C (1951-1980 Baseline)",
     xlab="", main = title,cex.main = 0.85)
  points(GISS_last_yr_frac, GISS_last_temp, col = "red", pch=19)
  last_pt <- paste( GISS_last_mo, ",", GISS_last_yr, "@ ", GISS_last_temp, "C")
  points(dec$dec_st, dec$dec_mean, type="s", col="blue")

## add legend
  legend(1880,0.9, c("Decade Mean Anomaly", "Monthly Anomaly" ,last_pt), col = c("blue", "grey", "red"),
       text.col = "black", lty = c(1,1,0),pch=c(0,0,16),pt.cex=c(0,0,1),
       merge = F, bg = "white", bty="o", cex = .75, box.col="white")

 out <- paste("GISS Monthly Temperature Anomaly \nData updated through:" , GISS_last_mo, ",", GISS_last_yr)
 cat(out, fill = 40)
 t_pos <- p_xmin + 0.5*(p_xmax-p_xmin)
 text(t_pos, -0.75, out, cex = 0.7, adj = 0)
 data_source <- "Data Source: http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt"

# Plot Annotation
 mtext(data_source,1,0, cex = 0.75, adj=0.5, outer=T)
 mtext("D Kelly O'Day - http://chartgraphs.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)