############## RClimate Script: JAXA_faily.R             ###################################
## Script stored on http://chartsgraphs.wordpress.com account for Users to source()       ##                                                                                       ## 
## 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    ?/??/??                           ##
############################################################################################
  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="/")
 j_date <- as.Date(j_date, "%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,j_data$JASIE)
  names(jd) <- c("dt", "yr_frac", "yr", "mo", "day","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 get all data for current month
  mo_series <- subset(jds, jds$mo==cur_mo)
  mo_name <- format(as.Date(as.character(mo_series$dt[1]), format = "%Y-%m-%d"), format = "%B")

## Create vectors for monthly statistics
  mo_cur_day <- as.numeric(8)
  cd_yr_frac <- as.numeric(8)
 
## Calc monthly changes by year
  for (y in 2003:2010) {
   mo_data <- subset(mo_series, mo_series$yr == y)   
   mo_cur_day [y-2002] <- mo_data[cur_day,6] 
   cd_yr_frac[y-2002] <- mo_data[cur_day,2]
      }
mo_yr <- seq(2003,2010)
mo_sum <- data.frame(cd_yr_frac,mo_cur_day) 

## Plot titles
 cur_day_title <- paste(cur_mo, cur_day,cur_yr,sep="/") 
 last_val <- signif(mo_data[cur_day,6]/10^6,5)
 title <- paste("Daily JAXA Arctic Sea Ice Extent\n Month of ",mo_name,
     "\n", cur_day_title, " @ ", last_val, " million km^2",sep="")
 y_lab <- expression(paste("Arctic SIE - Million K",m^2, sep=""))

## Basic Plot
  par(mfrow=c(1,1)); par(ask=F); par(oma=c(3,2,1,1)) ; par(mar=c(2,4,3,1)); par(ps=10)
  plot(mo_series$yr_frac, mo_series$sie/10^6, type="p", pch=16, col="blue", cex=0.7,
    main=title, xlab="", ylab=y_lab, las=1, cex.main=0.85)
  points(mo_sum[,1], mo_sum[,2]/10^6, col="red", pch=16, cex=0.7)
  abline(h=last_val, col = "green")

## Plot Legend
  cdv <- paste("Day ", cur_day, " of ", mo_name, " for each year", sep="")
  cdline <- paste(cur_day_title, " @ ", last_val, " m km^2", sep="")
  legend(2002.25, 9.7, c("Daily SIE", cdv,cdline), col = c("blue", "red", "green"),
       text.col = "black", lty = c(0,0,1), pch=c(16,16,NA), 
       merge = T, bg = "white", bty="n", cex=0.8)

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