CCISS.Predict_Mar21.R

#Hi Ben or whomever ends up working on this code! Assuming you already have the random forest model built,
#this code will take the new data, predict the subzones, use the edatopic overlap to determine likley site
#series, and produce summary tables. The first part of the code up to the function and loop section as line 175
#hasn't really been changed, but almost everything after that has been. For creating the summary statistics,
#there are some cutoffs/constants that might be changed or could even be adjustable on the web tool. I've 
#tried to comment these. Please don't hesitate to get in touch if anything doesn't makes sense. Thanks!

.libPaths("E:/R packages")
require (RGtk2)
require(plyr)
require (rChoiceDialogs)
require (data.table)
require(doBy)
require (utils)
require(labdsv)
require(tools )
require(svDialogs)
require(tcltk)
require(randomForest)
require(foreach)
require(dplyr)
require(reshape2)
require(reshape)
library(doParallel)
require(data.table)
#===============================================================================
# Clear environment and add in previous rF model 
#===============================================================================

rm(list=ls())

wd=tk_choose.dir()
setwd(wd)

###enter model file name
fname="BGCv10_2000Pt_Rnd_Normal_1961_1990MSY_RFmodelKiriFinal"
load(paste(fname,".Rdata",sep=""))

#choose data file
fplot=(file.choose())

Columns = c("GCM", "ID1", "ID2", "Latitude", "Longitude", "Elevation", "AHM", "bFFP",
            "CMD07","DD5_sp","EMT","Eref_sm","EXT","FFP","MCMT","MSP",
            "PPT07","PPT08", "PPT05","PPT06","PPT09", "SHM","TD","Tmax_sp","Tmin_at",
            "Tmin_sm","Tmin_wt", "PPT_at","PPT_wt", "PAS","eFFP",
            "Eref09","MAT","Tmin_sp","CMD")

Y1 <- fread(fplot, select = Columns, stringsAsFactors = FALSE, data.table = FALSE) #fread is faster than read.csv
fplot=basename(file_path_sans_ext(fplot))

gc()


colnames (Y1) [1:3] = c("Model", "SiteNo", "BGC")

Y1$Model <- as.factor (Y1$Model)
Y1$SiteNo <- as.factor(Y1$SiteNo)
Y1$BGC <-gsub(" ", "", Y1$BGC, fixed = TRUE)
Y1$BGC <- as.factor(Y1$BGC)

#####generate some addtional variables
####
Y1$PPT_MJ <- Y1$PPT05 + Y1$PPT06 # MaY/June precip
Y1$PPT_JAS <- Y1$PPT07 + Y1$PPT08 + Y1$PPT09 # July/Aug/Sept precip
Y1$PPT.dormant <- Y1$PPT_at + Y1$PPT_wt # for calculating spring deficit
Y1$CMD.def <- 500 - (Y1$PPT.dormant)# start of growing season deficit original value was 400 but 500 seems better
Y1$CMD.def [Y1$CMD.def < 0] <- 0 #negative values set to zero = no deficit
Y1$CMDMax <- Y1$CMD07
Y1$CMD.total <- Y1$CMD.def + Y1$CMD

##Specify variables in model
model = "5.1"
VarList = c("AHM", "bFFP","CMD.total","DD5_sp","EMT","Eref_sm","EXT","FFP","MCMT","MSP",
            "PPT_JAS","PPT_MJ","PPT06","SHM","TD","Tmax_sp","Tmin_at","Tmin_sm","Tmin_wt",
            "PAS","CMD.def","CMDMax","eFFP","Eref09","MAT","PPT07","Tmin_sp")
List = c("Model", "SiteNo", "BGC")
Y1save = Y1
Y1$BGC  <- as.factor(Y1$BGC)
Y1.sub=Y1[,names(Y1) %in% c(List,VarList)]

Y1.sub$BGC <- gsub(" ", "", Y1.sub$BGC, fixed = TRUE)

