############## RClimate Script: Compare JAXA_2007_2010.R ###################################
##                                                                                        ## 
## Download and process    JAXA daily SIE                                                 ##
## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ##
##                   http:chartsgraphs.wordpress.com    6/26/10                           ##
############################################################################################
  library(plotrix)
  options(digits=4)
  par(mfrow=c(1,1))
  par(ps=10); par(oma=c(0,0,0,0));par(mar=c(2,4,2,1))

  layout(matrix(c(1,2)), heights=c(2,1))
  #layout.show(2)

 link <- "http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv"
  j_data <- read.csv(link, header = F,
             colClasses = rep("numeric",3),
             comment.char = "#", na.strings = c(-9999),
             col.names = c("mo", "day", "yr", "JASIE") )

## Basic Date calcs
 j_date <- paste(j_data$mo,j_data$day, j_data$yr, sep="/")
 ## day of year
   day_val <- strptime(j_date, "%m/%d/%Y")  
   doy <- day_val$yday  +1 
 j_date <- as.Date(j_date, format="%m/%d/%Y")
 j_mo <- j_data$mo
 j_yr <- j_data$yr
 yr_frac <- as.numeric(j_data$yr + (j_mo-1)/12 + (j_data$day/30)/12)
 j_dy <- as.numeric(format(j_date,"%d"))

## data.frame of orig data & dates
  jd <- data.frame(j_date, yr_frac, j_yr,j_mo,j_dy,doy, j_data$JASIE/10^6)
  names(jd) <- c("dt", "yr_frac", "yr", "mo", "day","DOY", "sie")
 
## JAXA's source file inludes all days of year with NA's for upcomming days. Remove all NAs
   jds <- subset(jd, jd$sie != "NA")
   
## Get last day of jds data.frame
  no_row <- nrow(jds)
  cur_dt <- jds[no_row,1] 
  cur_day <- as.numeric(format(cur_dt, "%d"))
  cur_mo <-  as.numeric(format(cur_dt, "%m"))
  cur_yr <- as.numeric(format(cur_dt, "%y"))

## subset jds to make data.frames for all data for 2007 & 2010
 jd_2010 <- subset(jds, jds$yr == 2010)
  n_10 <-  nrow(jd_2010)
  peak_2010 <- max(jd_2010$sie, na.rm=T)

  peak_dt_2010 <-   format(jd_2010$dt[which(jd_2010$sie== peak_2010)], "%m/%d")
  peak_doy_2010 <-  jd_2010$DOY[which(jd_2010$sie== peak_2010)]

  last_dt <- format(jd_2010$dt[n_10], "%m/%d")
  last_doy <- jd_2010[n_10,6]
  last_2010 <- jd_2010[n_10,7]
  decr_2010 <- peak_2010 - last_2010

  jd_2007 <- subset(jds, jds$yr== 2007 & jds$DOY <= last_doy)
  n_07 <- nrow(jd_2007)
  peak_2007 <- max(jd_2007$sie, na.rm=T)
  peak_dt_2007 <-   format(jd_2007$dt[which(jd_2007$sie== peak_2007)], "%m/%d")
  peak_doy_2007 <-  jd_2007$DOY[which(jd_2007$sie== peak_2007)]

  last_2007 <- jd_2007[n_07,7]
  decr_2007 <- peak_2007 - last_2007

###  Day of year axis setup
## Set up basic day of year vectors (mon_names, 1st day of mon)
 mon_names <- c("Jan", "Feb", "Mar", "April", "May", "June", "July", "Aug", "Sept", "Oct","Nov","Dec")
 mon_doy <- c(1,32,60,91,121,151,182,213,244,274,305,335)
 mon_pos <- c(16, 46, 75, 106,135, 165, 200, 228, 255, 289, 320, 355)
 y_min <- (min(jd_2010$sie) )*0.9
 y_max <- (max(jd_2010$sie))* 1.1

# plot titles 
  main_title <- paste("Will JAXA 2010 Minimum Arctic Sea Ice Extent \nDrop Below Record 2007 Level?\n ",last_dt," Status", sep="")
  y_title <- "JAXA Arctic Sea Ice Extent - millions km^2"
