########################################################################### ## RClimate.R ## ## Collection of Climate data analysis functions to assist with ## ## climate data analysis. These functions developed and maintained by ## ## D Kelly O'Day, http://chartsgraphs.wordpress.com ## ## Data files and RClimate scripts stored at http://processtrends.com ## ## Users are asked to credit RClimate when using scripts or plots ## ## June 14, 2010 ## ## Updated 8/31/10 to include get_lota() and which_series(num) ## ## auto update system ## ########################################################################### #### yr_mon function ############################################## ## to get yr_mon from yr_frac # y_f is yr_frac vector func_yr_mon <- function(y_f) { yr <- as.integer(y_f) ## Each month is 1/12 or 0.083 of calandar year inc <- 1/12 mo <- ceiling((y_f-yr)/(inc)) mo_char <- formatC(mo,width=2,flag='0') yr_mon <- as.numeric(as.character(paste(yr, mo_char, sep="") )) return(yr_mon) } ##### GISS function ################################################## # GISS Monthly Temp Anomaly Trend - Land & Sea C - 1880 - latest month # Function to retrieve GISS monthly LOTA produce long format data.frame func_GISS <- function() { library(ggplot2) url <- c("http://data.giss.nasa.gov/gistemp/tabledata/GLB.Ts+dSST.txt") file <- c("GLB.Ts+dSST.txt") download.file(url, file) ## Cleanup File rows <- length(readLines(file)) - 12 # Read file as char vector, one line per row, Exclude first 7 rows lines <- readLines(file, n=rows)[9:rows] lines2 <- gsub("\\*{3,5}", " NA", lines, perl=TRUE) #Convert the character vector to a dataframe df <- read.table( textConnection(lines2), header=TRUE, colClasses = "character") closeAllConnections() # We are only interested in the montly data in first 13 columns df <- df[,1:13] # Convert all variables (columns) to numeric format df <- colwise(as.numeric) (df) # Remove rows where Year=NA from the dataframe df <- df [!is.na(df$Year),] # Convert from wide format to long format dfm <- melt(df, id.var="Year", variable_name="Month") GISS_mo_num <- unclass(dfm$Month) GISS_mo_frac <- as.numeric(( unclass(dfm$Month)-0.5)/12) GISS_yr_frac <- dfm$Year + GISS_mo_frac GISS_anom <- dfm$value/100 GISS_yr_mon <- func_yr_mon(GISS_yr_frac) dfm <- data.frame(dfm, GISS_mo_num, GISS_yr_frac, GISS_yr_mon, GISS_anom) dfm <- dfm[order(dfm$GISS_yr_frac), ] dfm <- dfm[!is.na(dfm$GISS_anom),] names(dfm) <- c("year", "month", "value","mo_num", "yr_frac", "yr_mon", "GISS") dfm<- dfm[,c(5,6,7)] return(dfm) } ############################################################### ## Download and process UAH Data File ############### func_UAH <- function() { UAH_link <- "http://vortex.nsstc.uah.edu/public/msu/t2lt/tltglhmam_5.3" UAH_file <- c("tltglhmam_5.3") download.file(UAH_link, UAH_file) library(reshape) ## The first 5 rows and the last row of the textfile contain text information, not data ## Find out the number of rows in a file, and exclude the last 1,it has text UAH_rows <- length(readLines(UAH_file)) - 1 ## Read file as char vector, one line per row, Exclude first 5 rows UAH_lines <- readLines(UAH_file, n=UAH_rows)[6:UAH_rows] ## Data Manipulation, R vector with data lines. ## Use regexp to replace all the occurences of **** with NA lines2 <- gsub("\\*{3,5}", " NA", UAH_lines, perl=TRUE) ##Convert the character vector to a dataframe UAH_df <- read.table( textConnection(lines2), header=F, colClasses = "character") closeAllConnections() ## We are only interested in the most recent montly data columns 1:3 UAH_df <- UAH_df[,1:3] ## Convert all variables (columns) to numeric format UAH_df <- colwise(as.numeric) (UAH_df) UAH_yr_frac <- UAH_df[1]+ +(UAH_df[2]-0.5)/12 UAH_yr_mon <- func_yr_mon(UAH_yr_frac$V1) UAH_df <- data.frame(UAH_yr_frac, UAH_yr_mon, UAH_df) #UAH_num <- nrow(UAH_df) # UAH_yf <- UAH_df[UAH_num,1] UAH_df <- UAH_df[,c(1,2,5)] names(UAH_df) <- c("yr_frac", "yr_mon", "UAH") UAH_df <- subset(UAH_df, UAH_yr_frac >= 1979) return(UAH_df) } ######################################################################################### ################ NOAA ################################################################### func_NOAA <- function() { #link_NOAA <- "ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.land_and_ocean.90S.90N.df_1901-2000mean.dat" link_NOAA <- "ftp://ftp.ncdc.noaa.gov/pub/data/anomalies/monthly.land_ocean.90S.90N.df_1901-2000mean.dat" NOAA_in <- read.table(link_NOAA, skip = 0, sep = "", dec=".", row.names = NULL, header = FALSE, as.is = T, colClasses = c(rep("numeric",3)), comment.char = "#", na.strings = c("*", "-",-99.9, -999.9), col.names = c("NOAA_yr", "NOAA_mo", "NOAA")) yr_frac <- NOAA_in$NOAA_yr + (NOAA_in$NOAA_mo-0.5)/12 # yr_frac yr_mon <- func_yr_mon(yr_frac) NOAA_in <- data.frame(yr_frac, yr_mon, NOAA_in$NOAA) NOAA_df <- subset(NOAA_in, NOAA_in$NOAA > -999) #NOAA_num <- nrow(NOAA_df) #NOAA_yf <- NOAA_df[NOAA_num,1] names(NOAA_df) <- c("yr_frac", "yr_mon", "NOAA") return(NOAA_df) } ########################################################################################## ############################################################################################ ########################## HADCRUT3 ######################################################### ## Download and process HADCRUT3 Data File ############### func_HAD <- function() { url_HAD <- c("http://www.cru.uea.ac.uk/cru/data/temperature/hadcrut3gl.txt") file <- c("hadcrut3gl.txt") download.file(url_HAD, file) ## Find out the number of rows in a file HAD_rows <- length(readLines(file)) ## HAD file has alternating data & observation count for each year ## Want to delete the 2nd row for year; Develop vector for rows tobe deleted HAD_pts <- seq(from = 1, to = HAD_rows, by=2) ## Read HAD data file as char vector, one line per row HAD_lines <- readLines(file, n=HAD_rows) ## Create wanted data by applying HAD_pts vector as index vector HAD_liones HAD_data <- HAD_lines[HAD_pts] ##Convert the HAD_data character vector to a dataframe HAD_wide <- read.table( textConnection(HAD_data), header=F, colClasses = "character") closeAllConnections() ## Convert HAD_df from text data.frame to numeric data.frame HAD_wide <- colwise(as.numeric) (HAD_wide) ## We are only interested in the montly data in first 13 columns; drop 1th col, annual value HAD_wide <- HAD_wide[,1:13] names(HAD_wide) <- c("year", "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec") ## Convert HAD_wide from wide format to long format with melt() HAD_dfm <- melt(HAD_wide, id.var="year", variable_name="Month") names(HAD_dfm) <- c("year","Month", "HAD") ## Establish month numbers for yr_frac calculation HAD_mo_num <- as.numeric(unclass(HAD_dfm$Month)) ## Create yr_frac vector yr_frac <- as.numeric(HAD_dfm$year) + ((HAD_mo_num-0.5)/12) yr_mon <- func_yr_mon(yr_frac) ## Create HAD_dfm data.frame & order by yr_frac HAD_df <- data.frame(yr_frac, yr_mon, HAD_dfm$HAD) HAD_df <- HAD_df[order(HAD_df$yr_frac), ] ## Remove all HAD_anom == 0.000 HAD_df <- subset(HAD_df, HAD_df$yr_frac >1880 & HAD_df$HAD != 0.000) names(HAD_df) <- c("yr_frac", "yr_mon", "HAD") return(HAD_df) } ################### RSS Data ################################################################ func_RSS <- function() { link_RSS<- "http://www.remss.com/data/msu/monthly_time_series/RSS_Monthly_MSU_AMSU_Channel_TLT_Anomalies_Land_and_Ocean_v03_2.txt" mo_RSS_in <- read.table(link_RSS, skip = 3, sep = "", dec=".", row.names = NULL, header = FALSE, as.is = T, colClasses = c(rep("numeric",11)), comment.char = "#", na.strings = c("*", "-",-99.9, -999.9), col.names = c("yr", "mo", "RSS", "s20_n20", "n20_n825", "s70_s20", "n60_n825", "R", "USA", "e_n825", "s70_e")) yr_frac <- mo_RSS_in$yr + (mo_RSS_in$mo-0.5)/12 # yr_frac simplifies calcs yr_mon <- func_yr_mon(yr_frac) ## Create new RSS_df data frame RSS_df <- data.frame(yr_frac,yr_mon, mo_RSS_in[3]) return(RSS_df) } ############################################################ get_lota <- function() { link <- "http://processtrends.com/Files/RClimate_consol_temp_anom_latest.csv" lota <- read.csv(link) return(lota) } ####################################################################################### which_series <- function(num) { st <- c(rep(1880,3),rep(1979,2)) series_names<- c("GISS", "HAD", "NOAA", "RSS", "UAH") d_col = num[1] + 2 # define lota column to plot series_id <- series_names[num] series_data <- lota[,c(1,d_col)] cat("\n", "You entered", series_id, "\n") st_yr <- st[num] series_data <- subset(series_data, series_data$yr_frac> st_yr) series_ts <- ts(series_data,start=st_yr, freq=12) ma_13 <- filter(series_ts, rep(1/13,13), sides=1) main_title <- paste(series_id, " Global Temperature Anomaly Trend - C", sep="") plot(series_data[,1], series_data[,2], type = "l", col = "grey", xlab="", las=1, main = main_title, ylab = "Temperature Anomaly - C" ) points(ma_13[,1],ma_13[,2], col = "red", type = "l") legend(1885,.6, c("Monthly Avg", "13 Month Moving avg"), col = c("grey", "red"), text.col = "black", lty = c(1,1), merge = F, bg = "white", bty="o", box.col = "grey", cex = .7) } #CALCULATE ANNUAL AVERAGE FROM MONTHLY TIME SERIES annavg<-function(x) { n<-length(x) m<-n-12*floor(n/12) if(m>0) x<-c(x, rep(NA,12-m)) years<-length(x)/12 x<-array(x,dim=c(12,years)) annavg<-apply(x,2,mean,na.rm=T) return(annavg) }