#############################################################
#										#
#	R SCRIPT FOR MODELING PROJECTED GLOBAL MANGROVE		#
#		DISTRIBUTIONS UNDER CLIMATE CHANGE			#
#										#
#	NOAH CHARNEY and SYDNE RECORD, APRIL 2012					#
#										#
#############################################################
setwd('C:/Users/PhantPhant/Dropbox/mangroves')
source('mangrove_source.R')
source('kernelScale.R')
library(shapefiles)
library(rgdal)
library(raster)
library(cwhmisc)
library(ncf)
library(gbm)
library(MASS)

#############################################
# Individual species distribution models
###################################################

###########################################
### Generate pseudo absence data sets and site x species and environment matrices
###########################################

#Clear workspace and set working directory
rm(list=ls())
setwd("c:.../mangroves") # Replace ... with working directory

# Read in GBIF presence/absence and weights (# of records per grid cell) data for entire world
GBIFwt <- read.csv("GBIF_final_site-spp_matrix.csv", header=TRUE)
GBIFpa <- read.csv("GBIF_final_binary.csv", header=TRUE)

# Read in current environmental data (contains all cells in coastline) for entire world
currentEnvWorld <- read.csv("env-values_0m.csv", header=TRUE)

# Crop environment data matrices to exclude far poleward parts of globe 
ymin <- -5.2e6
ymax <- 5.2e6
currentEnvCrop <- currentEnvWorld[currentEnvWorld[,'y'] > ymin ,]
currentEnvCrop <- currentEnvCrop[currentEnvCrop[,'y'] < ymax ,]

# Create a correlation matrix for environmental variables for subsequent variable selection process.
cormat <- cor(currentEnvCrop[,4:24])
write.csv(cormat, file=".../cormat_current.csv")

# Plot coastline and occurrence data showing crop lines
x.coast <- currentEnvCrop[,'x']
y.coast <- currentEnvCrop[,'y']
xy <- cbind(x.coast,y.coast)
plot(xy[sample(1:(dim(currentEnvCrop)[1]),10000),],cex=.2)
points(GBIFpa[,'x'],GBIFpa[,'y'], col='red')
abline(h=c(ymin,ymax)) 

# Sum species occurrence data to determine what species have a total number of cells with occurrences >50
sumPA <- apply(GBIFpa, 2, sum)
species <- c("Avicennia_germinans","Ceriops_tagal","Laguncularia_racemosa", "Lumnitzera_littorea", "Lumnitzera_racemosa", "Rhizophora_apiculata","Rhizophora_mangle", "Rhizophora_mucronata", "Rhizophora_racemosa", "Rhizophora_stylosa", "Sonneratia_alba")
length(species)

########### Crop future environmental matrices 
rise0 <- read.csv("env-values_0m.csv", header=TRUE)
rise0Crop <- rise0[rise0[,'y'] > ymin ,]
rise0Crop <- rise0Crop[rise0Crop[,'y'] < ymax ,]
write.csv(rise0Crop, file="future0m.csv")

rise1 <- read.csv("env-values_1m.csv", header=TRUE)
rise1Crop <- rise1[rise1[,'y'] > ymin ,]
rise1Crop <- rise1Crop[rise1Crop[,'y'] < ymax ,]
write.csv(rise1Crop_noNA, file="future1m.csv")

rise3 <- read.csv("env-values_3m.csv", header=TRUE)
rise3Crop <- rise3[rise3[,'y'] > ymin ,]
rise3Crop <- rise3Crop[rise3Crop[,'y'] < ymax ,]
write.csv(rise3Crop_noNA, file="future3m.csv")

rise6 <- read.csv("env-values_6m.csv", header=TRUE)
rise6Crop <- rise6[rise6[,'y'] > ymin ,]
rise6Crop <- rise6Crop[rise6Crop[,'y'] < ymax ,]
write.csv(rise6Crop_noNA, file="future6m.csv")

################################
# Generate presence site x spp-env matrix for each species to be modelled
for(i in 1:length(species)){

	#get index
	sp_name <- species[i]
	
	# Make environment matrix for sites where species i is present
	pres_row_index <- GBIFpa[ GBIFpa[,sp_name] > 0,'coast.index']
	pres_env_rows <- currentEnvCrop[match(pres_row_index, currentEnvCrop[,'X']),]
	
	# Create a vector of ones (indicating presences) to append to pres_env_rows matrix
	pres_vec <- as.matrix((rep(1, length(pres_row_index))))
	colnames(pres_vec) <- sp_name
	# Create a vector of weights
	GBIFwt_pres_rows <- GBIFwt[ match(pres_row_index, GBIFwt[,'coast.index']),]
	wt_vec <- as.matrix(GBIFwt_pres_rows[,sp_name])
	colnames(wt_vec) <- "weights"
	pres_env_mat <- data.frame(cbind(pres_env_rows, pres_vec, wt_vec))
	
	# write csv file
	write.csv(pres_env_mat_noNA, file=paste("C:...",sp_name,"_pres.csv",sep='')) 
	
}

# Pseudo-absence 1: all collections in the occurrence data set are used as absences, even where a given species was collected (Phillips 2009 approach)	
abs1_env_rows <- currentEnvCrop[match(GBIFpa[,'coast.index'], currentEnvCrop[,'X']),] 
write.csv(abs1_env_rows, file=".../pseudoabs1.csv")

# Pseudo-absence 2: select 500 random points from the entire study area (coastline) with an equal weight of presences versus background data (Fitzpatrick 2011 approach)
random500 <- sample(currentEnvCrop[,'X'], 500, replace=FALSE) 
abs2_env_rows <- currentEnvCrop[match(random500, currentEnvCrop[,'X']),]
# write csv for pseudo abs 2
write.csv(abs2_env_rows, file="C:.../pseudoabs2.csv")

# Pseudo-absence 3: select 1000 random points from the entire study area (coastline) with an equal weight of presences versus background data (Fitzpatrick 2011 approach)
random1000 <- sample(currentEnvCrop[,'X'], 1000, replace=FALSE) 
abs3_env_rows <- currentEnvCrop[match(random1000, currentEnvCrop[,'X']),]
# write csv for pseudo abs 3
write.csv(abs3_env_rows, file=".../pseudoabs3.csv")

# Pseudo-absence 4: select 10000 random points from the entire study area (coastline) with an equal weight of presences versus background data (Fitzpatrick 2011 approach)
random10000 <- sample(currentEnvCrop[,'X'], 10000, replace=FALSE) 
abs4_env_rows <- currentEnvCrop[match(random10000, currentEnvCrop[,'X']),]
# write csv for pseudo abs 4
write.csv(abs4_env_rows, file=".../pseudoabs4.csv")

#############################################################################
# Perform GBM variable importance analysis
#############################################################################
##clear everything first
rm(list=ls())

library(BIOMOD)
setwd("...") # Replace ... with working directory containing empty files with species names and corresponding folders for outputs
species <- list.files()