# Plot  
 plot(jd_2010$DOY, jd_2010$sie, type="l", col = "red",
   main=main_title, cex.main=0.85,ylab = y_title, las=1, xlim = c(1, last_doy+2), 
   ylim = c(y_min, y_max  ), axes=F, xaxs="i", yaxs="i", las=1, cex.main = 0.8)
# custom x & y axes
  axis(side = 1, at=mon_doy, labels=F, xaxs="i")
  axis(side=1, at= mon_pos, labels=mon_names, tick=F, line=F, xaxs="i")
  axis(side=2, at=NULL, labels=T, yaxs="i", las=1)
  abline(v=1)    # need full length y axis line
# Add 2007 trend & 2007-2010 peak sie
  points(jd_2007$DOY, jd_2007$sie, type="l", col = "blue")
  points(peak_doy_2010, peak_2010, type = "p", col = "red", pch=16, cex=1.25)
  points(peak_doy_2007, peak_2007, type = "p", col = "blue", pch=16, cex=1.25)
  points(last_doy, last_2007, col = "blue", type="p", pch=18,cex=1.25)
  points(last_doy, last_2010, col = "red", type="p", pch=18, cex=1.25)

# create peak notes for legend  
  peak_note_2007 <- paste("2007 Peak @ ", signif(peak_2007, 4)," on ", peak_dt_2007, sep="")
  peak_note_2010 <- paste("2010 Peak @ ", signif(peak_2010, 4)," on ", peak_dt_2010,  sep="")
# create last date sie notes for legend
  last_note_2007 <- paste(last_dt, "/2007 @ ", signif(last_2007,3), sep="")
  last_note_2010 <- paste(last_dt, "/2010 @ ", signif(last_2010,3), sep="")

# legend
legend(5, 12, c("2007 Trend", "2010 Trend",peak_note_2007, peak_note_2010, last_note_2007, last_note_2010), col = c("blue", "red","blue", "red", "blue", "red"),
       text.col = "black", lty = c(1,1,0,0,0,0),pch=c(0,0,16,16,18,18),pt.cex=c(0,0,1,1,1,1),
       merge = F, bg = "white", bty="o", box.col = "white", cex = .8)


delta <- last_2010-last_2007
delta_note <- paste(last_dt, " Delta: \n2010 - 2007 =\n ", signif(delta,3), " million km^2", sep="")
text(25, 14.7, delta_note, cex=0.75, adj=0.5)


## function to calc slope last x days

  tr_func <- function(n) {
   last_n_2007 <- subset(jd_2007, jd_2007$DOY > last_doy-n)
   last_n_2010 <- subset(jd_2010, jd_2010$DOY > last_doy-n)

   lm_2007 <-  lm( sie ~ DOY, data = last_n_2007)
   tr_2007 <- coef(lm_2007)[2]

   lm_2010 <-  lm( sie ~ DOY, data = last_n_2010)
   tr_2010 <- coef(lm_2010)[2]
   my_df <- data.frame(as.integer(n), signif(tr_2007, 3), signif(tr_2010,3), row.names=NULL)
   return(my_df)
   }

 x <- c(2,3,5,7,10,20,30,60,90)
tr_df <- data.frame()
 how_many <- length(x)
 for (i in 1:how_many) {
  y <- x[i]
  results <- tr_func(y)
  tr_df[i,1] <-  y                #results[1]
  tr_df[i,2] <-  signif(results[2], 3)
  tr_df[i,3] <-  signif(results[3] ,3)
  tr_df[i,4] <-  signif(results[3]/results[2],3)
 }

names(tr_df) <- c("Prev days", "2007", "2010", "2010/2007")


par(mar=c(0,0,0.5,0))
plot(x=c(1,2,3), y= c(1,2,3), type="n", axes=F,ann=F, xlim=c(0,3), ylim=c(0,3))


table_title <- paste(last_dt, ": Arctic SIE Decline Rates for Number of Previous Days", sep="")
addtable2plot(0.825,2.95,tr_df,bty="o",display.rownames=F,hlines=T,cex=0.95,
  title=table_title, xjust=0, yjust=0)