Y1.sub$BGC <- factor(Y1.sub$BGC, levels=unique(Y1.sub$BGC))

##Predict future subzones######
Y1.sub$BGC.pred <- predict(BGCmodel, Y1.sub[,-c(1:2)]) 
gc()

Y1.sub$PlotNo <- attr(Y1.sub, "row.names")
S1 <- subset(Y1, select=c("Model", "SiteNo", "Latitude", "Longitude", "Elevation"))
S2 <- subset(Y1.sub, select=c("BGC", "BGC.pred"))   

Y2.sub <- cbind.data.frame(S1, S2)

Y2.sub <- Y2.sub[,c("Model", "SiteNo", "BGC", "BGC.pred", "Latitude", "Longitude", "Elevation")]

Y2.sub$BGC.pred <- gsub(" ", "", Y2.sub$BGC.pred, fixed = TRUE)

Y2.sub$BGC.pred <- factor(Y2.sub$BGC.pred, levels=unique(Y2.sub$BGC.pred))

Y2.sub$BGC <- gsub(" ", "", Y2.sub$BGC, fixed = TRUE)

Y2x <- as.character(Y2.sub$Model)
Ystr <- strsplit(Y2x, "_")
Y4 <- matrix(unlist(Ystr), ncol=3, byrow=TRUE)
Y2.sub <- cbind(Y4, Y2.sub)
colnames(Y2.sub)[1:3]=c("GCM","Scenario", "FuturePeriod" )

Y3.sub = Y2.sub

Y3.sub$Num <- 1

Y3.sub$BGC.pred <- gsub(" ", "", Y3.sub$BGC.pred, fixed = TRUE)


#===============================================================================
#Build Summarys of BGC prediction ratios
#===============================================================================

m <- summaryBy(SiteNo ~ BGC + FuturePeriod, data=Y3.sub, FUN=c(length))

p <- summaryBy(SiteNo ~ BGC + BGC.pred + FuturePeriod, data=Y3.sub,  FUN=c(length))

BGCratio <- merge(m,p,by=c("BGC","FuturePeriod"))

BGCratio$FutureRatio <- BGCratio$SiteNo.length.y/BGCratio$SiteNo.length.x

BGCratio$BGC.pred <- gsub(" ", "", BGCratio$BGC.pred, fixed = TRUE)

BGCratio2 <- BGCratio[, c("BGC", "BGC.pred", "FutureRatio")]

BGCratio2$BGC <- as.character(BGCratio2$BGC)
BGCratio2$BGC <- gsub("[[:space:]]","",BGCratio2$BGC) ##just incase any spaces got added
BGCratio2$BGC.pred <- gsub("[[:space:]]","",BGCratio2$BGC.pred)
#===============================================================================
#
#          Import suitability and edatopic tables
#           
#===============================================================================

edatopename="EdatopicWithSpecialCurrent"
edatopename2=paste(wd,"/",edatopename,".csv",sep="")
E1 <-read.csv(edatopename2,stringsAsFactors=FALSE,na.strings=".")

E1 <- E1[,-c(5:6)]
E1 <- unique(E1)

###create list of focal BGCs & edatopic space
e1 <- as.list(unique(BGCratio2$BGC), all.names=FALSE)
edatopic1 <- E1[E1$MergedBGC %in% e1,]
edatopic1$Codes[edatopic1$Codes == ""] <- NA

###create list of predicted BGCs and edatopic space
e2 <- as.list (unique(BGCratio2$BGC.pred), all.names=FALSE)
edatopic2 <- E1[E1$MergedBGC %in% e2,]
edatopic2$Codes[edatopic2$Codes == ""] <- NA

#################Import and Build Tree species suitability##############
treesuit="TreeSppSuit_v10.5"
treesuit2=paste(wd,"/",treesuit,".csv",sep="")
S1 <- read.csv(treesuit2,stringsAsFactors=F,na.strings=".")
S1 <- unique(S1)

