Data files used with the following script were downloaded from: ftp:\/\/ftp.ncdc.noaa.gov\/pub\/data\/ghcn\/daily\/.<\/p>\n
\r\n# Author: Manuel Gimond\r\n# Date: 10\/29\/2010\r\n#\r\n# Purpose: Extract data from DAILY GLOBAL HISTORICAL CLIMATOLOGY NETWORK for\r\n# a select time period. Only stations that satisfy the minimum NA count threshold\r\n# for the time period are kept.\r\n# The script converts the aggregated monthly data to z values, then generates trend\r\n# values for each station using linear regression model. The original inspiration\r\n# for this script was to reproduce fig.1 of Balling and Goodrich's "Spatial analysis\r\n# of variations in precipitation intensity in the USA" (Theor Appl Climatol, Oct 2010)\r\n# Libraries\r\n# library(MASS) # Used if robust regression is chosen\r\n#\r\n# User variables\r\nmy.path = "ghcnd_hcn" # Directory where raw files are stored\r\nfile.ext = ".dly" # Extension of filenames\r\nmet.par = "PRCP" # Meteorological paramater to extract\r\nmet.conv = 0.1 # Factor to convert raw number to real number (e.g. precip = 0.1)\r\nna.cap = 0.05 # Percent of no data accepted for analysis (stored as a fraction)\r\nyr.st = 1975 # Start year (inclusive)\r\nyr.en = 2009 # End year (inclusive)\r\nrobust.z = F # Use robust z scores? True (T) or False (F)\r\n#\r\n######################### Functions #########################################\r\n# FUNCTION: Computes standardized z values using conventional or robust methods\r\n# Note: all data values are converted including year\r\n# Table format should follow:\r\n# |Year,sta1,sta2,sta3,...|\r\nZval = function(x.in,Robust) {\r\nx.out = x.in # Copy input data frame to output data frame\r\nx.rows = nrow(x.in)\r\nx.cols = ncol(x.in)\r\n#\r\nif(Robust){ # quartile-based estimators\r\n m.x = apply(x.in, 2, median, na.rm = T)\r\n iqr.x = apply(x.in, 2, IQR, na.rm = T) # Normalized IQR\r\n s.x = iqr.x\/(qnorm(0.75) - qnorm(0.25))\r\n}else{ # moments-based estimators\r\n m.x = apply (x.in, 2, mean, na.rm=T)\r\n s.x = apply (x.in, 2, sd, na.rm=T)\r\n}\r\n#\r\nfor (i in 1:x.cols) {\r\n if ( s.x[i] != 0){\r\n x.out[,i] = (x.in [,i] - m.x[i])\/s.x[i] # z-score\r\n }else{\r\n x.out[,i] = 0\r\n }\r\n}\r\nreturn(x.out)\r\n}\r\n#\r\n# FUNCTION: Computes linear trend for each station across\r\n# entire time period. Appends coefficients to the dataframe.\r\n# Choice is between robust regression (Robust=T) or non-robust\r\n# regression (Robust=F).\r\n# Table format should follow:\r\n# |Year,sta1,sta2,sta3,...|\r\nStaReg = function(x,Robust) {\r\na = "intercept"\r\nb = "slope"\r\nR = "adj_R2" # Adjusted R2\r\nsta.lst = names(x)\r\nfor (i in 2:length(sta.lst)){\r\n if(Robust){\r\n lm.res = rlm(get(sta.lst[i]) ~ get(sta.lst[1]), data=x, na.action=na.omit)\r\n }else{\r\n lm.res = lm(get(sta.lst[i]) ~ get(sta.lst[1]), data=x, na.action=na.omit)\r\n }\r\n a[i] = lm.res$coefficients[1] # Intercept\r\n b[i] = lm.res$coefficients[2] # Slope\r\n R[i] = summary(lm.res)$adj.r.squared\r\n}\r\n# Append coefficients to existing table\r\nx = rbind(x,a,b,R)\r\nreturn(x)\r\n}\r\n######################### END of Functions #########################################\r\n#\r\n# Internal variables\r\nfile.stat = data.frame(name = numeric(0), expected = numeric(0),values = numeric(0),\r\n null = numeric(0), frac=numeric(0)) # Create an empty list to hold data stats\r\nfile.lst = list.files(path=my.path,pattern=file.ext) # Get list of files\r\nyr.days = (yr.en - yr.st +1) * 12 * 31 # Number of expected records for each station\r\n#\r\n# Internal variables specific to precipiation data\r\nall.Total = as.data.frame(seq(yr.st,yr.en,1))\r\ncolnames(all.Total) = "Year"\r\nall.INTNS = as.data.frame(seq(yr.st,yr.en,1))\r\ncolnames(all.INTNS) = "Year"\r\nall.N50 = as.data.frame(seq(yr.st,yr.en,1))\r\ncolnames(all.N50) = "Year"\r\n#\r\nfor (i in 1:length(file.lst)){\r\n#for (i in 1:30){\r\nprint( c("Processing",file.lst[i]))\r\nfile.stat[i,1] = file.lst[i]\r\n# Read data into columns (the '-1' value removes the flags)\r\nsub1 = read.fwf(paste(my.path,'\/',file.lst[i],sep=""),na.strings = "-9999",\r\n width=c(11,4,2,4,rep(c(5,-1,-1,-1),31)))\r\n# Grab only meteo parameter of interest\r\nmet.data = sub1[sub1$V4==met.par,]\r\n# Grab only years of interest\r\nmet.data = met.data[met.data$V2>=yr.st & met.data$V2<=yr.en,]\r\n# Compute the percent NA (make sure to ignore the header columns -> e.g. use seq(-1...))\r\n# Note that for months with fewer than 31 days, missing days are also counted as NA. This\r\n# should not influence the results significantly.\r\n#\r\nmet.data.na = sum(as.numeric(is.na(met.data[,seq(-1,-4,-1)]))) # Total NA\r\nmet.data.num = sum(as.numeric(abs(is.na(met.data[,seq(-1,-4,-1)])-1)))\r\n# Some files have missing dates; we need to add the missing data to NA\r\nmet.data.na = (yr.days - (met.data.num + met.data.na)) + met.data.na\r\nmet.data.na.frac = met.data.na \/ yr.days\r\n# Output file stats to list\r\nfile.stat[i,] = c(file.lst[i],yr.days,met.data.num,met.data.na,met.data.na.frac)\r\nif (met.data.na.frac <= na.cap) {\r\nmet.data.conv = cbind(met.data[2],met.data[3],met.data[,seq(-1,-4,-1)] * met.conv)\r\n * met.conv # Apply conversion factor to raw data & remove superfluous columns\r\n# The following is parameter specific (precipitation in this example)\r\n# Note the use of 'seq' (e.g. [,seq(-1,-2,-1)]) to exclude columns (e.g. 1 through 2)\r\nmet.data.conv = cbind(met.data[2],met.data[3],met.data[,seq(-1,-4,-1)] * met.conv)\r\nprec.cnt.mon = cbind(met.data.conv[1],met.data.conv[2],\r\n rowSums((met.data.conv[,seq(-1,-2,-1)] > 0),na.rm=T))\r\nprec.cnt.yr = aggregate(prec.cnt.mon[,seq(-1,-2,-1)], by = prec.cnt.mon[1],\r\n FUN = "sum") # total yearly count of rain events\r\nprec.tot.mon = cbind(met.data.conv[1],met.data.conv[2],rowSums(met.data.conv[,seq(-1,-2,-1)],na.rm=T))\r\nprec.tot.yr = aggregate(prec.tot.mon[,seq(-1,-2,-1)], by = prec.tot.mon[1], FUN = "sum")\r\nprec.INTNS.yr = cbind(prec.tot.yr[1], prec.tot.yr[2] \/ prec.cnt.yr[2])\r\nprec.N50.mon = cbind(met.data.conv[1],met.data.conv[2],\r\n rowSums((met.data.conv[,seq(-1,-2,-1)] >= 50),na.rm=T))\r\nprec.N50.yr = aggregate(prec.N50.mon[,seq(-1,-2,-1)], by = prec.N50.mon[1],\r\n FUN = "sum") # total yearly count of rain events above 50mm\r\nprec.Total.yr = prec.tot.yr\r\ncolnames(prec.Total.yr) = c("Year",gsub(file.ext,"",file.lst[i]))\r\ncolnames(prec.INTNS.yr) = c("Year",gsub(file.ext,"",file.lst[i]))\r\ncolnames(prec.N50.yr) = c("Year",gsub(file.ext,"",file.lst[i]))\r\nall.Total = merge(all.Total, prec.Total.yr, by="Year", all=T)\r\nall.INTNS = merge(all.INTNS, prec.INTNS.yr, by="Year", all=T)\r\nall.N50 = merge(all.N50, prec.N50.yr, by="Year", all=T)\r\n}\r\n}\r\n#\r\n# Convert rain indices to z-scores (normalized value) for each station (i.e. by column)\r\nall.Total.z = Zval(all.Total,robust.z)\r\nall.INTNS.z = Zval(all.INTNS,robust.z)\r\nall.N50.z = Zval(all.N50,robust.z)\r\n#\r\n# Compute temporal trend for each station\r\nall.Total.coef = StaReg(all.Total,F)\r\nall.Total.z.coef = StaReg(all.Total.z,F)\r\nall.INTNS.z.coef = StaReg(all.INTNS.z,F)\r\nall.N50.z.coef = StaReg(all.N50.z,F)\r\n#\r\n# Rename the first column (Year) to original value if it was converted to a z-score\r\nall.Total.z.coef[,1] = all.Total.coef[,1]\r\nall.INTNS.z.coef[,1] = all.Total.coef[,1]\r\nall.N50.z.coef[,1] = all.Total.coef[,1]\r\n#\r\n# Transpose new data subset so that the stations are listed as records\r\nall.Total.coef.t = t(all.Total.coef)\r\nall.Total.z.coef.t = t(all.Total.z.coef)\r\nall.INTNS.t = t(all.INTNS)\r\nall.INTNS.z.coef.t = t(all.INTNS.z.coef)\r\nall.N50.t = t(all.N50)\r\nall.N50.z.coef.t = t(all.N50.z.coef)\r\n#\r\n# Rename header name and remove redundant header\r\ncolnames(all.Total.coef.t) = all.Total.coef.t[1,]\r\ncolnames(all.Total.z.coef.t) = all.Total.z.coef.t[1,]\r\ncolnames(all.INTNS.t) = all.INTNS.t[1,]\r\ncolnames(all.INTNS.z.coef.t) = all.INTNS.z.coef.t[1,]\r\ncolnames(all.N50.t) = all.N50.t[1,]\r\ncolnames(all.N50.z.coef.t) = all.N50.z.coef.t[1,]\r\nall.Total.coef.t = data.frame(all.Total.coef.t)[-1,]\r\nall.Total.z.coef.t = data.frame(all.Total.z.coef.t)[-1,]\r\nall.INTNS.t = data.frame(all.INTNS.t)[-1,]\r\nall.INTNS.z.coef.t = data.frame(all.INTNS.z.coef.t)[-1,]\r\nall.N50.t = data.frame(all.N50.t)[-1,]\r\nall.N50.z.coef.t = data.frame(all.N50.z.coef.t)[-1,]\r\n#\r\n# Export to CSV file\r\nwrite.table(all.Total.coef.t,file="precip_Total_coef.csv", sep=",",\r\n quote=F,col.names=NA,row.names=T,na="") # Set NA to blank\r\nwrite.table(all.Total.z.coef.t,file="precip_Total_z.csv", sep=",",\r\n quote=F,col.names=NA,row.names=T,na="") # Set NA to blank\r\n<\/pre>\n
\n<\/code><\/p>\nThe following figure shows the standardized regression coefficients from a Universal Kriging analysis of the Precip_total_coef_all_z.csv file (the values are in fact correlation coefficients). The time period spans 1975 – 2009.<\/p>\n
![\"\"](\"http:\/\/web.colby.edu\/mgimond\/files\/2010\/11\/1975_2007_Total_precip_z_regr.jpg\")
<\/a><\/p>\n","protected":false},"excerpt":{"rendered":"