##########################################################################################
#### R Script to consolidate PDO, NINO34, AMO  and GISS LOT anomaly data, generate a    ##
##              combined index and produce csv file                                     ##
##   KOD   Orig 4/21/10; Updated 5/21/10 to correct error in merge where 3 months were  ##
##   not being added to consolidated data.frame                                         ##
##########################################################################################
 library("plyr"); library("ggplot2")
###########PDO monthly data from U of Washington - JISAO #################################
  ## Download and process PDO Data File ###############
    p_url <- "http://jisao.washington.edu/pdo/PDO.latest"
    p_file <- c("pdo_latest.txt")
    download.file(p_url, p_file)
  ## Find number of rows in the file
    p_rows <- length(readLines(p_file)) -1
  ## Read file as  char vector, one line per row, Exclude first 32 rows 
    p_lines <- readLines(p_file, n=p_rows)[32:p_rows]
    p_lines2 <- p_lines[1:111]   # delete documentation lines @ end of file
  ##C onvert the character vector to a dataframe using fixed width format (fwf)
    p_df <- read.fwf(textConnection(p_lines2), widths = c(4,-3, rep(7,12)), 
       header=F, skip=0, colClasses = "numeric")
    closeAllConnections()
    names(p_df) <- c("Year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", 
           "Sep", "Oct", "Nov", "Dec")
  ## Convert from wide format to long format
    dfp <- melt(p_df, id.var="Year", variable_name="Month")
    names(dfp) <- c("yr", "mon", "pdo")
  ## dfp$Month is factor, Convert to month number & calc yr_frac
    pdo_mo_num <- as.character(unclass(dfp$mon)) 
    pdo_mo_frac <- (as.numeric(pdo_mo_num))/12  
    yr_frac <- dfp$yr + pdo_mo_frac
    yr_mo <- paste(dfp$yr , pdo_mo_num, sep="-")
    dfp <- data.frame(yr_frac,  yr_mo, dfp)
################ NINO34 Index #############################################################
    n_link <-   "http://climexp.knmi.nl/data/inino5.dat"
    nino_data <- read.csv(n_link,
             skip = 1, sep = "", dec=".", as.is = T,
             colClasses=rep("numeric",13),
             comment.char = "#", na.strings = c("*", "-","-99.9", "-999.90"),
             col.names = c("Yr","Jan", "Feb", "Mar", "Apr", "May", "June", "July", "Aug",
                         "Sept", "Oct", "Nov", "Dec"))
    dfn <- melt(nino_data, id.var=c("Yr"), measure.var=c(2:13))
    names(dfn) <- c("yr", "mon", "nino34")
    nino_mo_num <- as.character(unclass(dfn$mon))   
    yr_frac <- dfn$yr + (as.numeric(dfn$mon)-1)/12 + (15/30)/12    
    yr_mo <- paste(dfn$yr, nino_mo_num, sep="-")
    dfn <- data.frame(yr_frac, yr_mo,dfn)
    dfn <- dfn[order(dfn$yr_frac),]
    dfn <- subset(dfn, dfn$nino34 != "NA")
    dfn <- dfn[,c(1,2,5)]

################ AMO Index ###############################################################
    a_link <- "http://www.esrl.noaa.gov/psd/data/correlation/amon.us.long.data"
 ## Download and process AMO Data File ###############
    a_file <-    c("amo_latest.txt")
    download.file(a_link, a_file)
    a_rows <- length(readLines(a_file))
 ## Read file as  char vector, one line per row, Exclude first row 
    a_lines <- readLines(a_file, n=a_rows)
    num_a <- a_rows - 4
    a_lines_2 <- a_lines[2:num_a]
 ##Convert the character vector to a dataframe using fixed width format (fwf)
    a_df <- read.table(
     textConnection(a_lines_2), header=F, skip=0, colClasses = "numeric")
    closeAllConnections()
    names(a_df) <- c("Year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec")
  ## Convert from wide format to long format
    dfa <- melt(a_df, id.var="Year", variable_name="Month")
    names(dfa) <- c("yr", "mon", "amo")
  ## dfa$Month is factor, Convert to month number & calc yr_frac

   

    amo_mo_num <- as.character(unclass(dfa$mon))
    amo_mo_frac <- (as.numeric(amo_mo_num))/12  
    yr_mo <- paste(dfa$yr, amo_mo_num, sep="-")
    amo_mo_frac <- as.numeric( (amo_mo_num-0.5)/12   )
    yr_frac <- as.numeric(dfa$yr) + amo_mo_frac
    head(dfa)
   dfa <- data.frame(yr_frac, yr_mo, dfa[,3])