#===============================================================================
# Builds list of all BGCs, Future BGCs, and Site Series
#===============================================================================

BGClist = (unique(BGCratio2$BGC))
FuturePeriod.list <- as.list(unique(Y2.sub$FuturePeriod))
BGCfutures.list <- as.list(unique(Y3.sub$BGC.pred)) ### to use later on to limit the site series
BGCfocalE <- edatopic1[edatopic1$MergedBGC %in% Y3.sub$BGC  ,] ### extracts edatopic space for BGC focal of SiteNo
BGCfutureE <- edatopic2[edatopic2$MergedBGC %in% Y3.sub$BGC.pred  ,] #extracts edatopic info only for predicted BGCs
Y3.sub$SiteNo <- as.character(Y3.sub$SiteNo)
SiteNo.list = as.list(unique(Y3.sub$SiteNo))
Y3.sub1 <-Y3.sub
Y3.sub1$BGC <- gsub("[[:space:]]","",Y3.sub1$BGC)
Y3.sub1$BGC.pred <- gsub("[[:space:]]","",Y3.sub1$BGC.pred)
gc()

###OPTIONAL: Set up to run loops in parallel###
#require(doParallel)
#cl = makePSOCKcluster(6)
#registerDoParallel(cl)
#clusterEvalQ(cl, .libPaths("E:/R packages"))

#=======================================================================
#             LOAD FUNCTIONS USED IN LOOP AND FOR SUMMARY STATS-- NEW!!!
#=======================================================================

###Function to find the proportion edatopic overlap###
edaOverlap <- function(dat1, dat2){
  return(length(dat1[dat1 %in% dat2])/length(dat2))
}

######Functions used in suitability rules####
bifurc <- function(x){ ###Check for bifurcated model in raw data
  t1 <- x[1] + x[2]
  cutoff1 <- 0.85 - t1
  t2 <- x[3] + x[4]
  cutoff2 <- 0.85 - t2
  if(max(x) > 0.65 | cutoff1 < 0.2 | cutoff2 < 0.2){
    return(FALSE)
  }else if(x[4] >= cutoff1 | x[1] >= cutoff2){
    return(TRUE)
  }else{
    return(FALSE)
  }
}

bifurcTrend <- function(x){ ##percent improve, stable, and same - bifurcated model in summary data
  if(x[1] >= 25 & x[3] >= 25){
    return("*")
  }else{
    return("")
  }
}

newSuit <- function(x){ ##Calculate new suitability accounting for current suitability
  suitLookup <- data.frame(SuitCurr = c(1,2,3,5), Value = c(-0.5, -0.3, 0.3, 0.6)) ##VALUES FOR CURRENT SUIT COULD BE ADJUSTED
  val <- suitLookup$Value[match(x[5], suitLookup$SuitCurr)]
  return(val+(1*x[1]+2*x[2]+3*x[3]+5*x[4]))
}

newSuitnoCurrent <- function(x){ ##New suitability without current suitability
  return(1*x[1]+2*x[2]+3*x[3]+5*x[4])
}

###stepSum calculates the difference between suitability in each time period####
stepSum <- function(x){##x is dataframe with colums Period, Spp, CurrentSuit, NewSuit
  if(x[1] == 2025){
    return(as.numeric(x[3])-as.numeric(x[4]))
  }else if(x[1] == 2055){
    y <- as.numeric(rawDat[rawDat$Spp == x[2] & rawDat$FuturePeriod == 2025, "NewSuit"])
    return(y - as.numeric(x[4]))
  }else{
    y <- as.numeric(rawDat[rawDat$Spp == x[2] & rawDat$FuturePeriod == 2055, "NewSuit"])
    return(y - as.numeric(x[4]))
  }
}

