####  Arctic_SIE_max_min_melt_by_yr.R                   #############################################
##  D Kelly O'Day  http://chartsgraphs.wordpress.com                                               ##
## RClimate script to calculate annual melt inArctic sea Ice Extent                                ##
## Use NSIDC annual data to calc max & min ASIE for each year                                      ##
## Original: 10/3/10                                                                               ##
#####################################################################################################
## Calc: melt(yr) = Max_SIE(yr) - Min_SIE(yr)
 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 = c("NA",-9999),
             col.names = c("yr_frac", "yr", "mo", "ext", "area"))
 sie <- sie[,c(1,2,3,4)]

## Calculate max, min & seasonal melts for eaqch year
 ann_max <- tapply(sie$ext, sie$yr, max, na.rm=T)
 ann_min <- tapply(sie$ext, sie$yr, min, na.rm=T)
 ann_melt = ann_max-ann_min
 ann_df <- data.frame(seq(1979,2010,1), ann_max, ann_min, ann_melt)
 names(ann_df)<- c("yr", "Max", "Min", "Melt")

## Create 3 panel chart showing max, min & seasonal melt
plot_func<- function() {
 par(mfrow=c(3,1)); par(las=1)
 par(oma=c(4,2,4,1)); par(mar=c(2,5,0,1))
 plot(ann_df$yr, ann_df$Max, type="o", pch=16, col="blue", xlab="", ylab="Annual Maximum SIE",axes=F)
    axis(1, at=NULL, labels=F)
    axis(2, at=NULL, labels =T)
    box(col="black")
    grid(nx=NULL, ny=NA, col = "grey", lty= "solid")
    rect(1979,14.6,1987 , 15, col = "white", border = "white")
    text(1979,14.8, "a) Maximum Monthly SIE", adj=0,font=2)
 plot(ann_df$yr, ann_df$Min, type="o", pch=16, col="green",xlab="", ylab="Annual Minimum SIE", axes=F)
    axis(1, at=NULL, labels=F)
    axis(2, at=NULL, labels =T)
    box(col="black")
    grid(nx=NULL, ny=NA, col = "grey", lty= "solid")
    rect(1979,4.5,1987, 4.9, col = "white", border = "white")
    text(1979,4.7, "b) Minimum Monthly SIE",adj=0,font=2)
 plot(ann_df$yr, ann_df$Melt, type="o", pch=16, col="red", xlab="", ylab="Annual Melt (Max - Min)")
   grid(nx=NULL, ny=NA, col = "grey", lty= "solid")
   rect(1979,7.5,1987 ,8.2, col = "white", border = "white")
   text(1979,7.8,"c) Seasonal Melt = Max - Min", adj=0, font=2)
# Generate and add bottom footer KOD & System data notes & overall chart title
       mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,1, adj = 0, cex = 0.7, outer=T)
       mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.7, outer=T)
       mtext("Data Source: NSIDC monthly data files: ftp://sidads.colorado.edu/DATASETS/NOAA/G02135/", 1,0, outer=T,adj=0.5, cex=0.7)
  m_title <- "Annual Arctic SIE Summary from NSIDC Monthly Data\n a) Maximum, b) Minimum, c) Melt (millions sq km)"
  mtext(m_title, 3,1, adj = 0.5, cex = 0.8, outer=T,font=2)
}

plot_func()
tail(ann_df,10)