########### GISS Land & Ocean Anomaly Data #################################################
  ## GISS monthly data import script develodped by http://LearnR.wordpress.com   
    g_url <- c("http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt")
    g_file <- c("GLB.Ts+dSST.txt")
    download.file(g_url, g_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
     g_rows <- length(readLines(g_file)) - 12
  ## Read file as  char vector, one line per row, Exclude first 8 rows 
     g_lines <- readLines(g_file, n=g_rows)[9:g_rows]
  ## Data Manipulation, R vector 
  ## Use regexp to replace all the occurences of **** with NA
     g_lines2 <- gsub("\\*{3,5}", " NA", g_lines, perl=TRUE)
  ## Convert the character vector to a dataframe
     dfg <- read.table(
      textConnection(g_lines2), header=TRUE, colClasses = "character")
     closeAllConnections()
  ## Select monthly data in first 13 columns
     dfg <- dfg[,1:13]
  ## Convert all variables (columns) to numeric format
     dfg <- colwise(as.numeric) (dfg)
  ## Remove rows where Year=NA from the dataframe
     dfg <- dfg [!is.na(dfg$Year),] 
  ## Convert from wide format to long format
     dfmg <- melt(dfg, id.var="Year", variable_name="Month")
     GISS_mo_num <- unclass(dfmg$Month)
     GISS_mo_frac <- as.numeric((  unclass(dfmg$Month)-0.5)/12)
     GISS_yr_frac <- dfmg$Year + GISS_mo_frac
     GISS_anom <- dfmg$value/100
      yr_mo <- paste(dfmg$Year, GISS_mo_num, sep="-")


     dfmg <- data.frame(GISS_yr_frac, yr_mo,GISS_anom)
     dfmg <- dfmg[order(dfmg$GISS_yr_frac), ]
     dfmg <- dfmg[!is.na(dfmg$GISS_anom),]
     dfg <- data.frame(dfmg)
     names(dfg) <- c("yr_frac", "yr_mo","giss") 

############### Merge data frames  ###########################################################
  ## merge df's in sequence: dfp, dfn, dfa, fg
    pe_df <- merge(dfp, dfn,by="yr_mo")
    pea_df <- merge(pe_df, dfa, by="yr_mo")

    peag_df <- merge(pea_df, dfg, by="yr_mo")
    tail(peag_df,24)

    peag_df <- subset(peag_df, peag_df$yr_frac.x > 1900)
    peag_df <- subset(peag_df, peag_df$pdo != "NA")
    peag_df <- subset(peag_df, peag_df$nino34 != "NA")
    peag_df <- peag_df[order(peag_df$yr_frac), ]

peag_df <- peag_df[,c(1,2,5,7,9,11)]
names(peag_df) <- c("yr_mo", "yr_frac", "pdo", "nino34", "amo", "giss")
    tail(peag_df, 12)

############Calculate status indexes #########################################################
## Calculate pdo_sta, nino34_sta, amo_sta  indexes and combined pdo-nino34, pdo_amo indexes
   pdo_sta <-     ifelse(peag_df$pdo < 0,-1,1)
   nino34_sta <-  ifelse(peag_df$nino34 <0, -1, 1)
   amo_sta <-     ifelse(peag_df$amo <0, -1,1)
  ## pdo _score 
   pno_sc <- as.numeric()
   pn_en <-   ifelse(peag_df$pdo <= 0 & peag_df$nino34 <=0, 1,0)
   pn_ep <- ifelse(peag_df$pdo <=0 & peag_df$nino34 >0, 2,0)
   pp_en <- ifelse(peag_df$pdo >0 & peag_df$nino34 <=0,3,0)  
   pp_ep <- ifelse(peag_df$pdo>0 & peag_df$nino34 >0, 4,0)
   pn_score <- pn_en + pn_ep + pp_en+ pp_ep
   pn_an <-   ifelse(peag_df$pdo <= 0 & peag_df$amo <=0, 1,0)
   pn_ap <- ifelse(peag_df$pdo <=0 & peag_df$amo >0, 2,0)
   pp_an <- ifelse(peag_df$pdo >0 & peag_df$amo <=0,3,0)  
   pp_ap <- ifelse(peag_df$pdo>0 & peag_df$amo >0, 4,0)
   pa_score <- pn_an + pn_ap + pp_an+ pp_ap
   pn_an_en <- ifelse(peag_df$pdo <= 0 & peag_df$amo <=0 & peag_df$nino34 <=0, 1,0)
   pn_an_ep <- ifelse(peag_df$pdo <= 0 & peag_df$amo <=0 & peag_df$nino34 >0, 2,0)
   pn_ap_en <- ifelse(peag_df$pdo <=0  & peag_df$amo >0 & peag_df$nino34 <=0, 3,0)
   pn_ap_ep <- ifelse(peag_df$pdo <= 0 & peag_df$amo >0 & peag_df$nino34 >0, 4,0 )
   pp_an_en <- ifelse(peag_df$pdo > 0 & peag_df$amo <=0 & peag_df$nino34 <=0, 5,0)
   pp_an_ep <- ifelse(peag_df$pdo > 0 & peag_df$amo <=0 & peag_df$nino34 >0, 6,0)
   pp_ap_en <- ifelse(peag_df$pdo > 0 & peag_df$amo >0 & peag_df$nino34 <=0, 7,0)
   pp_ap_ep <- ifelse(peag_df$pdo > 0 & peag_df$amo >0 & peag_df$nino34 >0, 8,0 )

   pae_score <- pn_an_en + pn_an_ep + pn_ap_en + pn_ap_ep+pp_an_en + pp_an_ep + pp_ap_en + pp_ap_ep

### Consolidate data frame & write csv file #################################################
   peag_df <- data.frame(peag_df, pdo_sta, nino34_sta, amo_sta, pn_score, pa_score, pae_score)
   peag_df <- peag_df[order(peag_df$yr_frac), ]


 link_out <- "C:\\R_Home\\Charts & Graphs Blog\\RClimateTools\\PDO_ENSO_Combined\\pdo_enso_amo_comb_index_data.csv"
 write.csv(peag_df, file = link_out, row.names = F)
 library(RCurl)

#  Specify destination file path to send image file to Processtrends.com\images
 f_remote<- "http://processtrends.com/Files/pdo_enso_amo_comb_index_data.csv"
# http://processtrends.com/Files/pdo_enso_amo_comb_index_data.csv
# Upload to Processtrends.com with  call to  ftpUpload
     ftpUpload(link_out, f_remote, userpwd = "processt:5cathouse")