###Calculate percent of models improving, same, or declining####
modDir2 <- function(x, direction){#columns 1,2,3,X,Currentsuit
  curr <- x[5]
  if(curr == 5 | is.na(curr)){
    curr <- 4
  }
  x[5] <- 0
  x <- c(0,x)
  improve <- sum(x[1:curr])
  same <- x[curr + 1]
  decline <- sum(x[(curr+2):length(x)])
  if(direction == "Improve"){
    return(round(improve*100, digits = 0))
  }else if(direction == "Stable"){
    return(round(same*100, digits = 0))
  }else if(direction == "Decline"){
    return(round(decline*100, digits = 0))
  }else{
    return("error")
  }
}
#=======================================================================================  
#####Nested foreach loops to caclulate site series - now matches special SS together#####
#    To run loops in parallel, change the %do% in the first loop to %dopar%
#    Would also need to load all packages used inside loop onto workers using .packages command
#========================================================================================

SiteNo.suit <-  foreach(SNL = SiteNo.list, .combine = rbind) %do% {##  for each SiteNo in the data
  
  options(warn=2)
  
  cat("===========================================","\n")
  cat(SNL, "\n", "===========================================","\n")
  
  SiteFuture.suit <- foreach(i = FuturePeriod.list, .combine = rbind)  %do% {
    
    Y3.each <- Y3.sub1[Y3.sub1$SiteNo %in% SNL ,] ## extracts data for each site
    Y3.each <- Y3.each[Y3.each$FuturePeriod %in% i,] ##extracts data for each time period
    Y3.siteno <-as.list(unique (Y3.each$SiteNo))
    Y3.BGC <- as.list(unique(Y3.each$BGC))
    Y3.BGC.pred<- unique(Y3.each$BGC.pred)
    BGCfocalE <- edatopic1[edatopic1$MergedBGC %in% Y3.BGC  , ] ### extracts edatopic space for BGC focal of SiteNo
    BGCfutureE <- edatopic2[edatopic2$MergedBGC %in% Y3.BGC.pred  , ] #extracts edatopic info only for predicted BGCs
    ##Y3.SSlist = "SBSmc2/01"
    Y3.SSlist = as.list(unique(BGCfocalE$SS_NoSpace))
    
    FTS2 <-  foreach(SS = Y3.SSlist, .combine =rbind) %do% {      ##  for each site series for a SiteNo BGC
  
      options(warn=2)
      
      cat("===========================================","\n")
      
      cat(SS, "\n", "===========================================","\n")
      
      SSfocal <- BGCfocalE[BGCfocalE$SS_NoSpace %in% SS ,] ###find focal site series cells
      SSfocallist <- as.list(unique(SSfocal$SS_NoSpace))
      SSfocalBGClist <- as.list(unique(SSfocal$MergedBGC))
      SSfocalE <- as.list(unique(SSfocal$Edatopic))
      
      ##select site series only with some edatopic overlap with SSfocal
      
      SSfutureE <- BGCfutureE[BGCfutureE$Edatopic %in% SSfocalE,]
      futureZones <- unique(SSfutureE$MergedBGC)
      futureSS.names <- unique(SSfutureE$SS_NoSpace)
      if(length(SSfutureE$Edatopic) > 0){
        
        SSfocal <- BGCfocalE[BGCfocalE$SS_NoSpace %in% SS,]
        SSfuture <- BGCfutureE[BGCfutureE$SS_NoSpace %in% futureSS.names,]
        
        ##match site series within each projected subzone
        futureSS <- foreach(futSS = futureZones, .combine = rbind) %do% {
          dat <- SSfuture[SSfuture$MergedBGC == futSS,]
          if(any(!is.na(SSfocal$Codes)) & any(dat$Codes %in% SSfocal$Codes)){###Are there matchin special edatopic cells?
            fut <- dat[dat$Edatopic %in% SSfocal$Edatopic,]
            oldSp <- unique(SSfocal$Codes[!is.na(SSfocal$Codes)])
            newSp <- unique(fut$SS_NoSpace[match(oldSp, fut$Codes)])
            dat <- unique(dat[dat$SS_NoSpace == newSp, -c(4:5)]) #which ones have the new special edatope
            dat$overlap <- 1
            dat$revoverlap <- 1
          }else{
            dat <- dat[is.na(dat$Codes),]
            ##if not special, loop through each SS to calculate overlap
            dat <- foreach(x = unique(as.character(dat$SS_NoSpace)), .combine = rbind) %do% {
              dat1 <- dat[dat$SS_NoSpace == x,]
              Overlap <- edaOverlap(SSfocal$Edatopic, dat1$Edatopic)
              Revoverlap <- edaOverlap(dat1$Edatopic, SSfocal$Edatopic)
              dat1 <- unique(dat1[-c(4:5)])
              dat1$overlap <- Overlap
              dat1$revoverlap <- Revoverlap
              dat1 <- as.data.frame(dat1)
            }
          }
          dat
        }
        ##Multiply forwards and backwards overlap
        futureSS$alloverlap <- futureSS$overlap*futureSS$revoverlap
        
        ####add BGC probability to FutureSS list (this is basically the same as Will had)
        m2 <- summaryBy(SiteNo~BGC, data=Y3.each, FUN=c(length))
        p2 <- summaryBy(SiteNo~BGC + BGC.pred, data=Y3.each,  FUN=c(length))
        BGCratio3 <- merge (m2,p2, by.x = "BGC", by.y = "BGC" )
        BGCratio3$FutureRatio <- BGCratio3$SiteNo.length.y/BGCratio3$SiteNo.length.x
        summaryBy(FutureRatio ~ BGC, data=BGCratio3, FUN=c(sum)) #### summary to check that equals 1 for each BGC
        
        futureSS$BGCratio <- BGCratio3$FutureRatio[match(futureSS$MergedBGC, BGCratio3$BGC.pred)]
        
        ####Calculate the SS ratio
        SSoverlap <- summaryBy(alloverlap~MergedBGC, data=futureSS, id = 'SS_NoSpace', FUN=c(sum))
        futureSS$overlaptot<- SSoverlap$alloverlap.sum[match(futureSS$MergedBGC, SSoverlap$MergedBGC )]
        futureSS$SSratio <- futureSS$alloverlap/futureSS$overlaptot
        summaryBy(SSratio ~ MergedBGC, data=futureSS, FUN=c(sum)) #### for checking that SSratio sums to 100%
        summaryBy(SSratio ~ SS_NoSpace, data=futureSS, FUN=c(sum)) #
        
        ####Calculated the overall site series probability
        futureSS$SSprob <- (futureSS$BGCratio * futureSS$SSratio)
        futureSS$SSCurrent <- rep(SS,length(futureSS$SSprob))
        
        futureSS$FuturePeriod <- as.character(i)  
        futureSS$SiteNo <- as.character(SNL)
        
        futureSS <- as.data.frame(futureSS)
      }
      
    } #For each Site
  } #For each year
  
} # for all

