############## RClimate Script: JAXA_Daily_Sea_Ice_Extent.R             ##################
##                                                                                      ## 
## Download and process 2002 to 2012 daily   JAXA daily SIE                             ##
##  August 28, 2012                     by        Kelly O'Day                           ##
## http://chartgraphs.wordpress.com                                                     ##
##########################################################################################
  link <- "http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv"
  setwd("C:\\R_Home\\Charts & Graphs Blog\\RClimateTools\\Arctic_sea-ice_extent")
  what_date <- format(Sys.Date(), "%b %d, %Y")  # with month as a word

  title <- paste("IARC-JAXA Daily Arctic Sea Ice Extent\n",what_date) 
  par(oma=c(2,1,1,1)); par(mar=c(2,4,2,1))
  par(mgp=c(2.5,1,0))  # Number of margin lines for title, axis labels and axis line

# Data File: Month,Day,1980's Avg,1990's Avg,2000's Average,2002:2012
  j_full_data <- read.csv(link, header = F, skip=1, na.strings = c(-9999))
  j_data <- j_full_data[,c(1,2,6:16)] 
  colnames(j_data) <- c("mo", "day", "02", "03", "04", "05", "06", "07", "08", "09", "10", "11", "12")
   head(j_data)
  
# File has data for each day in 366 day year
# Establish Day of year
  for (i in 1:366)   j_data$yr_frac[i] <- i
# Convert ASIE to millions Km^2
  j_data[,c(3:13)] <- j_data[,c(3:13)]/1000000

###  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,366)
 mon_pos <- c(16, 46, 75, 106,135, 165, 200, 228, 255, 289, 320, 355)

# Loop through years 

for (j in 2002:2012)
    {
     png_name <- paste("asie",j,".png",sep="")
#  png(filename=png_name)
   which_yr <- j
   no_yrs <- j-2001
   col_off <- which_yr-1999
     
# Calc min asie for year
     min_asie <- min(j_data[,col_off], na.rm = T)  # must remove na's to get valid answer
     lab_asie <- round(min_asie,2)
     min_r <- which(j_data[,col_off] == min_asie)
     min_d <- j_data[min_r,2]
     min_m <- j_data[min_r,1]
     min_date <- paste(min_m,"/",min_d,"/",j, sep="")
  
  plot(j_data[,14],  type="n", col = "grey",axes=F, xlab="", 
     ylab="Arctic Sea Ice Extent - Millions Sq KM",
     ylim=c(0,15),xaxs="i", yaxs = "i",
     main=title)     
     text(20,1,"Data Source: http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv", cex = 0.8, adj=0,col = "black")
     mtext("D Kelly O'Day - http://chartsgraphs.wordpress.com", 1,0.5, adj = 0, cex = 0.8, outer=T)
  # 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)
     points(j_data[min_r,14], min_asie, col = "red",pch=19, cex = 2)
  # Add each preivous yr data series   
 for (n in 1:no_yrs)
 {
  points(j_data[,14], j_data[,2+n], col="grey", type="l",lwd=1)
  text(182,14,which_yr, col = "red", cex = 1.1)
    }
  points(j_data[,14], j_data[,no_yrs+2], col="red", type="l",lwd=2.5)
  text(182,14,which_yr, col = "red", cex = 1.1) 
  text(120,min_asie+0.35, min_date, col="red", cex=0.9)
  text(120,min_asie, lab_asie, col="red", cex=0.9)  
  
 # dev.off()
  }