############## RClimate Script: Arctic Sea Ice Extent Trend by Month  ######################
## Script stored on http://Processtrends.com for Users to source()                        ##                                                                                       ## 
## Source data from NSIDC organized into 12 monthly files                                 ##
## I have built and posted a consoldidated monthly file  available at this link           ##
## Data link:   http://processtrends.com/files/RClimate_NSIDC_sea_ice_extent.csv          ##
## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ##
##                   http:chartsgraphs.wordpress.com    2/9/10                            ##
############################################################################################

  par(las=1); par(ps=10); par(oma=c(3.5,2,1,1)); par(mar=c(2,4,2,1))
  par(bg="white")
  library(zoo)

 link_n <- "http://processtrends.com/files/RClimate_NSIDC_sea_ice_extent.csv"
 sie <- read.table(link_n, skip = 1, sep = ",", header = F,
             colClasses = rep("numeric", 5),
             comment.char = "#", na.strings = "NA",
             col.names = c("yr_frac", "yr", "mo", "ext", "area")   )
  sie <- sie[,c(1,2,3,4)]
  sie <- subset(sie, yr_frac >= 1979) # get 1st full year of data

## Determine month - yr of last reading
    lr <- nrow(sie)                        # Find number of data rows
    l_mo <- sie$mo[lr]                     # Find last month of data
    l_yr <- sie$yr[lr]                     # Find last year od data
    l_yr_frac <- sie[lr,1]
    l_sie <- signif(sie$ext[lr],4)       # Find last sie reading

## Create plot
plot_func <- function() {
  par(las=1); par(ps=10); par(oma=c(3.5,2,1,1)); par(mar=c(2,4,2,1))
  par(bg="white")
  plot_col<- "darkgrey"
  title = paste("Arctic Sea Ice Extent Trend by Month\nBased on NSIDC Monthly Data: 1/1979 to ",
     l_mo, "/", l_yr, sep="")
  y_title <- expression(paste("Monthly Arctic Sea Ice Extent - Millions k",m^2))
 plot(sie$yr_frac, sie$ext, col = "grey", type = "n", cex.lab = 0.85,
   xlim = c(1979, 2013), xlab = "",  ylab = y_title, main = title, cex.main=0.9)
 # loop through months of year to show individual monthly trends 
  for (i in 1:12) {
    mon <- subset( sie, sie$mo == i)
    points(mon$yr_frac, mon$ext, type="l", col = plot_col)
    row_num <- nrow(mon)
    text(l_yr+ 1.5, mon[row_num,4],i,adj=1,col="blue", cex=0.65)
    } 
   points(l_yr_frac, l_sie, col="red", pch=16)  # show most recnet mon value in red
   last_mo <- subset(sie, sie$mo==l_mo)         # subset mon to show entire series in red
   points(last_mo$yr_frac, last_mo$ext,type="l", col="red")
  # Legend: create mo names vector to label most recent month
    mo_names <- c("Jan", "Feb", "Mar", "Apr", "May", "June", "July", "Aug", "Sept",  "Oct", "Nov", "Dec")
   legend(1980, 6.5, c("Monthly SIE", paste(l_mo, " - Month number", sep=""),
     paste("Trend for ",mo_names[l_mo], sep=""), paste("Most Recent ", mo_names[l_mo]," SIE", sep="" )), 
       text.col = c("darkgrey", "blue","red", "red"), lty = c(1,0,1,0),pch=c(0,0,0,19),
       col=c("grey","grey", "red", "red"),pt.cex=c(0,0,0,1),
       merge = F, bg = "white", bty="n", cex = .7)

## Plot Annotations
  source_1 <- "Original Monthly Files: NSIDC Data Source files: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/"
  source_2 <- "ProcessTrends.com consolidated time series file:"
  source_3 <- link_n
  mtext(source_1, 1,0,  outer=T, adj=0.5, cex=0.7)
  mtext(source_2,1, .75, outer=T, adj=0.5, cex = 0.7)
  mtext(source_3, 1,1.5, outer=T, adj=0.5,cex = 0.7)
  mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,2.25, adj = 0, cex = 0.8, outer=T)
  mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 2.25, adj = 1, cex = 0.8, outer=T)
  }

     x11()     # turns grahics console on

#  Plot on graphics console
     plot_func()