#################End of FOREACH LOOP###################
#######################################################
#####SiteNo.suit now contains projected BGC units as SS_NoSpace######

####Merge with Spp suitability and create output tables#####
############################################################

comb <- function(...) {##Combine multiple dataframe in foreach loop
  mapply(FUN = rbind, ..., SIMPLIFY=FALSE)
}
#===================================================================================
#####Foreach loops to calculate summary statistics within each current BGC unit#####
######allOutput is a list with 2 dataframes#########################################
#===================================================================================
allOutput <- foreach(Site = unique(SiteNo.suit$SiteNo), .combine =  comb, .multicombine = TRUE) %:%
  foreach(SS = unique(SiteNo.suit$SSCurrent), .combine =  comb, .multicombine = TRUE) %do%{
    test <- SiteNo.suit[SiteNo.suit$SSCurrent == SS &SiteNo.suit$SiteNo == Site ,c(13,12,11,3,10)]
    SuitSx <- S1[,c(2:4)]
    SuitSx <- SuitSx[!duplicated(SuitSx[,1:2]),]
    comb <- merge(test, SuitSx, by.x = "SS_NoSpace", by.y = "Unit", all.x = TRUE) ##Merge SS with suitability table
    comb <- comb[complete.cases(comb),]
    comb <- comb[order(comb$SS_NoSpace, comb$FuturePeriod),]
    comb$Int <- interaction(comb$SS_NoSpace,comb$FuturePeriod)
    index <- which(comb$Int[-1] != comb$Int[-length(comb$Int)]) #index at each factor level
    index <- c(0,index)
    allSp <- unique(comb$Spp)
    allSp <- allSp[!is.na(allSp)]
    
    ####following loop creates rows for not suitable species in each projected Site Series###
    for(i in 1:(length(index) - 1)){
      Sp <- comb$Spp[(index[i]+1):index[i+1]]
      missing <- allSp[!allSp %in% Sp]
      new <- comb[rep(index[i+1], length(missing)),]
      new$Spp <- missing
      new$Suitability <- NA
      comb <- rbind(comb, new)
    }
    
    comb <- comb[order(comb$SS_NoSpace, comb$FuturePeriod),]
    #comb <- comb[!is.na(comb$Spp),]
    
    comb$Suitability[is.na(comb$Suitability)] <- "X"
    #comb <- comb[complete.cases(comb),]
    
    ####new data frame with proportion votes for each suit class##
    numVotes <- cast(comb, Spp + FuturePeriod + SSCurrent ~ Suitability, value = "SSprob", fun.aggregate = sum) ###votes for each suitability
    numVotes$Sum <- rowSums(numVotes[,c(4:7)]) ##check votes sum to 1 (ignoring rounding errors)
    
    ###merge in current suitability
    SuitSub <- SuitSx[SuitSx$Unit == unique(numVotes$SSCurrent),2:3]
    numVotes <- merge(numVotes, SuitSub, by = "Spp", all.x = TRUE)
    numVotes$Suitability[is.na(numVotes$Suitability)] <- 5
  
    ###Create statistics for raw data####
    rawDat <- as.data.frame(numVotes)
    rawDat$ModAgree <- apply(rawDat[,4:7],1,FUN = max)
    rawDat$Bifurcated <- apply(rawDat[,4:7],1,FUN = bifurc)
    rawDat$NewSuit <- apply(rawDat[,4:7],1,FUN = newSuitnoCurrent)
    rawDat <- rawDat[order(rawDat$Spp, rawDat$FuturePeriod),]
    rawDat$SuitDiff <- apply(rawDat[,c("FuturePeriod","Spp","Suitability","NewSuit")],1,stepSum)
    
    ##THE BELOW CUTOFFS COULD BE ADJUSTED
    rawDat$PeriodTraj <- ifelse(rawDat$SuitDiff > 0.2, "Improving",
                                ifelse(rawDat$SuitDiff > -0.2, "Stable",
                                       ifelse(rawDat$SuitDiff <= -0.2, "Declining", NA)))
    rawDat$Unsuit <- ifelse(rawDat$X > 0.5, "UNSUITABLE","OK")
    rawDat$modAgrClass <- ifelse(rawDat$ModAgree >= 0.8, "High",
                                 ifelse(rawDat$ModAgree >= 0.6, "Moderate","Low"))
    
    ##Feasibility of establishment using 2025 time period####
    Feas <- numVotes[numVotes$FuturePeriod == 2025,]
    Feas$FeasEstab <- Feas$`1`+Feas$`2`+(Feas$`3`*0.75) ###CONSTANT HERE COULD BE ADJUSTED
    Feas$NewSuit <- apply(Feas[,c(4:8)],1,newSuit)
    Feas$NewSuit <- round(Feas$NewSuit, digits = 0)
    Feas$NewSuit <- ifelse(Feas$NewSuit == 0, 1,
                           ifelse(Feas$NewSuit >= 4, "X", Feas$NewSuit))
    Feas$FofEClass <- ifelse(Feas$FeasEstab >= 0.8, "High", ##CUTOFFS COULD BE ADJUSTED
                             ifelse(Feas$FeasEstab >= 0.65, "Moderate",
                                    ifelse(Feas$FeasEstab >= 0.5, "Questionable","Not Feasible")))
    
    ####2055 data used for summary statistics####
    dat55 <- numVotes[numVotes$FuturePeriod == 2055, c(1,4:8)]
    dat55$Improve <- apply(dat55[,2:6],1,modDir2,direction = "Improve")
    dat55$Stable <- apply(dat55[,2:6],1,modDir2,direction = "Stable")
    dat55$Decline <- apply(dat55[,2:6],1,modDir2,direction = "Decline")
    dat55$Bifurc <- apply(dat55[,7:9],1,bifurcTrend)###test for bifurcation based on percent improve/decline
    dat55$NewSuit <- apply(dat55[,2:6],1,FUN = newSuitnoCurrent)
    dat55$SuitDiff <- round(dat55$Suitability - dat55$NewSuit, digits = 0)
    dat55$OverallTraj <- ifelse(dat55$SuitDiff >= 1, paste("Improving",dat55$Bifurc),
                                ifelse(dat55$SuitDiff >=0, paste("Stable",dat55$Bifurc),
                                       ifelse(dat55$SuitDiff >= -1, paste("Declining",dat55$Bifurc), paste("Strongly Declining",dat55$Bifurc))))
    dat55 <- dat55[order(dat55$Spp),]
    Feas <- Feas[order(Feas$Spp),]
    
    ###create summary output
    outputSum <- data.frame(Species = Feas$Spp, CurrentSuit = Feas$Suitability, CurrentSS = Feas$SSCurrent, Feasibility.Estab = Feas$FofEClass,
                            NewSuit = Feas$NewSuit, TrendMidRot = dat55$OverallTraj, Improve = dat55$Improve, Stable = dat55$Stable, Decline = dat55$Decline)
    outputSum$SiteNo <- Site
    outputSum$CurrentSuit[outputSum$CurrentSuit == 5] <- "X"
    
    #############Create raw output##################
    outputRaw <- rawDat[,c(1:3,4:8,11,13:15,10)]
    outputRaw[,4:9] <- round(outputRaw[,4:9], digits = 2)
    outputRaw[,4:7] <- outputRaw[,4:7]*100
    outputRaw$SiteNo <- Site
    
    list(outputRaw,outputSum)
  }###end of loop


