# Sean Demet # Department of Physical Geography and Ecosystem Analysis # Lund University # Masters Degree Thesis # Advisor: Veiko Lehsten # This analysis was performed using the R version 2.15.1 (2012-06-22) setwd("/Users/SeanGDemet/Data/Regions") library(xtable) ##### # Here I import the data for Model Development and Model Evaluation. # NOTE: The .txt files have been created using CDO. The Lon/Lat have been corrected to match, the regions have been labeled, and the 3hr meteorological variables have been changed to 24hr variables. ##### # Reading in .txt data for REGION 1 development and evaluation ##### read.table("y2003_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r1 read.table("y2004_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r1 read.table("y2005_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r1 read.table("y2006_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r1 read.table("y2007_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r1 # Here I rewrite day values be continuous y04r1[,"day"] <- y04r1[,"day"]+max(y03r1$day) y05r1[,"day"] <- y05r1[,"day"]+max(y04r1$day) y06r1[,"day"] <- y06r1[,"day"]+max(y05r1$day) y07r1[,"day"] <- y07r1[,"day"]+max(y06r1$day) # Creating the Region DEVreg1 <- rbind(y03r1,y04r1,y05r1,y06r1,y07r1) # Reading in data for model evaluation read.table("y2008_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r1 read.table("y2009_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r1 read.table("y2010_r1.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r1 # Here I rewrite day values be continuous y08r1[,"day"] <- y08r1[,"day"]+max(y07r1$day) y09r1[,"day"] <- y09r1[,"day"]+max(y08r1$day) y10r1[,"day"] <- y10r1[,"day"]+max(y09r1$day) # Creating the Region EVALreg1 <- rbind(y08r1,y09r1,y10r1) ##### # Reading in .txt data for REGION 2 development and evaluation ##### read.table("y2003_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r2 read.table("y2004_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r2 read.table("y2005_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r2 read.table("y2006_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r2 read.table("y2007_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r2 # Here I rewrite day values be continuous y04r2[,"day"] <- y04r2[,"day"]+max(y03r2$day) y05r2[,"day"] <- y05r2[,"day"]+max(y04r2$day) y06r2[,"day"] <- y06r2[,"day"]+max(y05r2$day) y07r2[,"day"] <- y07r2[,"day"]+max(y06r2$day) # Creating the Region DEVreg2 <- rbind(y03r2,y04r2,y05r2,y06r2,y07r2) # Reading in data for model evaluation read.table("y2008_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r2 read.table("y2009_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r2 read.table("y2010_r2.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r2 # Here I rewrite day values be continuous y08r2[,"day"] <- y08r2[,"day"]+max(y07r2$day) y09r2[,"day"] <- y09r2[,"day"]+max(y08r2$day) y10r2[,"day"] <- y10r2[,"day"]+max(y09r2$day) # Creating the Region EVALreg2 <- rbind(y08r2,y09r2,y10r2) ##### # Reading in .txt data for REGION 4 development and evaluation ##### read.table("y2003_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r4 read.table("y2004_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r4 read.table("y2005_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r4 read.table("y2006_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r4 read.table("y2007_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r4 # Here I rewrite day values be continuous y04r4[,"day"] <- y04r4[,"day"]+max(y03r4$day) y05r4[,"day"] <- y05r4[,"day"]+max(y04r4$day) y06r4[,"day"] <- y06r4[,"day"]+max(y05r4$day) y07r4[,"day"] <- y07r4[,"day"]+max(y06r4$day) # Creating the Region DEVreg4 <- rbind(y03r4,y04r4,y05r4,y06r4,y07r4) # Reading in data for model evaluation read.table("y2008_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r4 read.table("y2009_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r4 read.table("y2010_r4.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r4 # Here I rewrite day values be continuous y08r4[,"day"] <- y08r4[,"day"]+max(y07r4$day) y09r4[,"day"] <- y09r4[,"day"]+max(y08r4$day) y10r4[,"day"] <- y10r4[,"day"]+max(y09r4$day) # Creating the Region EVALreg4 <- rbind(y08r4,y09r4,y10r4) ##### # Reading in .txt data for REGION 12 development and evaluation ##### read.table("y2003_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r12 read.table("y2004_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r12 read.table("y2005_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r12 read.table("y2006_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r12 read.table("y2007_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r12 # Here I rewrite day values be continuous y04r12[,"day"] <- y04r12[,"day"]+max(y03r12$day) y05r12[,"day"] <- y05r12[,"day"]+max(y04r12$day) y06r12[,"day"] <- y06r12[,"day"]+max(y05r12$day) y07r12[,"day"] <- y07r12[,"day"]+max(y06r12$day) # Creating the Region DEVreg12 <- rbind(y03r12,y04r12,y05r12,y06r12,y07r12) # Reading in data for model evaluation read.table("y2008_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r12 read.table("y2009_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r12 read.table("y2010_r12.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r12 # Here I rewrite day values be continuous y08r12[,"day"] <- y08r12[,"day"]+max(y07r12$day) y09r12[,"day"] <- y09r12[,"day"]+max(y08r12$day) y10r12[,"day"] <- y10r12[,"day"]+max(y09r12$day) # Creating the Region EVALreg12 <- rbind(y08r12,y09r12,y10r12) ##### # Reading in .txt data for REGION 13 development and evaluation ##### read.table("y2003_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r13 read.table("y2004_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r13 read.table("y2005_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r13 read.table("y2006_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r13 read.table("y2007_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r13 # Here I rewrite day values be continuous y04r13[,"day"] <- y04r13[,"day"]+max(y03r13$day) y05r13[,"day"] <- y05r13[,"day"]+max(y04r13$day) y06r13[,"day"] <- y06r13[,"day"]+max(y05r13$day) y07r13[,"day"] <- y07r13[,"day"]+max(y06r13$day) # Creating the Region DEVreg13 <- rbind(y03r13,y04r13,y05r13,y06r13,y07r13) # Reading in data for model evaluation read.table("y2008_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r13 read.table("y2009_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r13 read.table("y2010_r13.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r13 # Here I rewrite day values be continuous y08r13[,"day"] <- y08r13[,"day"]+max(y07r13$day) y09r13[,"day"] <- y09r13[,"day"]+max(y08r13$day) y10r13[,"day"] <- y10r13[,"day"]+max(y09r13$day) # Creating the Region EVALreg13 <- rbind(y08r13,y09r13,y10r13) ##### # Reading in .txt data for REGION 14 development and evaluation ##### read.table("y2003_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r14 read.table("y2004_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r14 read.table("y2005_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r14 read.table("y2006_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r14 read.table("y2007_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r14 # Here I rewrite day values be continuous y04r14[,"day"] <- y04r14[,"day"]+max(y03r14$day) y05r14[,"day"] <- y05r14[,"day"]+max(y04r14$day) y06r14[,"day"] <- y06r14[,"day"]+max(y05r14$day) y07r14[,"day"] <- y07r14[,"day"]+max(y06r14$day) # Creating the Region DEVreg14 <- rbind(y03r14,y04r14,y05r14,y06r14,y07r14) # Reading in data for model evaluation read.table("y2008_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r14 read.table("y2009_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r14 read.table("y2010_r14.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r14 # Here I rewrite day values be continuous y08r14[,"day"] <- y08r14[,"day"]+max(y07r14$day) y09r14[,"day"] <- y09r14[,"day"]+max(y08r14$day) y10r14[,"day"] <- y10r14[,"day"]+max(y09r14$day) # Creating the Region EVALreg14 <- rbind(y08r14,y09r14,y10r14) ##### # Reading in .txt data for REGION 18 development and evaluation ##### read.table("y2003_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r18 read.table("y2004_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r18 read.table("y2005_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r18 read.table("y2006_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r18 read.table("y2007_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r18 # Here I rewrite day values be continuous y04r18[,"day"] <- y04r18[,"day"]+max(y03r18$day) y05r18[,"day"] <- y05r18[,"day"]+max(y04r18$day) y06r18[,"day"] <- y06r18[,"day"]+max(y05r18$day) y07r18[,"day"] <- y07r18[,"day"]+max(y06r18$day) # Creating the Region DEVreg18 <- rbind(y03r18,y04r18,y05r18,y06r18,y07r18) # Reading in data for model evaluation read.table("y2008_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r18 read.table("y2009_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r18 read.table("y2010_r18.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r18 # Here I rewrite day values be continuous y08r18[,"day"] <- y08r18[,"day"]+max(y07r18$day) y09r18[,"day"] <- y09r18[,"day"]+max(y08r18$day) y10r18[,"day"] <- y10r18[,"day"]+max(y09r18$day) # Creating the Region EVALreg18 <- rbind(y08r18,y09r18,y10r18) ##### # Reading in .txt data for REGION 19 development and evaluation ##### read.table("y2003_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r19 read.table("y2004_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r19 read.table("y2005_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r19 read.table("y2006_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r19 read.table("y2007_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r19 # Here I rewrite day values be continuous y04r19[,"day"] <- y04r19[,"day"]+max(y03r19$day) y05r19[,"day"] <- y05r19[,"day"]+max(y04r19$day) y06r19[,"day"] <- y06r19[,"day"]+max(y05r19$day) y07r19[,"day"] <- y07r19[,"day"]+max(y06r19$day) # Creating the Region DEVreg19 <- rbind(y03r19,y04r19,y05r19,y06r19,y07r19) # Reading in data for model evaluation read.table("y2008_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r19 read.table("y2009_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r19 read.table("y2010_r19.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r19 # Here I rewrite day values be continuous y08r19[,"day"] <- y08r19[,"day"]+max(y07r19$day) y09r19[,"day"] <- y09r19[,"day"]+max(y08r19$day) y10r19[,"day"] <- y10r19[,"day"]+max(y09r19$day) # Creating the Region EVALreg19 <- rbind(y08r19,y09r19,y10r19) ##### # Reading in .txt data for REGION 20 development and evaluation ##### read.table("y2003_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y03r20 read.table("y2004_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y04r20 read.table("y2005_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y05r20 read.table("y2006_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y06r20 read.table("y2007_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y07r20 # Here I rewrite day values be continuous y04r20[,"day"] <- y04r20[,"day"]+max(y03r20$day) y05r20[,"day"] <- y05r20[,"day"]+max(y04r20$day) y06r20[,"day"] <- y06r20[,"day"]+max(y05r20$day) y07r20[,"day"] <- y07r20[,"day"]+max(y06r20$day) # Creating the Region DEVreg20 <- rbind(y03r20,y04r20,y05r20,y06r20,y07r20) # Reading in data for model evaluation read.table("y2008_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y08r20 read.table("y2009_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y09r20 read.table("y2010_r20.txt",col.names=c("lon","lat","day","FRP","FRP0","year","swc","wind","lsp","cp","d","t","region")) -> y10r20 # Here I rewrite day values be continuous y08r20[,"day"] <- y08r20[,"day"]+max(y07r20$day) y09r20[,"day"] <- y09r20[,"day"]+max(y08r20$day) y10r20[,"day"] <- y10r20[,"day"]+max(y09r20$day) # Creating the Region EVALreg20 <- rbind(y08r20,y09r20,y10r20) ##### # Here I select FRP events from the datasets, add additional data (FRP0, FRPratio, convert temperature to Celsius, total precip, Angstrom Index, Relative Humidity, tp0, lp0, & cp0) # I perform this for each region and save them as .RData to be called later ##### # REGION 1 (DEVreg1 -> r01_DEV.RData) ##### region <- DEVreg1 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg1 save(FRPreg1, file="r01_DEV.RData") ##### # REGION 1 (EVALreg1 -> r01_EVAL.RData) ##### region <- EVALreg1 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg1e save(FRPreg1e, file="r01_EVAL.RData") ##### # REGION 2 (DEVreg2 -> r02_DEV.RData) ##### region <- DEVreg2 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg2 save(FRPreg2, file="r02_DEV.RData") ##### # REGION 2 (EVALreg2 -> r02_EVAL.RData) ##### region <- EVALreg2 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg2e save(FRPreg2e, file="r02_EVAL.RData") ##### # REGION 4 (DEVreg4 -> r04_DEV.RData) ##### region <- DEVreg4 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg4 save(FRPreg4, file="r04_DEV.RData") ##### # REGION 4 (EVALreg4 -> r04_EVAL.RData) ##### region <- EVALreg4 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg4e save(FRPreg4e, file="r04_EVAL.RData") ##### # REGION 12 (DEVreg12 -> r12_DEV.RData) ##### region <- DEVreg12 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP events FRPreg <- region[which(region$FRP>0),] #- Calculate FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg12 save(FRPreg12, file="r12_DEV.RData") ##### # REGION 12 (EVALreg12 -> r12_EVAL.RData) ##### region <- EVALreg12 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg12e save(FRPreg12e, file="r12_EVAL.RData") ##### # REGION 13 (DEVreg13 -> r13_DEV.RData) ##### region <- DEVreg13 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg13 save(FRPreg13, file="r13_DEV.RData") ##### # REGION 13 (EVALreg13 -> r13_EVAL.RData) ##### region <- EVALreg13 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg13e save(FRPreg13e, file="r13_EVAL.RData") ##### # REGION 14 (DEVreg14 -> r14_DEV.RData) ##### region <- DEVreg14 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg14 save(FRPreg14, file="r14_DEV.RData") ##### # REGION 14 (EVALreg14 -> r14_EVAL.RData) ##### region <- EVALreg14 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg14e save(FRPreg14e, file="r14_EVAL.RData") ##### # REGION 18 (DEVreg18 -> r18_DEV.RData) ##### region <- DEVreg18 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg18 save(FRPreg18, file="r18_DEV.RData") ##### # REGION 18 (EVALreg18 -> r18_EVAL.RData) ##### region <- EVALreg18 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg18e save(FRPreg18e, file="r18_EVAL.RData") ##### # REGION 19 (DEVreg19 -> r19_DEV.RData) ##### region <- DEVreg19 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg19 save(FRPreg19, file="r19_DEV.RData") ##### # REGION 19 (EVALreg19 -> r19_EVAL.RData) ##### region <- EVALreg19 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg19e save(FRPreg19e, file="r19_EVAL.RData") ##### # REGION 20 (DEVreg20 -> r20_DEV.RData) ##### region <- DEVreg20 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg20 save(FRPreg20, file="r20_DEV.RData") ##### # REGION 20 (EVALreg20 -> r20_EVAL.RData ##### region <- EVALreg20 #- Adding FRP0, lp0, cp0. region[,"FRP0"] <- c(0,region[1:(length(region[,"FRP"])-1),"FRP"]) region[,"lsp0"] <- c(0,region[1:(length(region[,"lsp"])-1),"lsp"]) region[,"cp0"] <- c(0,region[1:(length(region[,"cp"])-1),"cp"]) #- Select FRP Events FRPreg <- region[which(region$FRP>0),] #- Adding FRPratio FRPreg[,"ratio"] = FRPreg$FRP/FRPreg$FRP0 #- Converting temperature FRPreg[,"t"]<-FRPreg$t-273.15 #Kelvin to Celsius FRPreg[,"d"]<-FRPreg$d-273.15 #Kelvin to Celsius #- Adding Total Precipitation and tp0 FRPreg[,"tp"] = FRPreg$lsp+FRPreg$cp FRPreg[,"tp0"] = FRPreg$lsp0+FRPreg$cp0 #- Angstrom Index (and by necessity relative humidity) #-- Variables for relative humidity calculation a<-17.271 b<-237.7 #-- Relative Humidity calculation FRPreg[,"rh"]<-100*(exp((a*FRPreg$d)/(b+FRPreg$d))/exp((a*FRPreg$t)/(b+FRPreg$t))) #-- Angstrom Index calculation FRPreg[,"ai"]<-(FRPreg$rh/20)+((27-FRPreg$t)/10) FRPreg->FRPreg20e save(FRPreg20e, file="r20_EVAL.RData") ###### # Alternatively: I load the regions as .RData files for analysis and modeling ##### load("r01_DEV.RData") load("r02_DEV.RData") load("r04_DEV.RData") load("r12_DEV.RData") load("r13_DEV.RData") load("r14_DEV.RData") load("r18_DEV.RData") load("r19_DEV.RData") load("r20_DEV.RData") load("r01_EVAL.RData") load("r02_EVAL.RData") load("r04_EVAL.RData") load("r12_EVAL.RData") load("r13_EVAL.RData") load("r14_EVAL.RData") load("r18_EVAL.RData") load("r19_EVAL.RData") load("r20_EVAL.RData") region01name <- "Af" region02name <- "Am" region04name <- "Aw" region12name <- "Csa" region13name <- "Csb" region14name <- "Csc" region18name <- "Dfa" region19name <- "Dfb" region20name <- "Dfc" E <- "Equatorial" WT <- "Warm Temperate" S <- "Snow" ##### # Here I sort for FRPratio for both the DEVELOPMENT and EVALUATION data ##### region01<-subset(FRPreg1,FRPreg1$ratio!=Inf&FRPreg1$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,rh,wind,ai,region)) region02<-subset(FRPreg2,FRPreg2$ratio!=Inf&FRPreg2$ratio<2&FRPreg2$ratio>0,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,rh,wind,ai,region)) region04<-subset(FRPreg4,FRPreg4$ratio!=Inf&FRPreg4$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,rh,wind,ai,region)) region12<-subset(FRPreg12,FRPreg12$ratio!=Inf&FRPreg12$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region13<-subset(FRPreg13,FRPreg13$ratio!=Inf&FRPreg13$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region14<-subset(FRPreg14,FRPreg14$ratio!=Inf&FRPreg14$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region18<-subset(FRPreg18,FRPreg18$ratio!=Inf&FRPreg18$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region19<-subset(FRPreg19,FRPreg19$ratio!=Inf&FRPreg19$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region20<-subset(FRPreg20,FRPreg20$ratio!=Inf&FRPreg20$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region01e<-subset(FRPreg1e,FRPreg1e$ratio!=Inf&FRPreg1e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region02e<-subset(FRPreg2e,FRPreg2e$ratio!=Inf&FRPreg2e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region04e<-subset(FRPreg4e,FRPreg4e$ratio!=Inf&FRPreg4e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region12e<-subset(FRPreg12e,FRPreg12e$ratio!=Inf&FRPreg12e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region13e<-subset(FRPreg13e,FRPreg13e$ratio!=Inf&FRPreg13e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region14e<-subset(FRPreg14e,FRPreg14e$ratio!=Inf&FRPreg14e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region18e<-subset(FRPreg18e,FRPreg18e$ratio!=Inf&FRPreg18e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region19e<-subset(FRPreg19e,FRPreg19e$ratio!=Inf&FRPreg19e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) region20e<-subset(FRPreg20e,FRPreg20e$ratio!=Inf&FRPreg20e$ratio<2,select=c(FRP,FRP0,ratio,lon,lat,year,day,lsp,lsp0,cp,cp0,tp,tp0,t,d,swc,wind,rh,ai,region)) rm(FRPreg1, FRPreg2, FRPreg4, FRPreg12, FRPreg13, FRPreg14, FRPreg18, FRPreg19, FRPreg20) rm(FRPreg1e, FRPreg2e, FRPreg4e, FRPreg12e, FRPreg13e, FRPreg14e, FRPreg18e, FRPreg19e, FRPreg20e) ##### # Here I create the three zones ##### equatorial<-rbind(region01,region02,region04) warmtemperate<-rbind(region12,region13,region14) snow<-rbind(region18,region19,region20) ##### # Here I create the matrices needed for mapping the variables ##### comp01 <- region01[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp02 <- region02[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp04 <- region04[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp12 <- region12[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp13 <- region13[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp14 <- region14[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp18 <- region18[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp19 <- region19[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] comp20 <- region20[,c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")] world <- rbind(comp01,comp02,comp04,comp12,comp13,comp14,comp18,comp19,comp20) outfilepath <- "/Users/SeanGDemet/Data/Regions/global.txt" write.table(world, row.names=FALSE, col.names=FALSE, file = outfilepath) work <- read.table("global.txt",col.names=c("lon","lat","year","day","FRP","ratio","lsp","cp","tp","t","rh")) # Here I determine the amount of rain per fire event rain = ddply(region20,~lon+lat+year,function(x){sum(x$tp)}) median(rain$V1) # Here I create global datasets for each variable included in the analysis library(plyr) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$FRP),]}) global_frp <- max[,c(1:3,5)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_frp.txt" write.table(global_frp, row.names=FALSE, col.names=FALSE, file = outfilepath) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$ratio),]}) global_ratio <- max[,c(1:3,6)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_ratio.txt" write.table(global_ratio, row.names=FALSE, col.names=FALSE, file = outfilepath) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$lsp),]}) global_lsp <- max[,c(1:3,7)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_lsp.txt" write.table(global_lsp, row.names=FALSE, col.names=FALSE, file = outfilepath) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$cp),]}) global_cp <- max[,c(1:3,8)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_cp.txt" write.table(global_cp, row.names=FALSE, col.names=FALSE, file = outfilepath) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$tp),]}) global_tp <- max[,c(1:3,9)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_tp.txt" write.table(global_tp, row.names=FALSE, col.names=FALSE, file = outfilepath) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$t),]}) global_t <- max[,c(1:3,10)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_t.txt" write.table(global_t, row.names=FALSE, col.names=FALSE, file = outfilepath) max = ddply(work,~lon+lat+year,function(x){x[which.max(x$rh),]}) global_rh <- max[,c(1:3,11)] outfilepath <- "/Users/SeanGDemet/Data/Regions/global_rh.txt" write.table(global_rh, row.names=FALSE, col.names=FALSE, file = outfilepath) # Here I divide the global datasets by year read.table("global_frp.txt", col.names=c("lon","lat","year","FRP")) -> work frp2003 <- work[which(work$year==2003),] frp2004 <- work[which(work$year==2004),] frp2005 <- work[which(work$year==2005),] frp2006 <- work[which(work$year==2006),] frp2007 <- work[which(work$year==2007),] read.table("global_cp.txt", col.names=c("lon","lat","year","cp")) -> work cp2003 <- work[which(work$year==2003),] cp2004 <- work[which(work$year==2004),] cp2005 <- work[which(work$year==2005),] cp2006 <- work[which(work$year==2006),] cp2007 <- work[which(work$year==2007),] read.table("global_lsp.txt", col.names=c("lon","lat","year","lsp")) -> work lsp2003 <- work[which(work$year==2003),] lsp2004 <- work[which(work$year==2004),] lsp2005 <- work[which(work$year==2005),] lsp2006 <- work[which(work$year==2006),] lsp2007 <- work[which(work$year==2007),] read.table("global_ratio.txt", col.names=c("lon","lat","year","ratio")) -> work ratio2003 <- work[which(work$year==2003),] ratio2004 <- work[which(work$year==2004),] ratio2005 <- work[which(work$year==2005),] ratio2006 <- work[which(work$year==2006),] ratio2007 <- work[which(work$year==2007),] read.table("global_rh.txt", col.names=c("lon","lat","year","rh")) -> work rh2003 <- work[which(work$year==2003),] rh2004 <- work[which(work$year==2004),] rh2005 <- work[which(work$year==2005),] rh2006 <- work[which(work$year==2006),] rh2007 <- work[which(work$year==2007),] read.table("global_t.txt", col.names=c("lon","lat","year","t")) -> work t2003 <- work[which(work$year==2003),] t2004 <- work[which(work$year==2004),] t2005 <- work[which(work$year==2005),] t2006 <- work[which(work$year==2006),] t2007 <- work[which(work$year==2007),] read.table("global_tp.txt", col.names=c("lon","lat","year","tp")) -> work tp2003 <- work[which(work$year==2003),] tp2004 <- work[which(work$year==2004),] tp2005 <- work[which(work$year==2005),] tp2006 <- work[which(work$year==2006),] tp2007 <- work[which(work$year==2007),] # Here I write the Mapping Functions (One with log scale, one with normal scale) library(scales) library(ggplot2) library(maps) mapping <- function(data) { colnames(data) <- c("Longitude","Latitude","year","Variable") global <- map_data("world") p <- (ggplot(aes(x=Longitude,y=Latitude,fill=Variable),data=data) + geom_tile()) p <- p + geom_polygon(data=global,aes(x=long, y=lat, group=group), colour="black", fill="white", alpha=0) p <- p + scale_fill_gradientn(colours = rev(rainbow(10)), trans = "log") return(p) } mapping2 <- function(data) { colnames(data) <- c("Longitude","Latitude","year", "Variable") global <- map_data("world") p <- (ggplot(aes(x=Longitude,y=Latitude,fill=Variable),data=data) + geom_tile()) p <- p + geom_polygon(data=global,aes(x=long, y=lat, group=group), colour="black", fill="white", alpha=0) p <- p + scale_fill_gradientn(colours = c("blue","green","yellow","orange","red"), values = rescale(c(quantile(data$Variable))), guide = "colorbar", limits=c(round(min(data$Variable),digits=2),round(max(data$Variable),digits=2))) return(p) } # Here I generate the 35 Maps frpmap03 <- mapping2(frp2003) frpmap04 <- mapping2(frp2004) frpmap05 <- mapping2(frp2005) frpmap06 <- mapping2(frp2006) frpmap07 <- mapping2(frp2007) cpmap03 <- mapping2(cp2003) cpmap04 <- mapping2(cp2004) cpmap05 <- mapping2(cp2005) cpmap06 <- mapping2(cp2006) cpmap07 <- mapping2(cp2007) lspmap03 <- mapping2(lsp2003) lspmap04 <- mapping2(lsp2004) lspmap05 <- mapping2(lsp2005) lspmap06 <- mapping2(lsp2006) lspmap07 <- mapping2(lsp2007) ratiomap03 <- mapping2(ratio2003) ratiomap04 <- mapping2(ratio2004) ratiomap05 <- mapping2(ratio2005) ratiomap06 <- mapping2(ratio2006) ratiomap07 <- mapping2(ratio2007) rhmap03 <- mapping2(rh2003) rhmap04 <- mapping2(rh2004) rhmap05 <- mapping2(rh2005) rhmap06 <- mapping2(rh2006) rhmap07 <- mapping2(rh2007) tmap03 <- mapping2(t2003) tmap04 <- mapping2(t2004) tmap05 <- mapping2(t2005) tmap06 <- mapping2(t2006) tmap07 <- mapping2(t2007) tpmap03 <- mapping2(tp2003) tpmap04 <- mapping2(tp2004) tpmap05 <- mapping2(tp2005) tpmap06 <- mapping2(tp2006) tpmap07 <- mapping2(tp2007) # And here I add the titles for said maps frpmap03 <- frpmap03 + ggtitle("Maximum daily Fire Radiative Power (W/m2) for the nine selected regions in 2003") frpmap04 <- frpmap04 + ggtitle("Maximum daily Fire Radiative Power (W/m2) for the nine selected regions in 2004") frpmap05 <- frpmap05 + ggtitle("Maximum daily Fire Radiative Power (W/m2) for the nine selected regions in 2005") frpmap06 <- frpmap06 + ggtitle("Maximum daily Fire Radiative Power (W/m2) for the nine selected regions in 2006") frpmap07 <- frpmap07 + ggtitle("Maximum daily Fire Radiative Power (W/m2) for the nine selected regions in 2007") cpmap03 <- cpmap03 + ggtitle("Maximum Daily Convective Precipitation (m) for the nine selected regions in 2003") cpmap04 <- cpmap04 + ggtitle("Maximum Daily Convective Precipitation (m) for the nine selected regions in 2004") cpmap05 <- cpmap05 + ggtitle("Maximum Daily Convective Precipitation (m) for the nine selected regions in 2005") cpmap06 <- cpmap06 + ggtitle("Maximum Daily Convective Precipitation (m) for the nine selected regions in 2006") cpmap07 <- cpmap07 + ggtitle("Maximum Daily Convective Precipitation (m) for the nine selected regions in 2007") lspmap03 <- lspmap03 + ggtitle("Maximum Daily Large Scale Precipitation (m) for the nine selected regions in 2003") lspmap04 <- lspmap04 + ggtitle("Maximum Daily Large Scale Precipitation (m) for the nine selected regions in 2004") lspmap05 <- lspmap05 + ggtitle("Maximum Daily Large Scale Precipitation (m) for the nine selected regions in 2005") lspmap06 <- lspmap06 + ggtitle("Maximum Daily Large Scale Precipitation (m) for the nine selected regions in 2006") lspmap07 <- lspmap07 + ggtitle("Maximum Daily Large Scale Precipitation (m) for the nine selected regions in 2007") ratiomap03 <- ratiomap03 + ggtitle(expression(paste("Maximum Daily ",Delta,"FRP"," values for the nine selected regions in 2003",sep=""))) ratiomap04 <- ratiomap04 + ggtitle(expression(paste("Maximum Daily ",Delta,"FRP"," values for the nine selected regions in 2004",sep=""))) ratiomap05 <- ratiomap05 + ggtitle(expression(paste("Maximum Daily ",Delta,"FRP"," values for the nine selected regions in 2005",sep=""))) ratiomap06 <- ratiomap06 + ggtitle(expression(paste("Maximum Daily ",Delta,"FRP"," values for the nine selected regions in 2006",sep=""))) ratiomap07 <- ratiomap07 + ggtitle(expression(paste("Maximum Daily ",Delta,"FRP"," values for the nine selected regions in 2007",sep=""))) rhmap03 <- rhmap03 + ggtitle("Maximum Daily Relative Humidity (%) for the nine selected regions in 2003") rhmap04 <- rhmap04 + ggtitle("Maximum Daily Relative Humidity (%) for the nine selected regions in 2004") rhmap05 <- rhmap05 + ggtitle("Maximum Daily Relative Humidity (%) for the nine selected regions in 2005") rhmap06 <- rhmap06 + ggtitle("Maximum Daily Relative Humidity (%) for the nine selected regions in 2006") rhmap07 <- rhmap07 + ggtitle("Maximum Daily Relative Humidity (%) for the nine selected regions in 2007") tmap03 <- tmap03 + ggtitle("Maximum Daily Air Temperature (ºC) for the nine selected regions in 2003") tmap04 <- tmap04 + ggtitle("Maximum Daily Air Temperature (ºC) for the nine selected regions in 2004") tmap05 <- tmap05 + ggtitle("Maximum Daily Air Temperature (ºC) for the nine selected regions in 2005") tmap06 <- tmap06 + ggtitle("Maximum Daily Air Temperature (ºC) for the nine selected regions in 2006") tmap07 <- tmap07 + ggtitle("Maximum Daily Air Temperature (ºC) for the nine selected regions in 2007") tpmap03 <- tpmap03 + ggtitle("Maximum Daily Total Precipitation (m) for the nine selected regions in 2003") tpmap04 <- tpmap04 + ggtitle("Maximum Daily Total Precipitation (m) for the nine selected regions in 2004") tpmap05 <- tpmap05 + ggtitle("Maximum Daily Total Precipitation (m) for the nine selected regions in 2005") tpmap06 <- tpmap06 + ggtitle("Maximum Daily Total Precipitation (m) for the nine selected regions in 2006") tpmap07 <- tpmap07 + ggtitle("Maximum Daily Total Precipitation (m) for the nine selected regions in 2007") ##### # Here I generate the Correlation Tables # dvtableS | dvtableE | dvtableWT | ratiotable ##### dvtableE <- signif(cor(equatorial[,c("ratio","FRP","FRP0","lsp","lsp0","cp","cp0","tp","tp0","t","d", "swc", "wind","rh","ai")]), digits=2) dvtableWT <- signif(cor(warmtemperate[,c("ratio","FRP","FRP0","lsp","lsp0","cp","cp0","tp","tp0","t","d", "swc", "wind","rh","ai")]), digits=2) dvtableS <- signif(cor(snow[,c("ratio","FRP","FRP0","lsp","lsp0","cp","cp0","tp","tp0","t","d", "swc", "wind","rh","ai")]), digits=2) # This step isolates to the correlation coefficients for the FRPratio ONLY. ratiotable <- cbind(dvtableS[c(1,3:15),1],dvtableE[c(1,3:15),1],dvtableWT[c(1,3:15),1]) colnames(ratiotable) <- c("Equatorial","Warm Temperate","Snow") print(xtable(dvtableE[3:15,3:15]), floating.environment='sidewaystable') print(xtable(dvtableWT[3:15,3:15]), floating.environment='sidewaystable') print(xtable(dvtableS[3:15,3:15]), floating.environment='sidewaystable') print(xtable(ratiotable), floating=FALSE) ##### # Here I generate the Temporal AUTOCORRELATION of FRP Table for all zones. ##### r1 = acf(region01$FRP, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region02$FRP, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region04$FRP, lag.max=4, type="correlation", plot=FALSE) r4 = acf(equatorial$FRP, lag.max=4, type="correlation", plot=FALSE) acfFRP.E <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acfFRP.E) <- c(paste(region01name),paste(region02name),paste(region04name),"Equatorial") r1 = acf(region12$FRP, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region13$FRP, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region14$FRP, lag.max=4, type="correlation", plot=FALSE) r4 = acf(warmtemperate$FRP, lag.max=4, type="correlation", plot=FALSE) acfFRP.WT <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acfFRP.WT) <- c(paste(region12name),paste(region13name),paste(region14name),"Warm Temperate") r1 = acf(region18$FRP, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region19$FRP, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region20$FRP, lag.max=4, type="correlation", plot=FALSE) r4 = acf(snow$FRP, lag.max=4, type="correlation", plot=FALSE) acfFRP.S <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acfFRP.S) <- c(paste(region18name),paste(region19name),paste(region20name),"Snow") print(xtable(acfFRP.E), floating=FALSE) print(xtable(acfFRP.WT), floating=FALSE) print(xtable(acfFRP.S), floating=FALSE) # Here I determine the descriptive statistics of the ACF values for the whole data set. frp_acf1 <- c(acfFRP.E[2,],acfFRP.WT[2,],acfFRP.S[2,]) frp_acf2 <- c(acfFRP.E[3,],acfFRP.WT[3,],acfFRP.S[3,]) max(frp_acf1) min(frp_acf1) mean(frp_acf1) max(frp_acf2) min(frp_acf2) mean(frp_acf2) ##### # Here I generate the Temporal AUTOCORRELATION of LSP Table for all zones. ##### r1 = acf(region01$lsp, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region02$lsp, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region04$lsp, lag.max=4, type="correlation", plot=FALSE) r4 = acf(equatorial$lsp, lag.max=4, type="correlation", plot=FALSE) acflsp.E <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acflsp.E) <- c(paste(region01name),paste(region02name),paste(region04name),"Equatorial") r1 = acf(region12$lsp, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region13$lsp, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region14$lsp, lag.max=4, type="correlation", plot=FALSE) r4 = acf(warmtemperate$lsp, lag.max=4, type="correlation", plot=FALSE) acflsp.WT <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acflsp.WT) <- c(paste(region12name),paste(region13name),paste(region14name),"Warm Temperate") r1 = acf(region18$lsp, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region19$lsp, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region20$lsp, lag.max=4, type="correlation", plot=FALSE) r4 = acf(snow$lsp, lag.max=4, type="correlation", plot=FALSE) acflsp.S <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acflsp.S) <- c(paste(region18name),paste(region19name),paste(region20name),"Snow") print(xtable(acflsp.E), floating=FALSE) print(xtable(acflsp.WT), floating=FALSE) print(xtable(acflsp.S), floating=FALSE) # Here I determine the descriptive statistics of the ACF values for the whole data set. lsp_acf1 <- c(acflsp.E[2,],acflsp.WT[2,],acflsp.S[2,]) lsp_acf2 <- c(acflsp.E[3,],acflsp.WT[3,],acflsp.S[3,]) max(lsp_acf1) min(lsp_acf1) mean(lsp_acf1) max(lsp_acf2) min(lsp_acf2) mean(lsp_acf2) ##### # Here I generate the Temporal AUTOCORRELATION of CP Table for all zones. ##### r1 = acf(region01$cp, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region02$cp, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region04$cp, lag.max=4, type="correlation", plot=FALSE) r4 = acf(equatorial$cp, lag.max=4, type="correlation", plot=FALSE) acfcp.E <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acfcp.E) <- c(paste(region01name),paste(region02name),paste(region04name),"Equatorial") r1 = acf(region12$cp, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region13$cp, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region14$cp, lag.max=4, type="correlation", plot=FALSE) r4 = acf(warmtemperate$cp, lag.max=4, type="correlation", plot=FALSE) acfcp.WT <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acfcp.WT) <- c(paste(region12name),paste(region13name),paste(region14name),"Warm Temperate") r1 = acf(region18$cp, lag.max=4, type="correlation", plot=FALSE) r2 = acf(region19$cp, lag.max=4, type="correlation", plot=FALSE) r3 = acf(region20$cp, lag.max=4, type="correlation", plot=FALSE) r4 = acf(snow$cp, lag.max=4, type="correlation", plot=FALSE) acfcp.S <- cbind(r1$acf,r2$acf,r3$acf,r4$acf) colnames(acfcp.S) <- c(paste(region18name),paste(region19name),paste(region20name),"Snow") print(xtable(acfcp.E), floating=FALSE) print(xtable(acfcp.WT), floating=FALSE) print(xtable(acfcp.S), floating=FALSE) # Here I determine the descriptive statistics of the ACF values for the whole data set. cp_acf1 <- c(acfcp.E[2,],acfcp.WT[2,],acfcp.S[2,]) cp_acf2 <- c(acfcp.E[3,],acfcp.WT[3,],acfcp.S[3,]) max(cp_acf1) min(cp_acf1) mean(cp_acf1) max(cp_acf2) min(cp_acf2) mean(cp_acf2) ##### # Here I generate the cumulative HISTOGRAM of FRPratio for the three zones. ##### h <- hist(equatorial[,3], breaks=50, main=expression(paste("Histogram of ",Delta,"FRP"," for the Equatorial Zone",sep="")), xlab=expression(paste(Delta,"FRP",sep=""))) h$counts <- cumsum(h$counts) plot(h, main=expression(paste("Cumulative Histogram of ", Delta,"FRP"," for the Equatorial Zone",sep="")), xlab=expression(paste(Delta,"FRP",sep="")),sub=paste("Median:",format(median(equatorial$ratio),digits=3))) k <- hist(region20$ratio, breaks=50,main="Histogram of FRP/FRP0 for the E zone",xlab="FRP/FRP0 (W/m2)") k$counts <- cumsum(k$counts) plot(k, main=expression(paste("Cumulative Histogram of ", Delta,"FRP"," for the Dfc Region",sep="")), xlab=expression(paste(Delta,"FRP",sep="")),sub=paste("Median:",format(median(region20$ratio),digits=3))) i <- hist(warmtemperate$ratio, breaks=50,main="Histogram of FRP/FRP0 for the Warm Temperate Zone",xlab="FRP/FRP0 (W/m2)") i$counts <- cumsum(i$counts) plot(i, main=expression(paste("Cumulative Histogram of ", Delta,"FRP"," for the Warm Temperate Zone",sep="")), xlab=expression(paste(Delta,"FRP",sep="")),sub=paste("Median:",format(median(warmtemperate$ratio),digits=3))) j <- hist(snow$ratio, breaks=50,main="Histogram of FRP/FRP0 for the Snow Zone",xlab="FRP/FRP0 (W/m2)") j$counts <- cumsum(j$counts) plot(j, main=expression(paste("Cumulative Histogram of ", Delta,"FRP for the Snow Zone",sep="")), xlab=expression(paste(Delta,"FRP",sep="")),sub=paste("Median:",format(median(snow$ratio),digits=3))) #--- Simple Histogram preamble: Editing zone and zonename can provide histograms for any zone/region zone <- region20 zonename <- "Dfc" #--- par(mfrow=c(1,2)) teat <- hist(region20$day, breaks=1826, axes= FALSE, main="Occurrence of Fire Events over time for the Dfc region",xlab="", ylab="Frequency") axis(1, par(srt=45),at=c(0,182.5,365,548,731,913.5,1096,1278.5,1461,1643.5,1826),labels=c("Jan 2003","July 2003","Jan 2004","July 2004","Jan 2005","July 2005","Jan 2006","July 2006","Jan 2007","July 2007","Jan 2008"),las=2) axis(2) teatEz <- hist(equatorial$day, breaks=1826, axes= FALSE, main="Occurrence of Fire Events over time for the Equatorial zone",xlab="", ylab="Frequency") axis(1, par(srt=45),at=c(0,182.5,365,548,731,913.5,1096,1278.5,1461,1643.5,1826),labels=c("Jan 2003","July 2003","Jan 2004","July 2004","Jan 2005","July 2005","Jan 2006","July 2006","Jan 2007","July 2007","Jan 2008"),las=2) axis(2) teatWTz <- hist(warmtemperate$day, breaks=1826, axes= FALSE, main="Occurrence of Fire Events over time for the Warm Temperate zone",xlab="", ylab="Frequency") axis(1, par(srt=45),at=c(0,182.5,365,548,731,913.5,1096,1278.5,1461,1643.5,1826),labels=c("Jan 2003","July 2003","Jan 2004","July 2004","Jan 2005","July 2005","Jan 2006","July 2006","Jan 2007","July 2007","Jan 2008"),las=2) axis(2) teatSz <- hist(snow$day, breaks=1826, axes= FALSE, main="Occurrence of Fire Events over time for the Snow Zone",xlab="", ylab="Frequency") axis(1, par(srt=45),at=c(0,182.5,365,548,731,913.5,1096,1278.5,1461,1643.5,1826),labels=c("Jan 2003","July 2003","Jan 2004","July 2004","Jan 2005","July 2005","Jan 2006","July 2006","Jan 2007","July 2007","Jan 2008"),las=2) axis(2) hist(zone$t, breaks=50,main="Histogram of Air Temperature",sub=paste("Region:", zonename,"| Median:",format(median(zone$t),digits=3)),xlab="Air Temperature (°C)") hist(zone$tp,breaks=1000, xlim = range(0,0.01),main="Histogram of Total Precipitation",sub=paste("Region:", zonename,"| Max:", max(zone$tp),"| Median:",format(median(zone$tp),digits=3)),xlab="Total Precipitation (m)") hist(zone$rh,breaks=50,main="Histogram of Relative Humidity",sub=paste("Region:", zonename,"| Median:",format(median(zone$rh),digits=3)),xlab="Relative Humidity (%)") par(mfrow=c(3,1)) hist(zone$FRP, breaks=10000, xlim = range(0,0.01),main="Histogram of FRP",sub=paste("Region:", zonename,"| Max:", max(zone$FRP),"| Median:",format(median(zone$FRP),digits=3)),xlab=expression(paste("FRP (W/m"^2,")"))) hist(zone$d,breaks=50,main="Histogram of Dewpoint Temperature",sub=paste("Region:", zonename,"| Median:",format(median(zone$d),digits=3)),xlab="Dew Point Temperature (°C)") hist(zone$wind, breaks=50,main="Histogram of Wind",sub=paste("Region:", zonename,"| Median:",format(median(zone$wind),digits=3)),xlab="Wind (m/s)") hist(zone$swc, breaks=50,main="Histogram of Soil Water Content",sub=paste("Region:", zonename,"| Median:",format(median(zone$swc),digits=3)),xlab=expression(paste("Soil Water Content (m"^3,"/m"^3,")"))) hist(zone$lsp,breaks=1000,xlim = range(0,0.006),main="Histogram of Large Scale Precipitation",sub=paste("Region:", zonename,"| Max:", max(zone$lsp),"| Median:",format(median(zone$lsp),digits=3)),xlab="Large Scale Precipitation (m)") hist(zone$cp,breaks=1000, xlim = range(0,0.003),main="Histogram of Convective Precipitation",sub=paste("Region:", zonename,"| Max:", max(zone$cp),"| Median:",format(median(zone$cp),digits=3)),xlab="Convective Precipitation (m)") hist(zone$ai,breaks=50,main="Histogram of Angstrom Index",sub=paste("Region:", zonename,"| Median:",format(median(zone$ai),digits=3)),xlab="Angstrom Index") ##### # Here I generate the Rsquared Table ##### r2byvar <- function(var){ data <- c( format(summary(lm(region01[,3]~region01[,var]))$adj.r.squared, digits=3), format(summary(lm(region02[,3]~region02[,var]))$adj.r.squared, digits=3), format(summary(lm(region04[,3]~region04[,var]))$adj.r.squared, digits=3), format(summary(lm(equatorial[,3]~equatorial[,var]))$adj.r.squared, digits=3), format(summary(lm(region12[,3]~region12[,var]))$adj.r.squared, digits=3), format(summary(lm(region13[,3]~region13[,var]))$adj.r.squared, digits=3), format(summary(lm(region14[,3]~region14[,var]))$adj.r.squared, digits=3), format(summary(lm(warmtemperate[,3]~warmtemperate[,var]))$adj.r.squared, digits=3), format(summary(lm(region18[,3]~region18[,var]))$adj.r.squared, digits=3), format(summary(lm(region19[,3]~region19[,var]))$adj.r.squared, digits=3), format(summary(lm(region20[,3]~region20[,var]))$adj.r.squared, digits=3),format(summary(lm(snow[,3]~snow[,var]))$adj.r.squared, digits=3) ) return(data) } tabler2 <- cbind(r2byvar("FRP0"),r2byvar("lsp"),r2byvar("lsp0"),r2byvar("cp"),r2byvar("cp0"),r2byvar("tp"),r2byvar("tp0"),r2byvar("t"),r2byvar("d"),r2byvar("swc"),r2byvar("wind"),r2byvar("rh"),r2byvar("ai")) colnames(tabler2) <- c("FRP0","lsp","lsp0","cp","cp0","tp","tp0","t","d","swc","wind","rh","ai") rownames(tabler2) <- c(region01name, region02name, region04name, "Equatorial", region12name, region13name, region14name, "Warm Temperate", region18name, region19name, region20name, "Snow" print(xtable(tabler2), floating.environment='sidewaystable') ##### # Here I generate the SCATTERPLOTS # Preamble: Editing zone and zonename can provide histograms for any zone/region ##### zone <- region20 zonename <- "Dfc" #expression(paste("Maximum ",Delta,"FRP"," values for the nine selected regions in 2003",sep=""))) # 550 x 430 #-- par(mfrow=c(2,2)) #900 x 800 fit1 <-lm(zone$ratio~zone$FRP0) plot(zone$ratio~zone$FRP0,xlab=expression(paste("FRP"[0]," (W/m"^2,")")),ylab=expression(paste(Delta,"FRP")),main=expression(paste("FRP"[0]," vs. ",Delta,"FRP for ","Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$cp) plot(zone$ratio~zone$cp,xlab="Convective Precipitation (m)",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Convective Precipitation vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$t) plot(zone$ratio~zone$t,xlab="Air Temperature (ºC)",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Air Temperature vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$rh) plot(zone$ratio~zone$rh,xlab="Relative Humidity (%)",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Relative Humidity vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) par(mfrow=c(3,2)) #900 x 1100 fit1 <-lm(zone$ratio~zone$swc) plot(zone$ratio~zone$swc,xlab=expression(paste("Soil Water Content (m"^3,"/m"^3,")")),ylab=expression(paste(Delta,"FRP")),main=expression(paste("Soil Water Content vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$wind) plot(zone$ratio~zone$wind,xlab="Wind (m/s)",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Wind vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$lsp) plot(zone$ratio~zone$lsp,xlab="Large Scale Precipitation (m)",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Large Scale Precipitation vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$d) plot(zone$ratio~zone$d,xlab="Dew Point Temperature (ºC)",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Dew Point Temperature vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) fit1 <-lm(zone$ratio~zone$ai) plot(zone$ratio~zone$ai,xlab="Angstrom Index",ylab=expression(paste(Delta,"FRP")),main=expression(paste("Angstrom Index vs. ",Delta,"FRP for Dfc Region", sep="")), sub=paste("R2 is",format(summary(fit1)$adj.r.squared, digits=2))) #abline(fit1, col='red',lwd=2) ##### # Here I generate the BoxPlots for the data. # Preamble: Editing zone and zonename can provide histograms for any zone/region ##### zone <- region02 zonename <- "Am" #-- par(mfrow=c(2,re2)) tpfac <- cut(zone$FRP0,c(quantile(zone$FRP0))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5, main=expression(paste("Box Plot of FRP"[0]," with possible outliers")), ylab=expression(paste(Delta,"FRP")), xlab=expression(paste("FRP"[0]," (W/m"^2,")")),sub=paste("Region:",zonename,"| min:",format(min(zone$FRP0),digits=3),"| M:",format(median(zone$FRP0),digits=3),"| max:", format(max(zone$FRP0),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$cp,c(quantile(zone$cp))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main=expression("Box Plot of Convective Precipitation with possible outliers"), ylab=expression(paste(Delta,"FRP")),xlab="Convective Precipitation (m)",sub=paste("Region:",zonename,"| min:",format(min(zone$cp),digits=3),"| M:",format(median(zone$cp),digits=3),"| max:", format(max(zone$cp),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$t,c(quantile(zone$t))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main=expression("Box Plot of Air Temperature with possible outliers"), ylab=expression(paste(Delta,"FRP")),xlab="Air Temperature (ºC)",sub=paste("Region:",zonename,"| min:",format(min(zone$t),digits=3),"| M:",format(median(zone$t),digits=3),"| max:", format(max(zone$t),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$rh,c(quantile(zone$rh))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main=expression("Box Plot of Relative Humidity with possible outliers"), ylab=expression(paste(Delta,"FRP")),xlab="Relative Humidity (%)",sub=paste("Region:",zonename,"| min:",format(min(zone$rh),digits=3),"| M:",format(median(zone$rh),digits=3),"| max:", format(max(zone$rh),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) par(mfrow=c(3,2)) tpfac <- cut(zone$swc,c(quantile(zone$swc))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main="Box Plot of Soil Water Content with possible outliers", ylab=expression(paste(Delta,"FRP")),xlab=expression(paste("Soil Water Content (m"^3,"/m"^3,")")),sub=paste("Region:",zonename,"| min:",format(min(zone$swc),digits=3),"| M:",format(median(zone$swc),digits=3),"| max:", format(max(zone$swc),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$wind,c(quantile(zone$wind))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main="Box Plot of Wind Gusts with possible outliers", ylab=expression(paste(Delta,"FRP")),xlab="Wind Gusts (m/s)",sub=paste("Region:",zonename,"| min:",format(min(zone$wind),digits=3),"| M:",format(median(zone$wind),digits=3),"| max:", format(max(zone$wind),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$lsp,c(quantile(zone$lsp))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main="Box Plot of Large Scale Precipitation with poss. outliers", ylab=expression(paste(Delta,"FRP")),xlab="Large Scale Precipitation (m)",sub=paste("Region:",zonename,"| min:",format(min(zone$lsp),digits=3),"| M:",format(median(zone$lsp),digits=3),"| max:", format(max(zone$lsp),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$lsp0,c(quantile(zone$lsp0))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main="Box Plot of Large Scale Precipitation 0 with possible outliers", ylab=expression(paste(Delta,"FRP")),xlab="Large Scale Precipitation (m)",sub=paste("Region:",zonename,"| min:",format(min(zone$lsp0),digits=3),"| M:",format(median(zone$lsp0),digits=3),"| max:", format(max(zone$lsp0),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$d,c(quantile(zone$d))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main="Box Plot of Dew Point Temperature with possible outliers", ylab=expression(paste(Delta,"FRP")),xlab="Dew Point Temperature (ºC)",sub=paste("Region:",zonename,"| min:",format(min(zone$d),digits=3),"| M:",format(median(zone$d),digits=3),"| max:", format(max(zone$d),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) tpfac <- cut(zone$ai,c(quantile(zone$ai))) tpbox <- boxplot(zone$ratio~tpfac, range=1.5,main="Box Plot of Angstrom Index with possible outliers", ylab=expression(paste(Delta,"FRP")),xlab="Angstrom Index",sub=paste("Region:",zonename,"| min:",format(min(zone$ai),digits=3),"| M:",format(median(zone$ai),digits=3),"| max:", format(max(zone$ai),digits=3)), names = c("0-25%","25-50%","50-75%","75-100%")) lines(tpbox$stats[c(1),]) lines(tpbox$stats[c(2),]) lines(tpbox$stats[c(3),]) lines(tpbox$stats[c(4),]) lines(tpbox$stats[c(5),]) ##### # Functions needed for U-Test, in particular the tabular analysis. ##### # Function to determine Trend direction myfu <-function(data){ if(data > 0){ data = "+" }else if(data<0){ data = "-" }else if(data == 0){ data = "=" } } # Utest function for CP: function(region,"var") prows <- function(region,var,name){ utest.rI <- region[order(region[,paste(var)]),] q=round(length(region[,paste(var)])*0.05) r=round(length(region[,paste(var)])*0.25) s=round(length(region[,paste(var)])*0.50) t=round(length(region[,paste(var)])*0.75) u=round(length(region[,paste(var)])*0.95) utest.rI[c(1:q),"utest"] <- 1 utest.rI[c(q:r),"utest"] <- 2 utest.rI[c(r:s),"utest"] <- 3 utest.rI[c(s:t),"utest"] <- 4 utest.rI[c(t:u),"utest"] <- 5 utest.rI[c(u:length(region[,paste(var)])),"utest"] <- 6 i <- utest.rI[which(utest.rI[,"utest"] == 1),] ii <- utest.rI[which(utest.rI[,"utest"] == 2),] iii <- utest.rI[which(utest.rI[,"utest"] == 3),] iv <- utest.rI[which(utest.rI[,"utest"] == 4),] v <- utest.rI[which(utest.rI[,"utest"] == 5),] vi <- utest.rI[which(utest.rI[,"utest"] == 6),] ii.vi <- rbind(ii,iii,iv,v,vi) i.v <- rbind(i,ii,iii,iv,v) a = median(i[,"ratio"], na.rm=TRUE) b = median(ii[,"ratio"], na.rm=TRUE) c = median(iii[,"ratio"], na.rm=TRUE) d = median(iv[,"ratio"], na.rm=TRUE) e = median(v[,"ratio"], na.rm=TRUE) f = median(vi[,"ratio"], na.rm=TRUE) g = median(ii.vi[,"ratio"], na.rm=TRUE) h = median(i.v[,"ratio"], na.rm=TRUE) obj11 = g-a obj21 = b-a obj31 = c-b obj41 = d-c obj51 = e-d obj61 = f-e obj71 = f-h obj1 <- wilcox.test(i[,"ratio"],ii.vi[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj2 <- wilcox.test(i[,"ratio"],ii[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj3 <- wilcox.test(ii[,"ratio"],iii[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj4 <- wilcox.test(iii[,"ratio"],iv[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj5 <- wilcox.test(iv[,"ratio"],v[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj6 <- wilcox.test(v[,"ratio"],vi[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj7 <- wilcox.test(vi[,"ratio"],i.v[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) mat <- cbind(name,var, obj1[3], myfu(obj11), obj2[3], myfu(obj21), obj3[3], myfu(obj31), obj4[3], myfu(obj41), obj5[3],myfu(obj51), obj6[3], myfu(obj61), obj7[3], myfu(obj71)) return(mat) } # Utest function for ALL other variables: function(region,var) utest <- function(region,var,name){ utest.rI <- region[order(region[,paste(var)]),] utest.rI[,"utest"] <-cut(utest.rI[,paste(var)],c(quantile(utest.rI[,paste(var)], probs=c(0,0.05,0.25,0.50,0.75,0.95,1))), labels=FALSE) i <- utest.rI[which(utest.rI[,"utest"] == 1),] ii <- utest.rI[which(utest.rI[,"utest"] == 2),] iii <- utest.rI[which(utest.rI[,"utest"] == 3),] iv <- utest.rI[which(utest.rI[,"utest"] == 4),] v <- utest.rI[which(utest.rI[,"utest"] == 5),] vi <- utest.rI[which(utest.rI[,"utest"] == 6),] ii.vi <- rbind(ii,iii,iv,v,vi) i.v <- rbind(i,ii,iii,iv,v) a = median(i[,"ratio"], na.rm=TRUE) b = median(ii[,"ratio"], na.rm=TRUE) c = median(iii[,"ratio"], na.rm=TRUE) d = median(iv[,"ratio"], na.rm=TRUE) e = median(v[,"ratio"], na.rm=TRUE) f = median(vi[,"ratio"], na.rm=TRUE) g = median(ii.vi[,"ratio"], na.rm=TRUE) h = median(i.v[,"ratio"], na.rm=TRUE) obj11 = g-a obj21 = b-a obj31 = c-b obj41 = d-c obj51 = e-d obj61 = f-e obj71 = f-h obj1 <- wilcox.test(i[,"ratio"],ii.vi[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj2 <- wilcox.test(i[,"ratio"],ii[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj3 <- wilcox.test(ii[,"ratio"],iii[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj4 <- wilcox.test(iii[,"ratio"],iv[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj5 <- wilcox.test(iv[,"ratio"],v[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj6 <- wilcox.test(v[,"ratio"],vi[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) obj7 <- wilcox.test(vi[,"ratio"],i.v[,"ratio"], alternative="t",exact=FALSE, correct=TRUE) mat <- cbind(name,var, obj1[3], myfu(obj11), obj2[3], myfu(obj21), obj3[3], myfu(obj31), obj4[3], myfu(obj41), obj5[3],myfu(obj51), obj6[3], myfu(obj61), obj7[3], myfu(obj71)) return(mat) } # Star Function. For table presentation. addstars <- function(data){ lim <- (0.05/7) for (u in c(1:length(data[,1]))){ if (data[u,3]