############## RClimate Script: AMU Channel 5 Temperature Anomaly ########################## ## Script stored on http://chartsgraphs.wordpress.com account for Users to source() ## ## ## Download and process UAH's channel 5 data ## ## http://discover.itsc.uah.edu/amsutemps/data/amsu_daily_85N85S_ch05.r002.txt ## ## ## Developed by D Kelly O'Day to demonstrate use of source() function for climate data ## ## http:chartsgraphs.wordpress.com 5/9/10 ## ############################################################################################ library(ggplot2); library(zoo) ## Read source data file from UAH uah_link <- "http://discover.itsc.uah.edu/amsutemps/data/amsu_daily_85N85S_ch05.r002.txt" ch5 <- read.table(uah_link, skip = 27, sep = "", dec=".", as.is = T, row.names = NULL, header = F, colClasses=c("character", rep("numeric",16)), na.strings= c("-999.000","0.000"), col.names = c("doy","y_1998", "y_1999", "y_2000", "y_2001", "y_2002", "y_2003", "y_2004", "y_2005", "y_2006", "y_2007", "y_2008", "y_2009", "y_2010", "20_avg", "20_max", "20_min")) ## Build long vectors num_yrs <- 2010-1997 num_rows <- 366*num_yrs # Establish new vectors dc <- character(num_rows) # date_char t <- numeric(num_rows) # temp C avg <- numeric(num_rows) # 20 yr avg a <- numeric(num_rows) # anom - C based on 20 year avg n <- 0 # row counter for long vectors dmax_a <- numeric(366) # Build date and temp vectors with 2 loops # Outer loop goes across 14 data columns- years # Inner loop goes down 366 days/yr for(c in 2:14) { for (d in 1:366) { n<- n+1 # long vector row count dc[n] <- paste(ch5[d,1],"/", 1996+c, sep="") # date in char mode t[n] <- ch5[d,c] avg[n] <- ch5[d,15] # read 20 yr avg for day a[n] <- t[n] - avg[n] # calc anomaly } } dmax_a[1:366] <- ch5[,16]-ch5[,15] # max anom by day of year dd <- as.Date(dc, format = "%m/%d/%Y") # date in date mode yr <- as.numeric(format(dd, format="%Y")) mo <- as.numeric(format(dd, format="%m")) dy <- as.numeric(format(dd, format="%d")) yr_frac <- as.numeric(yr + (mo-1)/12 + (dy/30)/12) ## Consolidate long vectors into data.frame df <- data.frame(dd,yr_frac, t,avg,a,dmax_a) #Use ggplots's remove_missing() to eliminate rows with missing data in data.frame cdf <- remove_missing(df, na.rm = T) ## Determine month - yr of last reading nr <- nrow(cdf ) # Find number of data rows l_dd <- cdf$dd[nr] # Find date of last data point l_mon <- as.numeric(format(l_dd, format="%m")) l_day <- as.numeric(format(l_dd, format="%d")) l_mon_day <- paste(l_mon, "/",l_day, sep="") l_yr <- as.numeric(format(l_dd, format="%Y")) l_a <- signif(cdf$a[nr],5 ) # Find last anom l_yr_frac <- cdf$yr_frac[nr] l_m_n <- months(l_dd, abbreviate = F) l_m_num <- as.numeric(l_m_n) ## create zoo object to be able to use zoo() aggregate() functions ch5_daily <- with(cdf, zoo(a, dd)) ch5_mon <- aggregate(ch5_daily, as.yearmon, mean, na.rm = T) m <- as.numeric(format(as.Date(index(ch5_mon)), "%m")) m_df <- data.frame(m, ch5_mon) # data frame of monthly average anomalies ## Subset current month (cm) cm_df <- subset(m_df, m_df$m == l_mon) cm_min <- min(cm_df$ch5_mon) cm_avg <- mean(cm_df$ch5_mon) cm_max <- max(cm_df$ch5_mon, na.rm=T) ## subset current month first_day <- as.Date(paste(l_yr,l_mon,"1",sep="-")) d_in_m <- c(31,28,31,30,31,30,31,31,30,31,30,31) last_day <-as.Date(first_day + d_in_m[l_mon]-1) l_m_df <- subset(df, df$dd>= first_day & df$dd <= last_day) min_dd_lm <- min(l_m_df$dd, na.rm=T) max_dd_lm <- max(l_m_df$dd, na.rm=T) min_a_lm <- min(l_m_df$dmax_a, na.rm=T) max_a_lm <- max(l_m_df$dmax_a, na.rm=T) avg_lm <- mean(l_m_df$a, na.rm=T) ### Plot Notes l_dd_note <- paste(l_dd, " @ ",l_a, sep="") title <- paste("UAH Channel 5 Anomaly Trends: (14,000 ft/ 4.4 km/ 600mb)\n a) Monthly Avg (1998 - 2010) and\n b) Daily: ",l_m_n, ",", l_yr," through ",l_mon_day, sep="") fig_b_y_lab <- paste(format(first_day, "%b"), " Anomaly - C", sep="") max_mon_note <- paste("Previous max ", l_m_n, " @ ", signif(cm_max,3), sep="") avg_lm_note <- paste("Avg ", l_m_n, " MTD - ", signif(avg_lm,3), sep="") max_daily <- paste(format(first_day, "%b"), " Max Anomaly by Day", sep="") ## Plot ch 5 anomaly series out_link <- "C:\\R_Home\\Charts & Graphs Blog\\RClimateTools\\UAH_ch05\\UAH_ch5_trends_long_recent.png" # png(file=out_link,width=600,height=550, res=72) par(las=1);par(ps=11); par(las=1); par(oma=c(2.5,1,1.5,1)); par(mar=c(2,4,2.5,0)) plot(df$yr_frac, df$a, type = "n", xlab = "", ylab = "Temperature Anomaly - C (20 yr baseline period)", col = "grey", ylim=c(-1,1),main=title, cex.main=0.8 ) points(ch5_mon,col="blue", type="s") points(l_yr_frac, l_a, type = "p",pch=16, col = "red") text(1999,0.95,"a) Monthly Avg Trend", adj=0) ############################################################################################ # Generate and add bottom footer KOD , source, System date notes source_note <- paste("Data Source: ", uah_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) ############################################################################################# legend(1997.7,-0.5, c(l_dd_note, "Monthly Anomaly", avg_lm_note, max_daily, max_mon_note), col = c("red", "blue","green", "brown", "orange"), text.col = "grey", lty = c(0,1,1,1,1),pch=c(16,0,0,0,0),pt.cex=c(1,0,0,0,0), merge = F, bg = "white", bty="o", box.col="lightgrey", cex = .8) ### Create 2nd chart showing current month par(fig=c(0.55,.97,0.1, 0.525), new = T, las=1, par(ps=11)) plot(l_m_df$dd,l_m_df$a, type = "n", col = "red", xlab = "", ylim=c(min_a_lm-0.025,1),xlim=c(min_dd_lm-1, max_dd_lm+1), ylab=fig_b_y_lab, bty="o", xaxs = "i", yaxs = "i", axes = T, mgp=c(2.3,1,0) ) box(col = "darkgrey") abline(h=avg_lm, col="green") points(l_m_df$dd, l_m_df$dmax_a, type="l", col = "brown") abline(h=cm_max, col="orange") points(l_m_df$dd,l_m_df$a, type = "b", col = "red",pch=16, cex=0.3) points(l_dd, l_a,type="p",col="red", pch=16) text(min_dd_lm, 0.95, paste("b) ", l_m_n, " Daily Values",sep=""), adj=0)