############## RClimate Script: PDO Trend                ###################################
## Script stored on http://chartsgraphs.wordpress.com account for Users to source()       ##                                                                                       ## 
## Download and process     ????????                                                      ##
## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ##
## Pacific Decadal Oscillation - PDO data from U of Wash, JISAO                           ## 
## Monthly data file w/ Documentation strings at top & bottom of file                     ##
##                   http:chartsgraphs.wordpress.com    4/18/2010                         ##
############################################################################################

## Monthly data in columns
  library("plyr"); library("ggplot2")
#url <- "ftp://ftp.ncdc.noaa.gov/pub/data/ersst-v2/pdo.1854.latest.ts"
  url <- "http://jisao.washington.edu/pdo/PDO.latest"
## Download and process PDO Data File ###############
  file <- c("pdo_latest.txt")
  download.file(url, file)
## Read & Cleanup File
  # Find number of rows in the file
    rows <- length(readLines(file)) -1
# Read file as  char vector, one line per row, Exclude first 32 rows 
  lines <- readLines(file, n=rows)[32:rows]
  lines2 <- lines[1:111]   # delete documentation lines end of file

#Convert the character vector to a dataframe using fixed width format (fwf)
  df <- read.fwf(
  textConnection(lines2), widths = c(4,-3, rep(7,12)), header=F, skip=0, colClasses = "numeric")
  closeAllConnections()
 names(df) <- c("Year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")

# Convert from wide format to long format
  dfm <- melt(df, id.var="Year", variable_name="Month")
  names(dfm) <- c("Year", "Month", "pdo")
# dfm$Mointh is factor, Convert to month number & calc yr_frac
    pdo_mo_num <- unclass(dfm$Month)
    pdo_mo_frac <- as.numeric( (pdo_mo_num-0.5)/12   )
    yr_frac <- as.numeric(dfm$Year) + pdo_mo_frac
# build consolidated data.frame 
  dfm <- data.frame(yr_frac,  dfm)
  dfm <- dfm[order(dfm$yr_frac), ]
  dfm <- dfm[!is.na(dfm$pdo),]
# develop positive and negative pdo data.frames
  pos <- subset(dfm, dfm$pdo>= 0)
  neg <- subset(dfm,dfm$pdo<0)

# calc pdo moving average
  map <- 13
  ma_text <- paste(map, " Month Moving Avg", sep="")
  ma <- filter(dfm$pdo, rep(1/map,map), sides=1)
  ma <- ts(ma, start=c(1900,1), freq=12)

## Determine Month, Year for last reading
    c <- nrow(dfm)                    # Find number of data rows
    lmo <- dfm$Month[c]                           # Find last month of data
    lyr <- dfm$Year[c]                           # Find last year od data
    lpdo <- dfm$pdo[c]                         # Find last pdo reading
    lyr_frac <- dfm$yr_frac[c]

####### Plot pdo series
 par(las=1); par(ps=10); par(oma=c(2.5,1,1,1)); par(mar=c(2,4,2.5,1))
 main_title <- paste("Pacific Decadal Oscilliation (PDO)\n Univ of Washington, JISAO: Jan., 1990 to ",lmo, ", ", lyr, sep="")
     
 plot(pos$yr_frac, pos$pdo, type="n", col = "grey", main=main_title,cex.main=0.9,
   xlim = c(1900,2011), ylim=c(-3.5, 5), xlab="", ylab= "PDO Index")
 ## add grey lines for pdo shifts in 1925, 1947, 1977 
   abline(v=1925, col = "grey")
   abline(v=1947, col = "grey")
   abline(v=1977, col = "grey")
 ## add pos & neg pdo series and moving avg
   points(neg$yr_frac, neg$pdo, type="h", col = "lightblue")
   points(pos$yr_frac, pos$pdo, type="h", col = "pink")
   points(ma, type="l", col = "black", lwd=1.5)

  ## vertical line across top & label pdo phases
   points(x=c(1900, 2011), y=c(4,4), type="l", col = "grey")
   text(1936,4.5,"Warm\nPhase", cex=0.8, col = "black")
   text(1962,4.5,"Cool\n Phase", cex=0.8, col = "black")
   text(1995,4.5,"Warm\n Phase", cex=0.8, col = "black")
   text(2010, 4.5, "?", adj=1, cex=1.5, col = "red")
  
  ## highlight last pdo reading
    points(lyr_frac, lpdo, col = "red", type = "p", pch=16) 
    l_val_txt <-   paste(lmo,", ",lyr, ": ", lpdo,sep="")   # Note on last data point
   
  legend(1898,-2.1, c("Positive Monthly PDO", "Negative Monthly PDO", ma_text, l_val_txt ), col = c("pink", "lightblue","black", "red"),
       text.col = "darkgrey", lty = c(1,1,1,0),lwd=c(2), pch=c(0,0,0,16),pt.cex=c(0,0,0,1),
       merge = F, bg = "white", bty="o", box.col="white",cex = .725)

#############################################################################################
# Generate and add bottom footer KOD , source, System date notes
       source_note <- paste("Data Source: ", url)
       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)
       mtext(source_note, 1,0,adj=0.5, cex=0.8, outer=T)
#############################################################################################

tail(dfm,12)