write.csv(allOutput[[2]], "SummaryExample.csv")
write.csv(allOutput[[1]], "RawDataExample.csv")
















#===============================================================================
#OLD CODE
#         
#===============================================================================
treeShape <- function(temp){
  x <- temp[1,]
  for(i in 2:5){
    x <- cbind(x, temp[i,])
  }
  return(x)
}

Y1 <- Y1[-grep("Ensemble", Y1$GCM),]
Y1 <- Y1[-grep("rcp26", Y1$GCM),]
Y1 <- Y1[grep("2025", Y1$GCM),]

Y1 <- Y1[order(Y1$Latitude, Y1$Longitude),]
len <- length(Y1$PPT07)
Y1$ID1 <- rep(1:(len/90), each = 90)
Y1 <- Y1[Y1$Latitude >= 53,]


modDir <- function(x, direction){#columns 1,2,3,X,Currentsuit
  curr <- x[5]
  if(curr == 5 | is.na(curr)){
    curr <- 4
  }
  x[5] <- 0
  x <- c(0,x)
  improve <- sum(x[1:curr])
  same <- x[curr + 1]
  decline <- sum(x[(curr+2):length(x)])
  v <- c(improve,same,decline)
  names(v) <- c("Improve","Same","Decline")
  out <- paste(round((max(v)[1]*100), digits = 0),"% ",names(v)[v == max(v)[1]], sep = "")
  return(out[1])
}
names(FTS2)

