Results from Rweb

You are using Rweb1.04 on the server at www.nku.edu

Rweb:>   
Rweb:> ############## RClimate Script: Compare JAXA_2007_2010.R ################################### 
Rweb:> ##                                                                                        ##  
Rweb:> ## Download and process    JAXA daily SIE                                                 ## 
Rweb:> ## Developed by D Kelly O'Day to demonstrate use of source() function for climate data    ## 
Rweb:> ##                   http:chartsgraphs.wordpress.com    6/26/10                           ## 
Rweb:> ############################################################################################ 
Rweb:> library(plotrix) 
Rweb:> # ael 
Rweb:> library(reshape) 
Rweb:> # lea 
Rweb:> options(digits=4) 
Rweb:> par(mfrow=c(1,1)) 
Rweb:> par(ps=10); par(oma=c(0,0,0,0));par(mar=c(2,4,2,1)) 
Rweb:>  
Rweb:> layout(matrix(c(1,2)), heights=c(2,1)) 
Rweb:>  
Rweb:> #ael: this is the file, which is kept current (I pulled the data from here) 
Rweb:> #  link <- "http://www.ijis.iarc.uaf.edu/seaice/extent/plot.csv" 
Rweb:> #ael: updated 
Rweb:> #ael: link <-"/var/www/html/data/Rweb/RClimate/files/plot.csv.new" 
Rweb:> link <-"/var/www/html/data/Rweb/RClimate/files/plot_extent_n_v2.csv" 
Rweb:>  
Rweb:> j_data <- read.csv(link, header = T, 
+              colClasses = rep("numeric",17), 
+              comment.char = "#", na.strings = c(-9999), 
+              col.names = c("Month","Day","1980's Average","1990's Average","2000's Average","2002","2003","2004","2005","2006","2007","2008","2009","2010","2011","2012","2013","2014","2015","2016") ) 
Rweb:>  
Rweb:> #ael -- this is my approach to getting in all the data, using the melt command 
Rweb:> #       which I just discovered today, Wed Dec 17 19:26:38 EST 2014:  
Rweb:> jdl <- melt(j_data,id.var=c("Month","Day"),variable_name="Year") 
Rweb:> names(jdl) <- c("mo", "day","yr","JASIE") 
Rweb:> jdl$yr <- gsub(".s.Average", "", jdl$yr) 
Rweb:> jdl$yr <- gsub("X", "", jdl$yr) 
Rweb:> j_data <- jdl; 
Rweb:> #lea 
Rweb:>  
Rweb:> ## Basic Date calcs 
Rweb:> j_date <- paste(j_data$mo,j_data$day, j_data$yr, sep="/") 
Rweb:> ## day of year 
Rweb:> day_val <- strptime(j_date, "%m/%d/%Y")   
Rweb:> doy <- day_val$yday  +1  
Rweb:> j_date <- as.Date(j_date, format="%m/%d/%Y") 
Rweb:>  
Rweb:> j_mo <- j_data$mo; 
Rweb:> j_yr <- as.numeric(j_data$yr); 
Rweb:> yr_frac <- as.numeric(j_yr + (j_mo-1)/12 + (j_data$day/30)/12) 
Rweb:> # ael: needed to turn strings into numbers: 
Rweb:> j_dy <- as.numeric(format(j_date,"%d")) 
Rweb:>  
Rweb:> ###################################################################### 
Rweb:> # ael: just writing out the current plot.csv file, for other files to use: 
Rweb:> # I remove the 1980, 1990, and 2000 average data. So we start with daily data from 2002. 
Rweb:> # jd <- data.frame(j_mo,j_dy,j_yr,j_data$JASIE) 
Rweb:> # write.csv(jd, file = link_out, row.names = F) 
Rweb:> # lea 
Rweb:> ###################################################################### 
Rweb:>  
Rweb:> ## data.frame of orig data & dates 
Rweb:> jd <- data.frame(j_date, yr_frac, j_yr,j_mo,j_dy,doy, j_data$JASIE/10^6) 
Rweb:> names(jd) <- c("dt", "yr_frac", "yr", "mo", "day","DOY", "sie") 
Rweb:>  
Rweb:>   
Rweb:> ## JAXA's source file inludes all days of year with NA's for upcomming days. Remove all NAs 
Rweb:> jds <- subset(jd, jd$sie != "NA") 
Rweb:>     
Rweb:> ## Get last day of jds data.frame 
Rweb:> no_row <- nrow(jds) 
Rweb:> cur_dt <- jds[no_row,1]  
Rweb:> cur_day <- as.numeric(format(cur_dt, "%d")) 
Rweb:> cur_mo <-  as.numeric(format(cur_dt, "%m")) 
Rweb:> cur_yr <- as.numeric(format(cur_dt, "%y")) 
Rweb:>  
Rweb:> ## subset jds to make data.frames for all data for 2007 & 2010 
Rweb:> jd_2010 <- subset(jds, jds$yr == 2010) 
Rweb:> n_10 <-  nrow(jd_2010) 
Rweb:> peak_2010 <- max(jd_2010$sie, na.rm=T) 
Rweb:> peak_dt_2010 <-   format(jd_2010$dt[which(jd_2010$sie== peak_2010)], "%m/%d") 
Rweb:> peak_doy_2010 <-  jd_2010$DOY[which(jd_2010$sie== peak_2010)] 
Rweb:> pit_2010 <- min(jd_2010$sie, na.rm=T) 
Rweb:> pit_dt_2010 <-   format(jd_2010$dt[which(jd_2010$sie== pit_2010)], "%m/%d") 
Rweb:> pit_doy_2010 <-  jd_2010$DOY[which(jd_2010$sie== pit_2010)] 
Rweb:> last_dt <- format(jd_2010$dt[n_10], "%m/%d") 
Rweb:> last_doy <- jd_2010[n_10,6] 
Rweb:> last_2010 <- jd_2010[n_10,7] 
Rweb:> decr_2010 <- peak_2010 - last_2010 
Rweb:>  
Rweb:> jd_2007 <- subset(jds, jds$yr== 2007 & jds$DOY <= last_doy) 
Rweb:> n_07 <- nrow(jd_2007) 
Rweb:> peak_2007 <- max(jd_2007$sie, na.rm=T) 
Rweb:> peak_dt_2007 <-   format(jd_2007$dt[which(jd_2007$sie== peak_2007)], "%m/%d") 
Rweb:> peak_doy_2007 <-  jd_2007$DOY[which(jd_2007$sie== peak_2007)] 
Rweb:> pit_2007 <- min(jd_2007$sie, na.rm=T) 
Rweb:> pit_dt_2007 <-   format(jd_2007$dt[which(jd_2007$sie== pit_2007)], "%m/%d") 
Rweb:> pit_doy_2007 <-  jd_2007$DOY[which(jd_2007$sie== pit_2007)] 
Rweb:> last_2007 <- jd_2007[n_07,7] 
Rweb:> decr_2007 <- peak_2007 - last_2007 
Rweb:>  
Rweb:> jd_2012 <- subset(jds, jds$yr == 2012) 
Rweb:> n_2012 <-  nrow(jd_2012) 
Rweb:> peak_2012 <- max(jd_2012$sie, na.rm=T) 
Rweb:> peak_dt_2012 <-   format(jd_2012$dt[which(jd_2012$sie== peak_2012)], "%m/%d") 
Rweb:> peak_doy_2012 <-  jd_2012$DOY[which(jd_2012$sie== peak_2012)] 
Rweb:> pit_2012 <- min(jd_2012$sie, na.rm=T) 
Rweb:> pit_dt_2012 <-   format(jd_2012$dt[which(jd_2012$sie== pit_2012)], "%m/%d") 
Rweb:> pit_doy_2012 <-  jd_2012$DOY[which(jd_2012$sie== pit_2012)] 
Rweb:> last_2012 <- jd_2012[n_2012,7] 
Rweb:> decr_2012 <- peak_2012 - last_2012 
Rweb:>  
Rweb:> jd_2016 <- subset(jds, jds$yr == 2016) 
Rweb:> n_2016 <-  nrow(jd_2016) 
Rweb:> peak_2016 <- max(jd_2016$sie, na.rm=T) 
Rweb:> peak_dt_2016 <-   format(jd_2016$dt[which(jd_2016$sie== peak_2016)], "%m/%d") 
Rweb:> peak_doy_2016 <-  jd_2016$DOY[which(jd_2016$sie== peak_2016)] 
Rweb:> pit_2016 <- min(jd_2016$sie, na.rm=T) 
Rweb:> pit_dt_2016 <-   format(jd_2016$dt[which(jd_2016$sie== pit_2016)], "%m/%d") 
Rweb:> pit_doy_2016 <-  jd_2016$DOY[which(jd_2016$sie== pit_2016)] 
Rweb:> last_2016 <- jd_2016[n_2016,7] 
Rweb:> decr_2016 <- peak_2016 - last_2016 
Rweb:>  
Rweb:> jd_1980 <- subset(jds, jds$yr == 1980) 
Rweb:> n_1980 <-  nrow(jd_1980) 
Rweb:> peak_1980 <- max(jd_1980$sie, na.rm=T) 
Rweb:> peak_dt_1980 <-   format(jd_1980$dt[which(jd_1980$sie== peak_1980)], "%m/%d") 
Rweb:> peak_doy_1980 <-  jd_1980$DOY[which(jd_1980$sie== peak_1980)] 
Rweb:> pit_1980 <- min(jd_1980$sie, na.rm=T) 
Rweb:> pit_dt_1980 <-   format(jd_1980$dt[which(jd_1980$sie== pit_1980)], "%m/%d") 
Rweb:> pit_doy_1980 <-  jd_1980$DOY[which(jd_1980$sie== pit_1980)] 
Rweb:> last_1980 <- jd_1980[n_1980,7] 
Rweb:> decr_1980 <- peak_1980 - last_1980 
Rweb:>  
Rweb:> jd_1990 <- subset(jds, jds$yr == 1990) 
Rweb:> n_1990 <-  nrow(jd_1990) 
Rweb:> peak_1990 <- max(jd_1990$sie, na.rm=T) 
Rweb:> peak_dt_1990 <-   format(jd_1990$dt[which(jd_1990$sie== peak_1990)], "%m/%d") 
Rweb:> peak_doy_1990 <-  jd_1990$DOY[which(jd_1990$sie== peak_1990)] 
Rweb:> pit_1990 <- min(jd_1990$sie, na.rm=T) 
Rweb:> pit_dt_1990 <-   format(jd_1990$dt[which(jd_1990$sie== pit_1990)], "%m/%d") 
Rweb:> pit_doy_1990 <-  jd_1990$DOY[which(jd_1990$sie== pit_1990)] 
Rweb:> last_1990 <- jd_1990[n_1990,7] 
Rweb:> decr_1990 <- peak_1990 - last_1990 
Rweb:>  
Rweb:> ###  Day of year axis setup 
Rweb:> ## Set up basic day of year vectors (mon_names, 1st day of mon) 
Rweb:> mon_names <- c("Jan", "Feb", "Mar", "April", "May", "June", "July", "Aug", "Sept", "Oct","Nov","Dec") 
Rweb:> mon_doy <- c(1,32,60,91,121,151,182,213,244,274,305,335,365) 
Rweb:> mon_pos <- c(16, 46, 75, 106,135, 165, 200, 228, 255, 289, 320, 355) 
Rweb:>  
Rweb:> y_min <- (min(jd_2010$sie) )*0.6 
Rweb:> y_max <- (max(jd_2010$sie))* 1.2 
Rweb:>  
Rweb:> # plot titles  
Rweb:> main_title <- paste("JAXA 2007, 2010, 2012, and 2016 Sea Ice Extent\nWith averages for decades 1980 and 1990\n",last_dt," Status", sep="") 
Rweb:> y_title <- "JAXA Arctic Sea Ice Extent - millions km^2" 
Rweb:>  
Rweb:> # Plot   
Rweb:> 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) 
Rweb:>  
Rweb:> # custom x & y axes 
Rweb:>   axis(side=1, at=mon_doy, labels=F, xaxs="i") 
Rweb:>   axis(side=1, at=mon_pos, labels=mon_names, tick=F, line=F, xaxs="i") 
Rweb:>   axis(side=2, at=NULL, labels=T, yaxs="i", las=1) 
Rweb:>   abline(v=1)    # need full length y axis line 
Rweb:>  
Rweb:> # Add 2007 trend & 2007-2010 peak sie 
Rweb:>   points(jd_2007$DOY, jd_2007$sie, type="l", col = "blue") 
Rweb:>   points(peak_doy_2010, peak_2010, type = "p", col = "red", pch=16, cex=1.25) 
Rweb:>   points(peak_doy_2007, peak_2007, type = "p", col = "blue", pch=16, cex=1.25) 
Rweb:>   points(pit_doy_2007, pit_2007, col = "blue", type="p", pch=18,cex=1.25) 
Rweb:>   points(pit_doy_2007, pit_2010, col = "red", type="p", pch=18, cex=1.25) 
Rweb:>  
Rweb:> # Add 2012 trend 
Rweb:>   points(jd_2012$DOY, jd_2012$sie, type="l", col = "green") 
Rweb:>   points(peak_doy_2012, peak_2012, type = "p", col = "green", pch=16, cex=1.25) 
Rweb:>   points(pit_doy_2012, pit_2012, col = "green", type="p", pch=18,cex=1.25) 
Rweb:>  
Rweb:> # Add 2016 trend 
Rweb:>   points(jd_2016$DOY, jd_2016$sie, type="l", col = "yellow") 
Rweb:>   points(peak_doy_2016, peak_2016, type = "p", col = "yellow", pch=16, cex=1.25) 
Rweb:>   points(pit_doy_2016, pit_2016, col = "yellow", type="p", pch=18,cex=1.25) 
Rweb:>  
Rweb:> # Add 1980 trend 
Rweb:>   points(jd_1980$DOY, jd_1980$sie, type="l", col = "black") 
Rweb:>   points(peak_doy_1980, peak_1980, type = "p", col = "black", pch=16, cex=1.25) 
Rweb:>   points(pit_doy_1980, pit_1980, col = "black", type="p", pch=18,cex=1.25) 
Rweb:>  
Rweb:> # Add 1990 trend 
Rweb:>   points(jd_1990$DOY, jd_1990$sie, type="l", col = "grey") 
Rweb:>   points(peak_doy_1990, peak_1990, type = "p", col = "grey", pch=16, cex=1.25) 
Rweb:>   points(pit_doy_1990, pit_1990, col = "grey", type="p", pch=18,cex=1.25) 
Rweb:>  
Rweb:> # create peak notes for legend   
Rweb:>   peak_note_2007 <- paste("2007 Peak @ ", signif(peak_2007, 4)," on ", peak_dt_2007, sep="") 
Rweb:>   peak_note_2010 <- paste("2010 Peak @ ", signif(peak_2010, 4)," on ", peak_dt_2010,  sep="") 
Rweb:>   peak_note_2012 <- paste("2012 Peak @ ", signif(peak_2012, 4)," on ", peak_dt_2012,  sep="") 
Rweb:>   peak_note_2016 <- paste("2016 Peak @ ", signif(peak_2016, 4)," on ", peak_dt_2016,  sep="") 
Rweb:>   peak_note_1980 <- paste("1980 Peak @ ", signif(peak_1980, 4)," on ", peak_dt_1980,  sep="") 
Rweb:>   peak_note_1990 <- paste("1990 Peak @ ", signif(peak_1990, 4)," on ", peak_dt_1990,  sep="") 
Rweb:>  
Rweb:> # create pit notes for legend   
Rweb:>   pit_note_2007 <- paste("2007 Pit @ ", signif(pit_2007, 4)," on ", pit_dt_2007, sep="") 
Rweb:>   pit_note_2010 <- paste("2010 Pit @ ", signif(pit_2010, 4)," on ", pit_dt_2010,  sep="") 
Rweb:>   pit_note_2012 <- paste("2012 Pit @ ", signif(pit_2012, 4)," on ", pit_dt_2012,  sep="") 
Rweb:>   pit_note_2016 <- paste("2016 Pit @ ", signif(pit_2016, 4)," on ", pit_dt_2016,  sep="") 
Rweb:>   pit_note_1980 <- paste("1980 Pit @ ", signif(pit_1980, 4)," on ", pit_dt_1980,  sep="") 
Rweb:>   pit_note_1990 <- paste("1990 Pit @ ", signif(pit_1990, 4)," on ", pit_dt_1990,  sep="") 
Rweb:>  
Rweb:> # legend 
Rweb:>   legend(40, 13, c( 
+ 	      peak_note_2007, peak_note_2010, peak_note_2012, peak_note_2016, peak_note_1980, peak_note_1990 
+ 	      ),  
+ 	      col = c("blue","red","green","yellow","black","grey"),  
+ 	      text.col = "black",  
+     	      lty=   c(0, 0, 0, 0, 0, 0), 
+               pch=   c(16,16,16,16,16,16), 
+ 	      pt.cex=c( 1, 1, 1, 1, 1, 1), 
+  	      merge = F, bg = "white", bty="o", box.col = "white", cex = .8) 
Rweb:>  
Rweb:>   legend(40, 8, c("2007 Trend", "2010 Trend", "2012 Trend", "2016 Trend", "1980 Average", "1990 Average" 
+ 	      ),  
+ 	      col = c("blue","red","green","yellow","black","grey"),  
+ 	      text.col = "black",  
+     	      lty=   c(1,1,1,1,1,1), 
+               pch=   c(0,0,0,0,0,0), 
+ 	      pt.cex=c(0,0,0,0,0,0), 
+  	      merge = F, bg = "white", bty="o", box.col = "white", cex = .8) 
Rweb:>  
Rweb:>   legend(120, 8, c( 
+     pit_note_2007, pit_note_2010, pit_note_2012, pit_note_2016, pit_note_1980, pit_note_1990 
+ 	      ),  
+ 	      col = c("blue","red","green","yellow","black","grey"),  
+ 	      text.col = "black",  
+     	      lty=   c(0, 0, 0, 0, 0, 0), 
+               pch=   c(18,18,18,18,18,18), 
+ 	      pt.cex=c(1, 1, 1, 1, 1, 1), 
+  	      merge = F, bg = "white", bty="o", box.col = "white", cex = .8) 
Rweb:>  
Rweb:>   delta <- last_2010-last_2007 
Rweb:>   delta_note <- paste(last_dt, " Delta: \n2010 - 2007 =\n ", signif(delta,3), " million km^2", sep="") 
Rweb:>   text(340, 16, delta_note, cex=0.75, adj=0.5) 
Rweb:>  
Rweb:> ## function to calc slope last x days 
Rweb:>  
Rweb:>   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) 
+    } 
Rweb:>  
Rweb:>   x <- c(2,3,5,7,10,20,30,60,90) 
Rweb:>   tr_df <- data.frame() 
Rweb:>   how_many <- length(x) 
Rweb:>   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) 
+     } 
Rweb:>  
Rweb:> names(tr_df) <- c("Prev days", "2007", "2010", "2010/2007") 
Rweb:>  
Rweb:> par(mar=c(0,0,0.5,0)) 
Rweb:>  
Rweb:> 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)) 
Rweb:>  
Rweb:>  
Rweb:> table_title <- paste(last_dt, ": Arctic SIE Decline Rates for Number of Previous Days", sep="") 
Rweb:> addtable2plot(0.825,2.95,tr_df,bty="o",display.rownames=F,hlines=T,cex=0.95, 
+   title=table_title, xjust=0, yjust=0) 
Rweb:>  
Rweb:>  


Gif Images