RClimate Script: Pacific Decadal Oscillation Trend

This RClimate Script lets users retrieve and plot the monthly and moving average  Pacific Decadal Oscillation (PDO) data  from the University of Washington’s JISAO website. The script retrieves the PDO data from January, 1900 until latest month available at time script is run.  The trend chart shows the JISAO PDO trend and user selected moving average period.

Brief Introduction to Pacific Decadal Oscillation (PDO)

I have read so many blog posts about the PDO’s potential role in global warming that I have decided to learn about it for myself.  While I plan to examine the PDO – ENSO – global warming relationships in future posts, readers can get a preview of my inquiry from the posts that have gotten me started.

  1. First, I would read Tamino’s post  at Open Mind at this link. This is probably the best place to start if you are new to PDO.
  2. Next, I would go to the data source, the JISAO site to get a complete grounding in PDO.
  3. Third, I’d go to Bob Tisdale’s Misunderstandings About the PDO – Revisited post to get a clearer picture of what the PDO is and what it is not.

PDO Trend Chart

Here’s my R script based JISAO’s PDO trend by month chart  since 1900.

This chart shows the January, 1900 to March 2010 PDO trend, the 13 monthly moving average PDO, with the latest monthly value highlighted in red.

Here is the Data Link.

Here is the R Script


## Monthly data in columns

library("plyr"); library("ggplot2"); library(reshape)
<pre>
#url <- "ftp://ftp.ncdc.noaa.gov/pub/data/ersst-v2/pdo.1854.latest.ts"
#  url <- "http://jisao.washington.edu/pdo/PDO.latest"
  url <- "http://www.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))
# Read file as  char vector, one line per row, Exclude first 32 rows
  lines <- readLines(file, n=rows)[32:rows]
 ## Need trial & rerro to get last row of data versus documentation text
  lines2 <- lines[1:114]   # delete documentation lines end of file
 # tail(lines2,10)

#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
p_func <- function() {
 par(las=1); par(ps=12); 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., 1900 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, c("Positive Monthly PDO", "Negative Monthly PDO", ma_text, l_val_txt ), col = c("pink", "lightblue","black", "red"),
       text.col = "black", lty = c(1,1,1,0),lwd=c(2), pch=c(0,0,0,16),pt.cex=c(0,0,0,1),
       merge = T, bg = "white", bty="o", box.col="white",cex = .775)

#############################################################################################
# Generate and add bottom footer KOD , source, System date notes
       source_note <- paste("Data Source: ", url)
       mtext("D Kelly O'Day - http://chartsgraphs.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)
#############################################################################################

}

This R script includes several tools that may be helpful for R beginners, including:

  • download.file() to retrieve on-line file
  • readLines() to  handle messy file with mix of documentation and data lines
  • read.fwf() to handle ** characters mixed with numeric  data
  • melt() to convert from wide format data set to long format
  • subset() to separate positive and negative PDO series
  • filter() to calculate moving average for user specified period
  • abline() to add grid lines at PDO regime shift years
  • legend() to provide detailed legend
  • designation of last month in chart title, legend
  • designation of last point in series
About these ads

One response to “RClimate Script: Pacific Decadal Oscillation Trend

  1. Pingback: Checking Do-It-Yourself Climate Science | Climate Charts & Graphs

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s