## > names(FTS2)
##  [1] "BGC"        "Spp"        "Suit"       "Same"       "Imp1"       "Imp2"       "Dec1"       "Dec2"       "A1"        
## [10] "A2"         "A3"         "Not Suit"   "FSuit"      "Tbase"      "Trajectory"


str(SiteNo.suit)

FTS3 <- SiteNo.suit[ ,c('FuturePeriod', 'SiteNo','BGC', 'Spp', 'Suit', 'FSuit', 'Trajectory', 'Same', 'Imp1', 'Imp2', 'Dec1', 'Dec2', 'A1', 'A2', 'A3', 'Not Suit')]

FTS3

FTS3[, ]

str(FTS3)

setnames(FTS3, old = c('BGC','Suit', 'Imp1', 'Imp2', 'Dec1', 'Dec2', 'FSuit' ), new = c('Site Series','CurrentSuit', 'Improve1', 'Improve2', 'Decline1', 'Decline2', 'FutureSuit' ))

str(FTS3)

#FTS3$Trajectory <- as.numeric(FTS3$Trajectory)
FTS3[, 8:16] <- round(FTS3[, 8:16], digits=2)

SuitforMaps <- FTS3[,3:6]
write.csv(SuitforMaps, file = "PredSuit_AllMods_2085.csv")
write.csv(FTS3, file = "FTS3_2085.csv")
Suit.table <- Suit.table[grep("01|05|06", Suit.table$`Site Series`),]