# Loop through each species and pseudo-absence data set performing GBM to get variable importance
for(i in 1:length(species)){
sp.name <- species[i]
setwd("...") # Set directory to path with folders containing the site x spp and environment matrices you made in the previous step 

####### Pseudo-absence method 1 ######################
# Read in current species/environment matrix for presence localities
current <- read.csv(paste(".../",sp.name,"_pres.csv",sep=''), header=TRUE)

# Read in pseudo-absence data 1
pseudo1 <- read.csv(".../pseudoabs1.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo1)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo1)[1]))
colnames(abs_wts) <- "weights"
pseudo1 <- data.frame(cbind(pseudo1,abs, abs_wts))
data <- data.frame(rbind(current, pseudo1))

# Specify x,y coordinates for the sites
coords <- data[,2:3]

setwd(paste(".../",sp.name,"/pseudo1",sep=''))

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = data[,4:24], sp.name=sp.name)

system.time(
Models(GLM = FALSE,
       GBM = TRUE, No.trees = 3000, 
       GAM = FALSE,
       CTA = FALSE,
       ANN = FALSE, 
       SRE = FALSE, 
       FDA = FALSE, 
       MARS = FALSE, 
       RF = FALSE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Check variable importance output to see what variables to keep for the final models
    v1 <- VarImportance[[1]]
#####################################
rm(list=setdiff(ls(), c("current", "sp.name","v1", "species"))) 
setwd("...") # Set file path to folder containing site x species and environment matrices
####### Pseudo-absence method 2 ######################
# Read in pseudo-absence data 2
pseudo2 <- read.csv(".../pseudoabs2.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo2)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo2)[1]))
colnames(abs_wts) <- "weights"
pseudo2 <- data.frame(cbind(pseudo2,abs, abs_wts))
data <- data.frame(rbind(current, pseudo2))

# Specify x,y coordinates for the sites
coords <- data[,2:3]

setwd(paste(".../",sp.name,"/pseudo2",sep=''))

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = data[,4:24], sp.name=sp.name)

system.time(
Models(GLM = FALSE,
       GBM = TRUE, 
       GAM = FALSE, 
       CTA = FALSE, 
       ANN = FALSE, 
       SRE = FALSE, 
       FDA = FALSE, 
       MARS = FALSE, 
       RF = FALSE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Check variable importance output to see what variables to keep for the final models
    v2 <- VarImportance[[1]]
#####################################    
#####################################
rm(list=setdiff(ls(), c("current", "sp.name","v1","v2","species"))) 
setwd("...") # Set file path to folder containing site x species and environment matrices
####### Pseudo-absence method 2 ######################
# Read in pseudo-absence data 3
pseudo3 <- read.csv(".../pseudoabs3.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo3)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo3)[1]))
colnames(abs_wts) <- "weights"
pseudo3 <- data.frame(cbind(pseudo3,abs, abs_wts))
data <- data.frame(rbind(current, pseudo3))

# Specify x,y coordinates for the sites
coords <- data[,2:3]

#setwd(paste("c:/users/phantphant/dropbox/mangroves/sydne/sdms/",sp.name,"/pseudo1",sep=''))
setwd(paste("C:/Documents and Settings/srecord.FAS_DOMAIN/My Documents/Sydne/Mangroves/sdms_2012/Variable_Importance/",sp.name,"/pseudo3",sep=''))

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = data[,4:24], sp.name=sp.name)

system.time(
Models(GLM = FALSE,# TypeGLM = "poly", Test = "BIC", 
       GBM = TRUE, No.trees = 3000, 
       GAM = FALSE,# Spline = 3, 
       CTA = FALSE, #CV.tree = 50, 
       ANN = FALSE, #CV.ann = 5, 
       SRE = FALSE, #quant=0.025, 
       FDA = FALSE, 
       MARS = FALSE, 
       RF = FALSE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Check variable importance output to see what variables to keep for the final models
    v3 <- VarImportance[[1]]
#####################################
####### Pseudo-absence method 4 ######################
rm(list=setdiff(ls(), c("current", "sp.name","v1","v2","v3","species"))) 

#setwd("C:/users/phantphant/dropbox/mangroves/sydne/sdms/")
setwd("c:/documents and settings/srecord.fas_domain/my documents/my dropbox/mangroves/sydne/sdms")
# Read in pseudo-absence data 4
pseudo4 <- read.csv("data_for_sdms/pseudoabs4.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo4)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo4)[1]))
colnames(abs_wts) <- "weights"
pseudo4 <- data.frame(cbind(pseudo4,abs, abs_wts))
data <- data.frame(rbind(current, pseudo4))

# Specify x,y coordinates for the sites
coords <- data[,2:3]

#setwd(paste("c:/users/phantphant/dropbox/mangroves/sydne/sdms/",sp.name,"/pseudo1",sep=''))
setwd(paste("C:/Documents and Settings/srecord.FAS_DOMAIN/My Documents/Sydne/Mangroves/sdms_2012/Variable_Importance/",sp.name,"/pseudo4",sep=''))

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = data[,4:24], sp.name=sp.name)

system.time(
Models(GLM = FALSE,# TypeGLM = "poly", Test = "BIC", 
       GBM = TRUE, No.trees = 3000, 
       GAM = FALSE,# Spline = 3, 
       CTA = FALSE, #CV.tree = 50, 
       ANN = FALSE, #CV.ann = 5, 
       SRE = FALSE, #quant=0.025, 
       FDA = FALSE, 
       MARS = FALSE, 
       RF = FALSE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Check variable importance output to see what variables to keep for the final models
    v4 <- VarImportance[[1]]
######

output <- rbind(v1, v2, v3, v4)
output_mu <- apply(output, 2, mean)
modVars <- names(sort(output_mu, decreasing=T))
write.csv(modVars, file=paste("C:/Documents and Settings/srecord.FAS_DOMAIN/My Documents/Sydne/Mangroves/sdms_2012/Variable_Importance/",sp.name,"/",sp.name,"TopVars.csv",sep=''))

rm(list=setdiff(ls(), c("current","sp.name","species"))) 
}

####################################################
# Loop through each species with >50 occurrence records at 2.5 minute resolution and fit species distribution models using BIOMOD
####################################################

##clear everything first
rm(list=ls())

library(BIOMOD)
species <- c("Avicennia_bicolor", "Avicennia_germinans", "Avicennia_marina", "Ceriops_australis", "Ceriops_tagal", "Kandelia_candel", "Laguncularia_racemosa", "Lumnitzera_littorea", "Lumnitzera_racemosa", "Nypa_fruticans", "Rhizophora_apiculata", "Rhizophora_harrisonii", "Rhizophora_mangle", "Rhizophora_mucronata", "Rhizophora_racemosa", "Rhizophora_stylosa", "Sonneratia_alba", "Sonneratia_caseolaris")

path <- # specify folder directory path
setwd(path)
explanatory_mat <- read.csv("explanatory_mat.csv",header=TRUE) # explanatory_mat contains the list of explanatory variables for each species as selected by the variable importance procedure

for(i in 1:length(species)){
sp.name <- species[i]
# Explanatory variables selected by variable importance procedure
explanatory_cols <- as.matrix(explanatory_mat[i,3:7])

####### Pseudo-absence method 1 ######################
# Read in current species/environment matrix for presence localities
current <- read.csv(paste("data_for_sdms/",sp.name,"_pres.csv",sep=''), header=TRUE)

# Read in pseudo-absence data 1
pseudo1 <- read.csv("data_for_sdms/pseudoabs1.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo1)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo1)[1]))
colnames(abs_wts) <- "weights"
pseudo1 <- data.frame(cbind(pseudo1,abs, abs_wts))
data <- data.frame(rbind(current, pseudo1))

# Specify x,y coordinates for the sites
coords <- data[,2:3]

# Explanatory variables selected by variable importance procedure
explanatory <- data[,explanatory_cols
setwd(paste("...",sp.name,"/pseudo1",sep='')) # Set directory ... to write BIOMOD output files to

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = explanatory, sp.name=sp.name)

system.time(
Models(GLM = TRUE, TypeGLM = "poly", Test = "BIC", 
       GBM = TRUE, No.trees = 3000, 
       GAM = TRUE, Spline = 3, 
       CTA = TRUE, CV.tree = 50, 
       ANN = TRUE, CV.ann = 5, 
       SRE = TRUE, quant=0.025, 
       FDA = TRUE, 
       MARS = TRUE, 
       RF = TRUE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Get evaluation results
Eval1 <- Evaluation.results.TSS[1]
#####################################
# Remove junk objects
rm(list=setdiff(ls(), c("explanatory_mat","current", "sp.name","Eval1", "explanatory_cols", "path", "species"))) 


####### Pseudo-absence method 2 ######################
# Read in pseudo-absence data 2
pseudo2 <- read.csv("data_for_sdms/pseudoabs2.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo2)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo2)[1]))
colnames(abs_wts) <- "weights"
pseudo2 <- data.frame(cbind(pseudo2,abs, abs_wts))
data <- data.frame(rbind(current, pseudo2))

# Specify x,y coordinates for the sites
coords <- data[,2:3]

# Explanatory variables
explanatory <- data[,explanatory_cols]

setwd(paste("...",sp.name,"/pseudo1",sep='')) # Set directory ... to write BIOMOD output files to 

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = explanatory, sp.name=sp.name)

system.time(
Models(GLM = TRUE, TypeGLM = "poly", Test = "BIC", 
       GBM = TRUE, No.trees = 3000, 
       GAM = TRUE, Spline = 3, 
       CTA = TRUE, CV.tree = 50, 
       ANN = TRUE, CV.ann = 5, 
       SRE = TRUE, quant=0.025, 
       FDA = TRUE, 
       MARS = TRUE, 
       RF = TRUE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Get evaluation results
Eval2 <- Evaluation.results.TSS[1]
#####################################    
#####################################
# Remove junk objects
rm(list=setdiff(ls(), c("explanatory_mat","path","species","current", "sp.name","Eval1","Eval2","explanatory_cols"))) 

####### Pseudo-absence method 3 ######################
# Read in pseudo-absence data 3
pseudo3 <- read.csv("data_for_sdms/pseudoabs3.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo3)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo3)[1]))
colnames(abs_wts) <- "weights"
pseudo3 <- data.frame(cbind(pseudo3,abs, abs_wts))
data <- data.frame(rbind(current, pseudo3))

# Specify x,y coordinates for the sites
coords <- data[,2:3]
explanatory <- data[,explanatory_cols]
setwd(paste("...",sp.name,"/pseudo1",sep='')) # Set directory ... to write BIOMOD output files to

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = explanatory, sp.name=sp.name)

system.time(
Models(GLM = TRUE, TypeGLM = "poly", Test = "BIC", 
       GBM = TRUE, No.trees = 3000, 
       GAM = TRUE, Spline = 3, 
       CTA = TRUE, CV.tree = 50, 
       ANN = TRUE, CV.ann = 5, 
       SRE = TRUE, quant=0.025, 
       FDA = TRUE, 
       MARS = TRUE, 
       RF = TRUE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

#Get evaluation results
Eval3 <- Evaluation.results.TSS[1] 

#####################################
####### Pseudo-absence method 4 ######################
# Remove junk objects
rm(list=setdiff(ls(), c("explanatory_mat","path","species","current", "sp.name","Eval1","Eval2","Eval3", "explanatory_cols","species"))) 

# Read in pseudo-absence data 4
pseudo4 <- read.csv("data_for_sdms/pseudoabs4.csv", header=TRUE)
abs <- as.matrix(rep(0, dim(pseudo4)[1]))
colnames(abs) <- sp.name
abs_wts <- as.matrix(rep(1, dim(pseudo4)[1]))
colnames(abs_wts) <- "weights"
pseudo4 <- data.frame(cbind(pseudo4,abs, abs_wts))
data <- data.frame(rbind(current, pseudo4))

# Specify x,y coordinates for the sites
coords <- data[,2:3]
explanatory  <- data[,explanatory_cols]

setwd(paste("...",sp.name,"/pseudo1",sep='')) # Set directory ... to write BIOMOD output files to

# Specify the initial state
Initial.State(Response = data[,44], Explanatory = explanatory, sp.name=sp.name)

system.time(
Models(GLM = TRUE, TypeGLM = "poly", Test = "BIC", 
       GBM = TRUE, No.trees = 3000, 
       GAM = TRUE, Spline = 3, 
       CTA = TRUE, CV.tree = 50, 
       ANN = TRUE, CV.ann = 5, 
       SRE = TRUE, quant=0.025, 
       FDA = TRUE, 
       MARS = TRUE, 
       RF = TRUE, 
       NbRunEval = 10,
       DataSplit = 70, Yweights=data[,45], Roc = T, Optimized.Threshold.Roc = T,
       Kappa = T, TSS=T, KeepPredIndependent = T, VarImport=10, NbRepPA=0)
)

# Get Evaluation results
Eval4 <- Evaluation.results.TSS[1]
#######

mu1 <- mean(as.matrix(Eval1[1][[1]][3]))
mu2 <- mean(as.matrix(Eval2[1][[1]][3]))
mu3 <- mean(as.matrix(Eval3[1][[1]][3]))
mu4 <- mean(as.matrix(Eval4[1][[1]][3]))

output <- as.vector(c(mu1, mu2,mu3,mu4))
write.csv(output, file=paste("C:/Documents and Settings/srecord.FAS_DOMAIN/My Documents/Sydne/Mangroves/sdms_2012/reduced/",sp.name,"eval.csv",sep=''))

rm(list=setdiff(ls(), c("explanatory_mat","path","species","current")))

}

###################################
# Project ensemble forecasts of models for the 500 random pseudo-absence approach that had the highest True Skill Statistic score
###################################
# Clear workspace to remove junks
rm(list=ls())

species <- c("Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
explanatory_mat <- read.csv(".../explanatory_mat.csv",header=TRUE) # load matrix of explanatory variables selected for each species according to the variable importance procedure
# Load species x environment matrices for each future sea level rise scenario (0m or no change, 1m, 3m, 4m)
norise0 <- read.csv(".../current0m.csv",header=TRUE)
rise0 <- read.csv(".../future0m.csv",header=TRUE)
rise1 <- read.csv(".../future1m.csv",header=TRUE)
rise3 <- read.csv(".../future3m.csv",header=TRUE)
rise6 <- read.csv(".../future6m.csv",header=TRUE)

for(i in 1:length(species)){

sp.name <- species[i]
explanatory_cols_current <- data.frame(subset(explanatory_mat, V1==sp.name[i]))
explanatory_cols_current <- as.matrix(explanatory_cols_current[,3:7])
currentEnv <- norise0[,explanatory_cols_current]
futureEnv0 <- rise0[,explanatory_cols_current]
futureEnv1 <- rise1[,explanatory_cols_current]
futureEnv3 <- rise3[,explanatory_cols_current]
futureEnv6 <- rise6[,explanatory_cols_current]

setwd(paste("...",sp.name,"/pseudo2",sep='')) # specify directory containing fitted BIOMOD outputs from pseudoabsence method 2
load(paste(sp.name,"_run.Rdata",sep=''))

# Current projections and ensemble forecast

Projection(Proj=currentEnv, Proj.name="Current",GLM = T, GBM = T, GAM = T,
 	CTA = T, ANN = T, SRE = T, FDA =T, MARS = T,
        RF = T, BinRoc = F, BinKappa = F, BinTSS = T, FiltRoc = F,
        FiltKappa = F, FiltTSS = F, repetition.models=T)

 
Ensemble.Forecasting(ANN = T, CTA = T, GAM = T, GBM = T, GLM = T, MARS = T,
        FDA = T, RF = T, SRE = T, Proj.name="Current", weight.method="TSS",
        decay = "proportional", PCA.median = F, binary = T, bin.method = "TSS", Test = TRUE,
        repetition.models=F, final.model.out=T)


# Future 0m sea-level rise
Projection(Proj=futureEnv0, Proj.name="Future0m",GLM = T, GBM = T, GAM = T,
 	CTA = T, ANN = T, SRE = T, FDA =T, MARS = T,
        RF = T, BinRoc = F, BinKappa = F, BinTSS = T, FiltRoc = F,
        FiltKappa = F, FiltTSS = F, repetition.models=T)


  
Ensemble.Forecasting(ANN = T, CTA = T, GAM = T, GBM = T, GLM = T, MARS = T,
        FDA = T, RF = T, SRE = T, Proj.name="Future0m", weight.method="TSS",
        decay = "proportional", PCA.median = F, binary = T, bin.method = "TSS", Test = TRUE,
        repetition.models=F, final.model.out=T)


# Future 1m sea-level rise
Projection(Proj=futureEnv1, Proj.name="Future1m",GLM = T, GBM = T, GAM = T,
 	CTA = T, ANN = T, SRE = T, FDA =T, MARS = T,
        RF = T, BinRoc = F, BinKappa = F, BinTSS = T, FiltRoc = F,
        FiltKappa = F, FiltTSS = F, repetition.models=T)


  
Ensemble.Forecasting(ANN = T, CTA = T, GAM = T, GBM = T, GLM = T, MARS = T,
        FDA = T, RF = T, SRE = T, Proj.name="Future1m", weight.method="TSS",
        decay = "proportional", PCA.median = F, binary = T, bin.method = "TSS", Test = TRUE,
        repetition.models=F, final.model.out=T)


# Future 3m sea-level rise
Projection(Proj=futureEnv3, Proj.name="Future3m",GLM = T, GBM = T, GAM = T,
 	CTA = T, ANN = T, SRE = T, FDA =T, MARS = T,
        RF = T, BinRoc = F, BinKappa = F, BinTSS = T, FiltRoc = F,
        FiltKappa = F, FiltTSS = F, repetition.models=T)


  
Ensemble.Forecasting(ANN = T, CTA = T, GAM = T, GBM = T, GLM = T, MARS = T,
        FDA = T, RF = T, SRE = T, Proj.name="Future3m", weight.method="TSS",
        decay = "proportional", PCA.median = F, binary = T, bin.method = "TSS", Test = TRUE,
        repetition.models=F, final.model.out=T)


# Future  m sea-level rise
Projection(Proj=futureEnv0, Proj.name="Future0m",GLM = T, GBM = T, GAM = T,
 	CTA = T, ANN = T, SRE = T, FDA =T, MARS = T,
        RF = T, BinRoc = F, BinKappa = F, BinTSS = T, FiltRoc = F,
        FiltKappa = F, FiltTSS = F, repetition.models=T)


  
Ensemble.Forecasting(ANN = T, CTA = T, GAM = T, GBM = T, GLM = T, MARS = T,
        FDA = T, RF = T, SRE = T, Proj.name="Future0m", weight.method="TSS",
        decay = "proportional", PCA.median = F, binary = T, bin.method = "TSS", Test = TRUE,
        repetition.models=F, final.model.out=T)
}

################################
# Obtain True Skill Statistic results from projected model outputs
################################
# Clear workspace
rm(list=ls())

species <- c("Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
# Create empty matrix for TSS score results for each ensemble projection
TSSmat <- matrix(NA, nrow=length(species), ncol=5)
rownames(TSSmat) <- species
colnames(TSSmat) <- c("Current", "Future_0m", "Future_1m", "Future_3m", "Future_6m")

# loop through species to get TSS results and insert those results into the output matrix 'TSSmat'
for(i in 1:length(species)){
setwd(paste(".../",species[i],sep=''))

# Load current ensemble forecast results
load("Proj.Current/consensus_current_results")
temp <- consensus_Current_results[1][1:6][[1]][4]
TSSmat[i,1] <- temp$test.results[6]
#Remove junk objects
rm("consensus_Current_results")
rm("temp")

# Load future 0m (no change in sea leve) ensemble forecast results
load("Proj.Future0m/consensus_Future0m_results")
temp <- consensus_Future0m_results[1][1:6][[1]][4]
TSSmat[i,2] <- temp$test.results[6]
#Remove junk objects
rm("consensus_Future0m_results")
rm("temp")

# Load future 1m sea-level rise ensemble forecast results
load("Proj.Future1m/consensus_Future1m_results")
temp <- consensus_Future1m_results[1][1:6][[1]][4]
TSSmat[i,3] <- temp$test.results[6]
#Remove junk objects
rm("consensus_Future1m_results")
rm("temp")

# Load future 3m sea-level rise ensemble forecast results
load("Proj.Future3m/consensus_Future3m_results")
temp <- consensus_Future3m_results[1][1:6][[1]][4]
TSSmat[i,4] <- temp$test.results[6]
#Remove junk objects
rm("consensus_Future3m_results")
rm("temp")

# Load future 6m sea-level rise ensemble forecast results
load("Proj.Future6m/consensus_Future6m_results")
temp <- consensus_Future6m_results[1][1:6][[1]][4]
TSSmat[i,5] <- temp$test.results[6]
#Remove junk objects
rm("consensus_Future6m_results")
rm("temp")

}

# Write TSSmat results to a csv file, set directory path ...
write.csv(TSSmat, file=".../TSS_results.csv")

#####################################################
# Obtain binary results from projections for mapping 
######################################################
# clear workspace
rm(list=ls())
path <- "..." # specify directory containing projected model outputs

# Make output matrices that contain x,y coordinates for each current/future sea-level rise scenario. 

species <- c("Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
currentEnv <- read.csv(".../current0m.csv",header=TRUE)
currentCoords <- as.matrix(currentEnv[,c("x","y")]) 
current_out <- matrix(NA, nrow=dim(currentCoords)[1], ncol=length(species)+2)
colnames(current_out) <- c("x","y","Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
current_out[,1:2] <- currentCoords

futureEnv0 <- read.csv(".../future0m.csv",header=TRUE)
future0Coords <- as.matrix(futureEnv0[,c("x","y")]) 
future0_out <- matrix(NA, nrow=dim(future0Coords)[1], ncol=length(species)+2)
colnames(future0_out) <- c("x","y","Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
future0_out[,1:2] <- future0Coords

futureEnv1 <- read.csv(".../future1m.csv",header=TRUE)
future1Coords <- as.matrix(futureEnv1[,c("x","y")]) 
future1_out <- matrix(NA, nrow=dim(future1Coords)[1], ncol=length(species)+2)
colnames(future1_out) <- c("x","y","Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
future1_out[,1:2] <- future1Coords

futureEnv3 <- read.csv(".../future3m.csv",header=TRUE)
future3Coords <- as.matrix(futureEnv3[,c("x","y")]) 
future3_out <- matrix(NA, nrow=dim(future3Coords)[1], ncol=length(species)+2)
colnames(future3_out) <- c("x","y","Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
future3_out[,1:2] <- future3Coords

futureEnv6 <- read.csv(".../future6m.csv",header=TRUE)
future6Coords <- as.matrix(futureEnv6[,c("x","y")]) 
future6_out <- matrix(NA, nrow=dim(future6Coords)[1], ncol=length(species)+2)
colnames(future6_out) <- c("x","y","Lumnitzera_littorea", "Rhizophora_apiculata", "Rhizophora_mucronata", "Rhizophora_racemosa", "Sonneratia_alba", "Rhizophora_stylosa", "Lumnitzera_racemosa", "Ceriops_tagal", "Avicennia_marina", "Rhizophora_mangle", "Laguncularia_racemosa", "Avicennia_germinans")
future6_out[,1:2] <- future6Coords

# Loop through each modeled species and get binary results from projection. Note that continuous model outputs are converted to binary using the TSS scrore.
for(i in 1:12){
# Set working directory to match output for each species
setwd(paste(path,species[i],sep=''))

# Load results of ensemble projected for current climate (0m change in sea-level)
load("Proj.Current/Total_consensus_Current_bin")
current_out[,i+2] <- as.matrix(Total_consensus_Current_Bin[,,6])
# Remove junk objects
rm("Total_consensus_Current_Bin")

# Load results of ensemble projected for future climate and 0m change in sea-level
load("Proj.Future0m/Total_consensus_Future0m_bin")
future0_out[,i+2] <- as.matrix(Total_consensus_Future0m_Bin[,,6])
#remove junk objects
rm("Total_consensus_Future0m_Bin")

# Load results of ensemble projected for future climate and 1m change in sea-level
load("Proj.Future1m/Total_consensus_Future1m_bin")
future1_out[,i+2] <- as.matrix(Total_consensus_Future1m_Bin[,,6])
# Remove junk objects
rm("Total_consensus_Future1m_Bin")

# Load results of ensemble projected for future climate and 3m change in sea-level
load("Proj.Future3m/Total_consensus_Future3m_bin")
future3_out[,i+2] <- as.matrix(Total_consensus_Future3m_Bin[,,6])
# remove junk objects
rm("Total_consensus_Future3m_Bin")

# Load results of ensemble projected for future climate and 6m change in sea-level
load("Proj.Future6m/Total_consensus_Future6m_bin")
future6_out[,i+2] <- as.matrix(Total_consensus_Future6m_Bin[,,6])
#remove junk objects
rm("Total_consensus_Future6m_Bin")

}

# Write outputs to csv files. Specify directory ...
write.csv(current_out, file=".../current_binary.csv")
write.csv(future0_out, file=".../future0_binary.csv")
write.csv(future1_out, file=".../future1_binary.csv")
write.csv(future3_out, file=".../future3_binary.csv")
write.csv(future6_out, file=".../future6_binary.csv")

# Quick plots of outputs
par(mfrow=c(2,3))
level.plot(as.matrix(current_out[,7]), current_out[,1:2])
level.plot(as.matrix(future0_out[,7]), future0_out[,1:2])
level.plot(as.matrix(future1_out[,7]), future1_out[,1:2])
level.plot(as.matrix(future3_out[,7]), future3_out[,1:2])
level.plot(as.matrix(future6_out[,7]), future6_out[,1:2])
####################################################
# Community Distribution Models
########################################################

###############################################
#	PROCESSED DATA VIEWING STEPS 
#		(earlier processing steps below)
###############################################
#
#	Make Richness data
#

names(GBIF[-(1:4)])
GBIF <- read.csv('GBIF_FINAL_site-spp_matrix.csv')
GBIF.binary <- GBIF
for(i in 1:(dim(GBIF)[1])){
	for(j in 5:34){
		GBIF.binary[i,j] <- min(GBIF[i,j],1)
	}
}	
write.csv(GBIF.binary[,-1],'GBIF_FINAL_binary.csv')

#	read in richness

GBIF.binary <- read.csv('GBIF_FINAL_binary.csv')
head(GBIF.binary)
richness <- apply(GBIF.binary[,5:34],1,sum)
richmat <- cbind(GBIF.binary[,2:4],richness)
#hist(richmat[,'richness'])
sum(richness>=3)
richrows<-richmat[,'richness'] >= 3
#points(richmat[richrows,'x'],richmat[richrows,'y'],col='orange')
#plot(GBIF.binary[,'x'],GBIF.binary[,'y'],cex=.4)
richmat.3sp.bin <- richmat[richmat[,'richness'] >=3,]
richmat.3sp.bin[,'richness'] <- 1
dim(richmat.3sp.bin)
write.csv(richmat.3sp.bin,'richness_3sp_bin.csv')

delete zeros
richmat<-richmat[richmat[,'richness']>0 , ]
write.csv(richmat,'richmat.csv')

#	make richness data matrix for modelling

num.absences <- 2000
env <-read.csv('env-values_0m_updated.csv')
names(env)[1] <- 'coast.index'
ymin <- -5.2e6
ymax <- 5.2e6
env <- subset(env,y > ymin)
env <- subset(env,y < ymax)

all.coast.index <- env[,'coast.index']
pres.index <- richmat[,'coast.index']
abs.index <- all.coast.index[!all.coast.index %in% pres.index]
abs.sample.index <- sample(abs.index, num.absences)

pres.data <- cbind(richness=richmat[,'richness'],
			env[ match(pres.index,env[,'coast.index']) , ]
			)
abs.data <- cbind(richness=0,
			env[ match(abs.sample.index,env[,'coast.index']) , ]
		)
abs.data.full <-cbind(richness=0,
			env[ match(abs.index,env[,'coast.index']) , ]
		)

pres.data <- pres.data[ !is.na(apply(pres.data,1,sum)),]
abs.data <- abs.data[ !is.na(apply(abs.data,1,sum)),]

all.data <- rbind(pres.data,abs.data)
all.data.full <- rbind(pres.data,abs.data.full)

# train data on eastern half of world, then test against west
# and visa versa

all.data.west <- subset(all.data,x<2.5e6)
all.data.east <- subset(all.data,x>2.5e6)
all.data.eastwest <- all.data

all.data.full.west <- subset(all.data.full,x<2.5e6)
all.data.full.east <- subset(all.data.full,x>2.5e6)
all.data.full.eastwest <- all.data.full



all.data <- all.data.east

###############################################3
#
#	MODEL SELECTION 
#			for richness
#
#################################################
#		
#
first.env <-5
last.env <- 25
data.names <- names(all.data)

pred.names <- data.names[5:25]
pred.names <- names(env)[4:24]

#messing around
#plot(data.holdout[,'richness']~predicted.gbm)
#plot(data.holdout[,'richness']~predicted.pois)

#input
model.names = c('gbm','gbm.red','gbm.drop','glm','step','p_red','dropcorr')
abs.nums = c(0,500,1000,2000,10000, length(abs.index)) #c(0,0)
iter.names <- paste( c(0,0,500,500,1000,1000,2000,2000,10000,10000, length(abs.index),length(abs.index)) , rep( c('E-W','W-E') , 5) ) 

LL.array <- array(dim = c(10,length(iter.names),length(model.names)),	dimnames=list(NULL,iter.names,model.names))

num.array <- LL.array
var.array <- array(dim=c(10,length(iter.names),length(pred.names),length(model.names)), 
			dimnames=list(NULL,iter.names,pred.names,model.names))
source('kernelscale.r')
system.time(
for(i in 1:10){
for(j in 1:6){

	pres.data <- cbind(richness=richmat[,'richness'],
			env[ match(pres.index,env[,'coast.index']) , ]
			)
	abs.data.full <-cbind(richness=0,
				env[ match(abs.index,env[,'coast.index']) , ]
			)
	data.full <- rbind(pres.data,abs.data.full)
	
	if(abs.nums[j]==0){
		data.mat <- pres.data
	}else{
		abs.sample.index <- sample(abs.index, abs.nums[j])
		abs.data <- cbind(richness=0,
				env[ match(abs.sample.index,env[,'coast.index']) , ]
			)
		
		pres.data <- pres.data[ !is.na(apply(pres.data,1,sum)),]
		abs.data <- abs.data[ !is.na(apply(abs.data,1,sum)),]
	
		data.mat <- rbind(pres.data,abs.data)
	}
	data.west <- subset(data.mat,x<2.5e6)
	data.east <- subset(data.mat,x>2.5e6)

	data.full.west <- subset(data.full,x<2.5e6)
	data.full.east <- subset(data.full,x>2.5e6)

	#run multimod.func.rescale() for doing resaled data, only works on GBM and GLM full
	#run multimod.func() for examining different ways of subsetting variables, does not rescale
	multimod.EW <- multimod.func.rescale(pred.names = pred.names, 
				data.train=data.east,
				data.holdout=data.full.west,
				cellsize=5e5, n.trees=3000)
	multimod.WE <- multimod.func.rescale(pred.names = pred.names, 
				data.train=data.west,
				data.holdout=data.full.east,
				cellsize=5e5, n.trees=3000)
	
	LL.array[i,2*j-1,] <-multimod.EW$LL.vec
	LL.array[i,2*j,] <-multimod.WE$LL.vec
	num.array[i,2*j-1,] <-multimod.EW$num.vec
	num.array[i,2*j,] <-multimod.WE$num.vec
	var.array[i,2*j-1,,] <-multimod.EW$var.mat
	var.array[i,2*j,,] <-multimod.WE$var.mat
	
	save(var.array,file='var_array_10.R')
	save(LL.array,file='LL_array_10.R')
	save(num.array,file='num_array_10.R')
	print(paste('i=',i))
	print(paste('j=',j))

}}
)



#view model comparisons
load('LL_array_8.r')
load('num_array_4.r')
load('var_array_4.r')
#play around with boxplotting different slices...
boxplot(t(LL.array[1,c(2,4,6,8),]+LL.array[2,c(2,4,6,8),]+LL.array[2,c(2,4,6,8),])/3)

###############################################
#
#	CONSTRUCT FINAL MODEL
#			& Predict to future
#
###############################################


env0m <- read.csv('env-values_0m_updated.csv')
env1m <- read.csv('env-values_1m_updated.csv')
env3m <- read.csv('env-values_3m_updated.csv')
env6m <- read.csv('env-values_6m_updated.csv')

#
#give future scenarios env columns the same as current column names
#
env <- env3m
curr.names <- names(env)[1:25]
fut.names <- names(env)[c(1:4,26:45, 25)]
fut.m <- matrix(nrow=dim(env)[1],ncol=25, dimnames=list(NULL,curr.names))
for(i in 1:25){
	cur.temp <- curr.names[i]
	fut.temp <- fut.names[i]
	fut.m[,cur.temp] <- env[,fut.temp]
}
dimnames(fut.m)
not.na.rows <- !is.na(apply(fut.m,1,sum))
fut.m <- data.frame(fut.m)
fut.3m <- fut.m[not.na.rows,]
fut.3m <- subset(fut.3m,y<5.2e6)
fut.3m <- subset(fut.3m,y>-5.2e6)

cur.0m <- env0m[,1:25]
not.na.rows <- !is.na(apply(cur.0m,1,sum))
cur.0m <- cur.0m[not.na.rows,]
cur.0m <- data.frame(cur.0m)
cur.0m <- subset(cur.0m,y<5.2e6)
cur.0m <- subset(cur.0m,y>-5.2e6)


#difference in percent of NAs in different coastlines
#	(expect less NAs with more sea level rise, b/c bioclim layers
#	are only inland
dim(fut.m)
dim(fut.0m)[1]/dim(env0m)[1]
dim(fut.1m)[1]/dim(env1m)[1]
dim(fut.3m)[1]/dim(env3m)[1]
dim(fut.6m)[1]/dim(env6m)[1]

#
#	make final model
#

#Loop through 10 times
system.time(
for(q in 1:100){
num.absences <- 2000
all.coast.index <- cur.0m[,'X']
abs.index <- all.coast.index[!all.coast.index %in% pres.index]
abs.sample.index <- sample(abs.index, num.absences)
abs.data <- cbind(richness=0,
			cur.0m[ match(abs.sample.index,cur.0m[,'X']) , ]
		)

richmat <- read.csv('richmat.csv')
pres.index <- richmat[,'coast.index']
pres.data <- cbind(richness=richmat[,'richness'],
			cur.0m[ match(pres.index,cur.0m[,'X']) , ]
			)

all.data <- rbind(pres.data,abs.data)


final.model.gbm <- gbm(make.formula(curr.names[5:25]), 
				data=all.data, 
				distribution='poisson', 
				n.trees=3000)

#
#	Make final predictions
#

predicted.gbm.cur0m <- predict.gbm(final.model.gbm,
					newdata=cur.0m, 
					type='response',
					n.trees=3000)
x <- cur.0m[,'x']
y <- cur.0m[,'y']
pred.gbm.mat.cur <- cbind(x,y,value=predicted.gbm.cur0m)
predicted.gbm.fut3m <- predict.gbm(final.model.gbm,
					newdata=fut.3m, 
					type='response',
					n.trees=3000)
x <- fut.3m[,'x']
y <- fut.3m[,'y']
pred.gbm.mat.fut3m <- cbind(x,y,value=predicted.gbm.fut3m)
#save
write.csv(pred.gbm.mat.cur, paste('predicted_GBM_2000abs_cur0m_iter',q,'.csv',sep=''))
write.csv(pred.gbm.mat.fut3m, paste('predicted_GBM_2000abs_fut3m_iter', q, '.csv', sep=''))
print(q)
}#end for loop over q
)


####################################################
#
#	MAP PREDICTIONS
#		& make figures
#
######################################################
GBIF <- read.csv('GBIF_FINAL_binary.csv')
curr <-read.csv('sydne/SDMs/output2/current_binary.csv')
fut0 <-read.csv('sydne/SDMs/output2/future0_binary.csv')
fut1 <-read.csv('sydne/SDMs/output2/future1_binary.csv')
fut3 <-read.csv('sydne/SDMs/output2/future3_binary.csv')
fut6 <-read.csv('sydne/SDMs/output2/future6_binary.csv')

################################
#
#	FIGURE 1 
#		histogram of latitudes
#
###################################

x <- curr[,'x']
y <- curr[,'y']
latlon<-project(xy=cbind(x,y),inv=TRUE,proj='+proj=goode')
hist.out<-hist(abs(latlon[,2]))
barplot(hist.out$counts*4.318,
		main='', xlab='Absolute value of latitude', 
		ylab='length of coastline (km)'
#		names.arg=hist.out$mids
		)
axis(side=1,at=c(0,6,12,18,24),labels=c(0,10,20,30,40))

#
#
#	INDIVIDUAL SPECIES MODELS
#	
#

#
#Crop individual predictions to meaningful values
#
sp.totals <- apply(GBIF,2,sum)
sp.names<-names(sp.totals[order(sp.totals,decreasing=TRUE)][5:16])
load('cropsmat.r')

for(i in 1:length(sp.names)){
	sp.temp <- sp.names[i]
	x.min <- crops.mat[sp.temp,'xmin']
	x.max <- crops.mat[sp.temp,'xmax']
	curr[curr[,'x'] > x.max , sp.temp] <- 0
	curr[curr[,'x'] < x.min , sp.temp] <- 0
	fut0[fut0[,'x'] > x.max , sp.temp] <- 0
	fut0[fut0[,'x'] < x.min , sp.temp] <- 0
	fut1[fut1[,'x'] > x.max , sp.temp] <- 0
	fut1[fut1[,'x'] < x.min , sp.temp] <- 0
	fut3[fut3[,'x'] > x.max , sp.temp] <- 0
	fut3[fut3[,'x'] < x.min , sp.temp] <- 0
	fut6[fut6[,'x'] > x.max , sp.temp] <- 0
	fut6[fut6[,'x'] < x.min , sp.temp] <- 0
}

#
#	Make coast points for plotting
#
#make coords on a 0 to 1 scale for plotting on same axes as "image"
#'image' function spills centers first and last cells at 0 and 1, so they spill over
# half a cell below 0 and above 1, so using range.x and denominator is ok

coast.plot.coords <- cbind(x=curr[,'x'],y=curr['y'])
x<-coast.plot.coords[,'x']
y<-coast.plot.coords[,'y']
coast.plot.coords[,'x'] <- (x + 2e7)/(4e7)
coast.plot.coords[,'y'] <- (y + 5.2e6)/(1.04e7)

#
#	make summary table for species
#

scenarios <- c('curr','fut0','fut1','fut3','fut6')
col.names <- c('latmin','latmax','latmean','latsd','sum','diff','perc.diff')
summary.array <- array(dim=c(length(scenarios),length(sp.names), length(col.names)),
				dimnames=list(scenarios,sp.names,col.names)
		)

for(i in 1:length(scenarios)){
	scenario.temp <- get(scenarios[i])
	xy <- cbind(scenario.temp[,'x'],scenario.temp[,'y'])
	latlon<-project(xy=xy,inv=TRUE,proj='+proj=goode')
	scenario.temp <- cbind(lon=latlon[,1],lat=latlon[,2],scenario.temp)

	curr.temp <- curr
	xy.curr <- cbind(curr.temp[,'x'],curr.temp[,'y'])
	latlon.cur<-project(xy=xy.curr,inv=TRUE,proj='+proj=goode')
	curr.temp <- cbind(lon=latlon.cur[,1],lat=latlon.cur[,2],curr.temp)



	for(j in 1:length(sp.names)){
			sp.temp <- sp.names[j]
			x.min <- crops.mat[sp.temp,'xmin']
			x.max <- crops.mat[sp.temp,'xmax']
			scenario.temp[scenario.temp[,'x']<x.min,sp.temp] <- 0
			scenario.temp[scenario.temp[,'x']>x.max,sp.temp] <- 0
			curr.temp[curr.temp[,'x']<x.min,sp.temp] <- 0
			curr.temp[curr.temp[,'x']>x.max,sp.temp] <- 0
			curr.rows.temp <- curr.temp[,sp.temp] >0

			rows.temp <- scenario.temp[,sp.temp] > 0
			summary.array[scenarios[i],sp.temp,'latmin'] <- min(abs(scenario.temp[rows.temp,'lat']))
			summary.array[scenarios[i],sp.temp,'latmax'] <- max(abs(scenario.temp[rows.temp,'lat']))
			summary.array[scenarios[i],sp.temp,'latmean'] <- mean(abs(scenario.temp[rows.temp,'lat']))
			summary.array[scenarios[i],sp.temp,'latsd'] <- sd(abs(scenario.temp[rows.temp,'lat']))
			summary.array[scenarios[i],sp.temp,'sum'] <- sum(rows.temp)
			summary.array[scenarios[i],sp.temp,'diff'] <- summary.array[scenarios[i],sp.temp,'sum'] - 
											summary.array[scenarios[1],sp.temp,'sum']
			summary.array[scenarios[i],sp.temp,'perc.diff'] <- summary.array[scenarios[i],sp.temp,'diff'] / 
											summary.array[scenarios[1],sp.temp,'sum']

			#make maps of change for each species summarized at coarse scale
			if(TRUE){

				x <- scenario.temp[rows.temp,'x']
				y <- scenario.temp[rows.temp,'y']
				value <- scenario.temp[rows.temp, sp.temp]
				pred.mat <- cbind(x,y,value)

				x <- curr.temp[curr.temp[,sp.temp] > 0 ,'x']
				y <- curr.temp[curr.temp[,sp.temp] >  0 ,'y']
				value <- curr.temp[curr.temp[,sp.temp]>0, sp.temp]
				curr.mat <- cbind(x,y,value)

				cellsize <- 1000000
				pred.worldmap <- make.worldmap(pred.mat, cellsize=cellsize,
						max.x=2e7, min.x = -2e7, max.y = 5.2e6, min.y = -5.2e6
				)
				curr.worldmap <- make.worldmap(curr.mat,cellsize=cellsize,
						max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6
				)

				if(i==1 & j==1){
					world.change.array <- array(dim=c(length(scenarios),length(sp.names),dim(curr.worldmap)),
						dimnames=list(scenarios,sp.names,NULL,NULL)
					)
				}

				#amount of change in species coastline occupancy
				world.change.array[scenarios[i],sp.names[j],,] <- pred.worldmap - curr.worldmap



#				image.change(world.out.fut=pred.worldmap,
#					world.out.cur=curr.worldmap,	
#					coast.plot.coords=coast.plot.coords,
#					main=paste(sp.temp,scenarios[i])
#				)
#
#				waitReturn(ask=TRUE)
			}
	}
}

save(world.change.array, file='world_change_array.r')

#
#	Turn summary array into table
#

summary.mat <- matrix(nrow=length(scenarios)*length(sp.names), ncol=(length(col.names)+3),
				dimnames=list(NULL,c('species','sp.code','scenario',col.names))
			)
split.names <- strsplit(sp.names,'_')
sp.codes <- sp.names
for(i in 1:length(sp.names)){
	sp.codes[i] <- paste(substr(split.names[[i]],1,2),collapse='')
}
summary.mat[,'species'] <- rep.collapse(sp.names,length(scenarios))
summary.mat[,'sp.code'] <- rep.collapse(sp.codes,length(scenarios))
summary.mat[,'scenario'] <- rep(scenarios,length(sp.names))
for(i in 1:(dim(summary.array)[1])){
	for(j in 1:(dim(summary.array)[2])){
		for(k in 1:(dim(summary.array)[3])){
			f.row <- (1:(dim(summary.mat)[1]))[ summary.mat[,'species' ]==dimnames(summary.array)[[2]][j] &
									summary.mat[,'scenario']==dimnames(summary.array)[[1]][i] ]
			summary.mat[ f.row , dimnames(summary.array)[[3]][k] ] <- summary.array[i,j,k]
		}
	}
}
write.csv(summary.mat,'summary_mat.csv')
save(summary.array,file='summary_array.R')

################################
#
#	FIGURE 2
#		boxplot of latitude of each species
#
#############################

plot(1,type='n',xlim=c(0,85),ylim= c(0,55),xaxt='n',xlab='',ylab='Latitude')

#this is the focal x position for each vertical bar
f.x <- 1:60 + rep.collapse( seq(0,22,2),5)
breaks <- c(0,get.steps(9,400,1.3))
scenarios.num <- rep(c('c',0,1,3,6),12)

for(i in 1:60){ 
	diff.temp <- as.numeric(summary.mat[i,'perc.diff'])*100
	if(diff.temp < breaks[2] & diff.temp > -breaks[2]){
		f.color <- gray(0.7)
	}else{
		if(diff.temp > 0){
			f.color <- pos.col[ ceiling(length(pos.col)*findInterval(diff.temp,breaks) / length(breaks))]
		}else{
			f.color <- neg.col[  ceiling(length(neg.col)*findInterval(-diff.temp,breaks) / length(breaks))]
		}
	}
	lines(x=c(f.x[i],f.x[i]),y=as.numeric(summary.mat[i,c('latmin','latmax')]),lwd=2, col=f.color)
	lines(x=c(f.x[i],f.x[i]),y=c( as.numeric(summary.mat[i,'latmean']) + as.numeric(summary.mat[i,'latsd']),
						as.numeric(summary.mat[i,'latmean']) - as.numeric(summary.mat[i,'latsd'])
					),
					lwd=8,
					col=f.color)
	lines(x=c(f.x[i]-.4,f.x[i]+.4),y=rep(as.numeric(summary.mat[i,'latmean']),2),lwd=4)
	text(x=f.x[i], y = as.numeric(summary.mat[i,'latmax'])+4, scenarios.num[i])
	if(summary.mat[i,'scenario']=='curr'){
		mtext(substr(summary.mat[i,'sp.code'],1,6), side=1, at=f.x[i],adj=0) #las=2 for perpendicular labels
	}
}

plot.scalebar(breaks,neg.col,pos.col,units='percent change')



###############################
#
#	FIGURE 3
#		plot change in amount of coastline for each species
#
########################################

load('world_change_array.r')
#this is to get color range for each species, so they are on same scale

#Mapping color breaks
		breaks <- c(0,get.steps(5,3000,2))

		#make color scheme
		neg.col=c(rainbow(200,alpha=1)[130:110],
				gray(.8)
			)
		pos.col=c(	gray(.8),
				rainbow(200,alpha=1)[30:10]
			)

map.aspect.r <- dim(world.change.array)[4] / dim(world.change.array)[3]
plot.mat <- t(matrix(1:12,nrow=3, ncol=4))
plot.widths <- rep(11,3) #set plot width here in cm (11 cm works on this little comp)
plot.heights <- rep(plot.widths[1]*map.aspect.r,4)
par(mai=c(0,0,0,0))
layout(mat=plot.mat,widths=lcm(plot.widths),heights=lcm(plot.heights),respect=TRUE)
#resize plot device window to full screen here

for(j in 1:12){ 
#12 species
	for(i in 4){
	#scenario 4 is future 3m 

		#get focal species and plot coastlie
		f.spec <- dimnames(world.change.array)[[2]][j]
		f.code <- paste(substr(strsplit(f.spec,'_')[[1]],1,2),collapse='')
		plot(coast.plot.coords, cex=.1, 
			#main=paste(scenarios[i],f.spec), 
			xlim = c(0,1), 
			ylim=c(0,1),bg='transparent',xaxt='n',yaxt='n',xlab='',ylab=''
			)
		#make seperate matrices for areas where coastline 
		#	increases(pos), decreases(neg), or stays near zero(neut)
		neut.abs <- pos.abs <- neg.abs <- world.change.array[i,f.spec,,]
		pos.abs[pos.abs < 0] <- NA
		neg.abs[neg.abs >= 0] <- NA
		neut.abs <- abs(neut.abs)

		neut.abs[neut.abs >= breaks[2]] <- NA
		if(sum(neut.abs,na.rm=TRUE)!=0){		
			image(neut.abs,col=gray(.8),add=TRUE)
		}
		for(k in 3:length(breaks)){
			pos.layer.temp <- pos.abs
			pos.layer.temp[pos.layer.temp < breaks[k-1] | pos.layer.temp >= breaks[k]] <- NA
			neg.layer.temp <- neg.abs
			neg.layer.temp[neg.layer.temp > -breaks[k-1] | neg.layer.temp <= -breaks[k]] <- NA					

			pos.col.temp <- pos.col[ ceiling(length(pos.col) * k/length(breaks))]
			neg.col.temp <- neg.col[ ceiling(length(neg.col) * k/length(breaks))]
			if(sum(pos.layer.temp,na.rm=TRUE)!=0){		
				image(pos.layer.temp, col=pos.col.temp,add=TRUE)
			}
			if(sum(neg.layer.temp,na.rm=TRUE)!=0){		
				image(neg.layer.temp, col=neg.col.temp,add=TRUE)
			}	

		}
		text(-0.04,.12,f.code,cex=3,pos=4)
	}
}

################################################
#
#	FIGURE 4
#		Richness Plots
#
################################################
 	
xrange <- 4e7
yrange <- 10.4e6
map.aspect.r <- yrange / xrange
plot.mat <- matrix(1:4,nrow=4, ncol=1)
plot.widths <- rep(14,1)#set plot width here in cm
plot.heights <- rep(plot.widths[1]*map.aspect.r,4)
#par(mai=c(.1,.1,.5,.1))
par(mai=c(0,0,0,0))
layout(mat=plot.mat,widths=lcm(plot.widths),heights=lcm(plot.heights),respect=TRUE)
#resize plot device window to full screen here



#
#	plot GBIF richness, observed
#
cellsize <- 5e5

richmat <- read.csv('richmat.csv')
GBIF.worldmap <- make.worldmap(richmat,cellsize=cellsize,value.name='richness',
				max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6
			)

plot(coast.plot.coords, cex=.2,xlim = c(0,1), ylim=c(0,1),bg='transparent',
			xaxt='n',yaxt='n',xlab='',
			ylab='a) Observed Species Density',
			main='')
text(-0.02,.9,'a',cex=3,pos=4)

x <- c(0,0)
y <- c(23.5,-23.5)
latlon <- project(xy=cbind(x,y),inv=FALSE,proj='+proj=goode')
latlon[,2] <- c(2616008, -2616008)
abline(h=(0.5+c(latlon[1,2]/2e7,-latlon[1,2]/2e7)),lty=2,col=gray(.5))


GBIF.col=c(rainbow(200,alpha=.8)[40:1])
breaks <- get.steps(1,300,2)
for(i in 1:length(breaks)){
	layer.temp <- GBIF.worldmap
	layer.temp[layer.temp < breaks[i] | layer.temp >= breaks[i+1]] <- NA
	image(layer.temp,col=GBIF.col[ ceiling(length(GBIF.col)*i/length(breaks)) ],add=TRUE)
}

#plot.scalebar(breaks,GBIF.col,GBIF.col,GBIF.col[1],'Species richness x grid cells',TRUE)


#	
#	Richness
#		plot richness based on sum of future projections of 12 species


rich.curr <- cbind(x=curr[,'x'],y=curr[,'y'],value = apply(curr[,4:15],1,sum,na.rm=TRUE))
rich.fut3 <- cbind(x=fut3[,'x'],y=fut3[,'y'],value = apply(fut3[,4:15],1,sum,na.rm=TRUE))

x<-rich.curr[,'x']
y<-rich.curr[,'y']
latlons <- project(xy=cbind(x,y),inv=TRUE,proj='+proj=goode')
weighted.mean(x=abs(latlons[,2]),w=rich.curr[,'value'])

x<-rich.fut3[,'x']
y<-rich.fut3[,'y']
latlons <- project(xy=cbind(x,y),inv=TRUE,proj='+proj=goode')
weighted.mean(x=abs(latlons[,2]),w=rich.fut3[,'value'])


curr.rich.worldmap<-make.worldmap(value.mat = rich.curr ,cellsize=cellsize,
			max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6)
fut3.rich.worldmap<-make.worldmap(value.mat = rich.fut3 ,cellsize=cellsize,
			max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6)
plot(coast.plot.coords, cex=.2, xlim = c(0,1), ylim=c(0,1),bg='transparent',
		xaxt='n',yaxt='n',xlab='',
		ylab='b) Composite Projection',
		main='')
text(-0.02,.9,'b',cex=3,pos=4)


x <- c(0,0)
y <- c(23.5,-23.5)
latlon<-project(xy=cbind(x,y),inv=FALSE,proj='+proj=goode')
latlon[,2] <- c(2616008, -2616008)
#WTF is going on with project()???? why does it sometimes give a warning
#but other times, with the EXACT SAME CODE, works fine?????????
abline(h=(0.5+c(latlon[1,2]/2e7,-latlon[1,2]/2e7)),lty=2,col=gray(.5))

#make color scheme
neg.col=c(rainbow(200,alpha=.8)[130:110],
		'transparent'
		)
pos.col=c(	'transparent',
		rainbow(200,alpha=.8)[30:10]
		)
world.change <- (fut3.rich.worldmap - curr.rich.worldmap)/mean(curr.rich.worldmap,na.rm=TRUE)
#hist(world.change)
max(abs(world.change[world.change!=0]),na.rm=TRUE)
#
#	Set BREAKS here for colors for rest of the 3 richness plots
#
breaks <- c(0,0,get.steps(.01,30,2))

#make seperate matrices for areas where coastline 
#increases(pos), decreases(neg),
neut.abs <- pos.abs <- neg.abs <- world.change
		pos.abs[pos.abs < 0] <- NA
		neg.abs[neg.abs >= 0] <- NA
		neut.abs <- abs(neut.abs)

		neut.abs[neut.abs >= breaks[2]] <- NA
		if(sum(neut.abs,na.rm=TRUE)!=0){		
			image(neut.abs,col='transparent',add=TRUE)
		}
	for(k in 3:length(breaks)){
		pos.layer.temp <- pos.abs
		pos.layer.temp[pos.layer.temp < breaks[k-1] | pos.layer.temp >= breaks[k]] <- NA
		neg.layer.temp <- neg.abs
		neg.layer.temp[neg.layer.temp > -breaks[k-1] | neg.layer.temp <= -breaks[k]] <- NA

		pos.col.temp <- pos.col[ ceiling(length(pos.col) * k/length(breaks))]
		neg.col.temp <- neg.col[ ceiling(length(neg.col) * k/length(breaks))]
		if(sum(pos.layer.temp,na.rm=TRUE)!=0){		
			image(pos.layer.temp, col=pos.col.temp,add=TRUE)
		}
		if(sum(neg.layer.temp,na.rm=TRUE)!=0){		
			image(neg.layer.temp, col=neg.col.temp,add=TRUE)
		}

	}

#plot.scalebar(breaks*100,neg.col,pos.col,'transparent',units='Relative change in species density (%)')
#

#
#	plot predicted richness based on binary models
#

cellsize <- 5e5
source('kernelscale.r')
rich.curr <- read.csv('current_binary.csv')
rich.fut3 <- read.csv('future3_binary.csv')
head(rich.curr)

x<-rich.curr[,'x']
y<-rich.curr[,'y']
latlons <- project(xy=cbind(x,y),inv=TRUE,proj='+proj=goode')
weighted.mean(x=abs(latlons[,2]),w=rich.curr[,'richness'])

x<-rich.fut3[,'x']
y<-rich.fut3[,'y']
latlons <- project(xy=cbind(x,y), inv=TRUE, proj='+proj=goode')
weighted.mean(x=abs(latlons[,2]),w=rich.fut3[,'richness'])

curr.rich.worldmap<-make.worldmap(value.mat = rich.curr ,cellsize=cellsize,
			value.name='richness',
			max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6)
fut3.rich.worldmap<-make.worldmap(value.mat = rich.fut3 ,cellsize=cellsize,
			value.name='richness',
			max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6)
plot(coast.plot.coords, cex=.2, xlim = c(0,1), ylim=c(0,1),bg='transparent',
		xaxt='n',yaxt='n',xlab='',
		ylab='c) Binary Projection', 
		main='')
text(-0.02,.9,'c',cex=3,pos=4)


x <- c(0,0)
y <- c(23.5,-23.5)
latlon <- project(xy=cbind(x,y),inv=FALSE,proj='+proj=goode')
latlon[,2] <- c(2616008, -2616008)
abline(h=(0.5+c(latlon[1,2]/2e7,-latlon[1,2]/2e7)),lty=2,col=gray(.5))

#make color scheme
neg.col=c(rainbow(200,alpha=.8)[130:110],
		'transparent'
		)
pos.col=c(	'transparent',
		rainbow(200,alpha=.8)[30:10]
		)
world.change <- (fut3.rich.worldmap - curr.rich.worldmap)/mean(curr.rich.worldmap,na.rm=TRUE)
#hist(world.change)
max(abs(world.change[world.change!=0]),na.rm=TRUE)

#make seperate matrices for areas where coastline 
#increases(pos), decreases(neg),
neut.abs <- pos.abs <- neg.abs <- world.change
		pos.abs[pos.abs < 0] <- NA
		neg.abs[neg.abs >= 0] <- NA
		neut.abs <- abs(neut.abs)

		neut.abs[neut.abs >= breaks[2]] <- NA
		if(sum(neut.abs,na.rm=TRUE)!=0){		
			image(neut.abs,col='transparent',add=TRUE)
		}
	for(k in 3:length(breaks)){
		pos.layer.temp <- pos.abs
		pos.layer.temp[pos.layer.temp < breaks[k-1] | pos.layer.temp >= breaks[k]] <- NA
		neg.layer.temp <- neg.abs
		neg.layer.temp[neg.layer.temp > -breaks[k-1] | neg.layer.temp <= -breaks[k]] <- NA

		pos.col.temp <- pos.col[ ceiling(length(pos.col) * k/length(breaks))]
		neg.col.temp <- neg.col[ ceiling(length(neg.col) * k/length(breaks))]
		if(sum(pos.layer.temp,na.rm=TRUE)!=0){		
			image(pos.layer.temp, col=pos.col.temp,add=TRUE)
		}
		if(sum(neg.layer.temp,na.rm=TRUE)!=0){		
			image(neg.layer.temp, col=neg.col.temp,add=TRUE)
		}

	}

#plot.scalebar(breaks*4.318,neg.col,pos.col,'transparent',units='richness times change in coastline (km)')

#
#	plot predicted richness based on GBM continuous model
#

cellsize <- 5e5
rich.curr <- read.csv('predicted_GBM_2000abs_cur0m.csv')
rich.fut3 <- read.csv('predicted_GBM_2000abs_fut3m.csv')

x<-rich.curr[,'x']
y<-rich.curr[,'y']
latlons <- project(xy=cbind(x,y),inv=TRUE,proj='+proj=goode')
weighted.mean(x=abs(latlons[,2]),w=rich.curr[,'value'])
#16.8

x<-rich.fut3[,'x']
y<-rich.fut3[,'y']
latlons <- project(xy=cbind(x,y),inv=TRUE,proj='+proj=goode')
weighted.mean(x=abs(latlons[,2]),w=rich.fut3[,'value'])
#17.2

curr.rich.worldmap<-make.worldmap(value.mat = rich.curr ,cellsize=cellsize,
			max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6)
fut3.rich.worldmap<-make.worldmap(value.mat = rich.fut3 ,cellsize=cellsize,
			max.x=2e7, min.x = -2e7 ,max.y=5.2e6, min.y=-5.2e6)
plot(coast.plot.coords, cex=.2, xlim = c(0,1), ylim=c(0,1),bg='transparent',
		xaxt='n',yaxt='n',xlab='',
		ylab='d) Continuous Projection',
		main='')
text(-0.02,.9,'d',cex=3,pos=4)


x <- c(0,0)
y <- c(23.5,-23.5)
latlon<-project(xy=cbind(x,y),inv=FALSE,proj='+proj=goode')
latlon[,2] <- c(2616008, -2616008)
abline(h=(0.5+c(latlon[1,2]/2e7,-latlon[1,2]/2e7)),lty=2,col=gray(.5))

#make color scheme
neg.col=c(rainbow(200,alpha=.8)[130:110],
		'transparent'
		)
pos.col=c(	'transparent',
		rainbow(200,alpha=.8)[30:10]
		)
world.change <- (fut3.rich.worldmap - curr.rich.worldmap)/mean(curr.rich.worldmap,na.rm=TRUE)

#hist(world.change)
max(abs(world.change[world.change!=0]),na.rm=TRUE)

#make seperate matrices for areas where coastline 
#increases(pos), decreases(neg),
neut.abs <- pos.abs <- neg.abs <- world.change
		pos.abs[pos.abs < 0] <- NA
		neg.abs[neg.abs >= 0] <- NA
		neut.abs <- abs(neut.abs)

		neut.abs[neut.abs >= breaks[2]] <- NA
		if(sum(neut.abs,na.rm=TRUE)!=0){		
			image(neut.abs,col='transparent',add=TRUE)
		}
	for(k in 3:length(breaks)){
		pos.layer.temp <- pos.abs
		pos.layer.temp[pos.layer.temp < breaks[k-1] | pos.layer.temp >= breaks[k]] <- NA
		neg.layer.temp <- neg.abs
		neg.layer.temp[neg.layer.temp > -breaks[k-1] | neg.layer.temp <= -breaks[k]] <- NA

		pos.col.temp <- pos.col[ ceiling(length(pos.col) * k/length(breaks))]
		neg.col.temp <- neg.col[ ceiling(length(neg.col) * k/length(breaks))]
		if(sum(pos.layer.temp,na.rm=TRUE)!=0){		
			image(pos.layer.temp, col=pos.col.temp,add=TRUE)
		}
		if(sum(neg.layer.temp,na.rm=TRUE)!=0){		
			image(neg.layer.temp, col=neg.col.temp,add=TRUE)
		}

	}

#plot.scalebar(breaks*4.318,neg.col,pos.col,'transparent',units='richness times change in coastline (km)')