SiteLocation = c("SiteNo","Latitude", "Longitude", "Elevation")
FTS.location = unique(Y1[,names(Y1) %in% c(SiteLocation)])
FTS4 <- merge(FTS.location, FTS3, by="SiteNo")
write.csv(FTS3, file = "SuitabilityNewMoriceSA.csv")

siteNames <- data.frame("Name" = c("Babine03","Babine05","Babine02","FultonFSR","Granisle","Hannay","Maxan","MoriceNew1","MoriceNew2"),
                        "SiteNo" = c(1,2,3,4,5,6,7,8,9))

Suit.table <- merge(siteNames, Suit.table, by = "SiteNo")
#===============================================================================
# Step 17: write formatted Tree species suitability table to file
#         
#===============================================================================
##Create examples for Bryce
Suit.table <- FTS3[order(FTS3$SiteNo, FTS3$`Site Series`,FTS3$Spp, FTS3$FuturePeriod),]
Suit.table <- Suit.table[,c(3,4,1,5,8:16)]
write.csv(Suit.table, file = "Suitability_9_StudyAreas.csv")

fname <- paste(fplot,"_TreeSuitTraject2",model,".csv",sep= "")

## write.csv(FTS3, file=paste(fplot,"_TreeSuitTraject2",model,".csv",sep=""))

#write.csv(FTS4, file=paste(fplot,"_TreeSuitTraject2",model,".csv",sep=""))   ## Need to output FTS4 rather than FTS3! 
write.csv(FTS4, file=paste(Tile.pick,"_TreeSuitTraject2",model,".csv",sep=""))   ## Need to output FTS4 rather than FTS3! 

## FTS4 contains Latitude, Longitude and Elevation, 
## whereas FTS3 does not! 


cat("\n\n","formatted Tree species suitability table is written to following file:", "\n\n",fname, "\n\n")


################################################################################

qsort <- function(x) {
  n <- length(x)
  if (n == 0) {
    x
  } else {
    p <- sample(n, 1)
    smaller <- foreach(y=x[-p], .combine=c) %:% when(y <= x[p]) %do% y
    larger  <- foreach(y=x[-p], .combine=c) %:% when(y >  x[p]) %do% y
    c(qsort(smaller), x[p], qsort(larger))
  }
}
qsort(runif(50))