####Mega_sdm_habs####

#??megasdm
# 
library(devtools)
#install_github("brshipley/megaSDM", force = TRUE)
#Using Bio-ORACLE variables and MARSPEC bathymetry only (for now)
#Environmental layers have been pre-processed: removed all depth cells below 500m 

getwd()
#setwd("C:/Users/Desktop/Bogi/HAB_SDMs")

# 1. Importing and processing Bio-ORACLE layers.
library(raster)


# 1.2 After exporting, save and exit R studio, come back, and import the reprojected layers

#Import reprojected rasters ; Stack according to scenario and time frame ; Assign layer names 

#PRESENT ENVIRONMENT
rPresentCurrentVelocityMax <- raster("reproj_Present_CurrentVelocityMax.tif")
rPresentCurrentVelocityMean <- raster("reproj_CurrentVelocityMean.tif")
rPresentCurrentVelocityMin <- raster("reproj_CurrentVelocityMin.tif")
rPresentCurrentVelocityRange <- raster("reproj_CurrentVelocityRange.tif")
rPresentSalinityMax <- raster("reproj_SalinityMax.tif")
rPresentSalinityMean <- raster("reproj_SalinityMean.tif")
rPresentSalinityMin <- raster("reproj_SalinityMin.tif")
rPresentSalinityRange <- raster("reproj_SalinityRange.tif")
rPresentTemperatureMax <- raster("reproj_TemperatureMax.tif")
rPresentTemperatureMean <- raster("reproj_TemperatureMean.tif")
rPresentTemperatureMin <- raster("reproj_TemperatureMin.tif")
rPresentTemperatureRange <- raster("reproj_TemperatureRange.tif")
rjBathymetry <- raster("reproj_Bathymetry.tif")

projection(rPresentCurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentCurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentCurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentCurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentSalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentSalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentSalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentSalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentTemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentTemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(rPresentTemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"

projection(rjBathymetry) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"

env_present <- stack(rPresentCurrentVelocityMax,rPresentCurrentVelocityMean,rPresentCurrentVelocityMin,
                     rPresentCurrentVelocityRange,rPresentSalinityMax,rPresentSalinityMean,
                     rPresentSalinityMin,rPresentSalinityRange,rPresentTemperatureMax,rPresentTemperatureMean,
                     rPresentTemperatureMin,rPresentTemperatureRange,rjBathymetry)
names(env_present)
names(env_present) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_present)


# #Checking raster correlation plots
# library("ENMTools")
# raster.cor.plot(env_present)

# 


#2050_26 ENVIRONMENT
r2050_26CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2050R26.tif")
r2050_26CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2050R26.tif")
r2050_26CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2050R26.tif")
r2050_26CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2050R26.tif")
r2050_26SalinityMax <- raster("reproj_SalinityMax2050R26.tif")
r2050_26SalinityMean <- raster("reproj_SalinityMean2050R26.tif")
r2050_26SalinityMin <- raster("reproj_SalinityMin2050R26.tif")
r2050_26SalinityRange <- raster("reproj_SalinityRange2050R26.tif")
r2050_26TemperatureMax <- raster("reproj_TemperatureMax2050R26.tif")
r2050_26TemperatureMean <- raster("reproj_TemperatureMean2050R26.tif")
r2050_26TemperatureMin <- raster("reproj_TemperatureMin2050R26.tif")
r2050_26TemperatureRange <- raster("reproj_TemperatureRange2050R26.tif")

projection(r2050_26CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_26TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2050_26 <- stack(r2050_26CurrentVelocityMax,r2050_26CurrentVelocityMean,r2050_26CurrentVelocityMin,
                     r2050_26CurrentVelocityRange,r2050_26SalinityMax,r2050_26SalinityMean,
                     r2050_26SalinityMin,r2050_26SalinityRange,r2050_26TemperatureMax,r2050_26TemperatureMean,
                     r2050_26TemperatureMin,r2050_26TemperatureRange,rjBathymetry)
names(env_2050_26)
names(env_2050_26) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2050_26)


#2050_45 ENVIRONMENT
r2050_45CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2050R45.tif")
r2050_45CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2050R45.tif")
r2050_45CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2050R45.tif")
r2050_45CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2050R45.tif")
r2050_45SalinityMax <- raster("reproj_SalinityMax2050R45.tif")
r2050_45SalinityMean <- raster("reproj_SalinityMean2050R45.tif")
r2050_45SalinityMin <- raster("reproj_SalinityMin2050R45.tif")
r2050_45SalinityRange <- raster("reproj_SalinityRange2050R45.tif")
r2050_45TemperatureMax <- raster("reproj_TemperatureMax2050R45.tif")
r2050_45TemperatureMean <- raster("reproj_TemperatureMean2050R45.tif")
r2050_45TemperatureMin <- raster("reproj_TemperatureMin2050R45.tif")
r2050_45TemperatureRange <- raster("reproj_TemperatureRange2050R45.tif")

projection(r2050_45CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_45TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2050_45 <- stack(r2050_45CurrentVelocityMax,r2050_45CurrentVelocityMean,r2050_45CurrentVelocityMin,
                     r2050_45CurrentVelocityRange,r2050_45SalinityMax,r2050_45SalinityMean,
                     r2050_45SalinityMin,r2050_45SalinityRange,r2050_45TemperatureMax,r2050_45TemperatureMean,
                     r2050_45TemperatureMin,r2050_45TemperatureRange,rjBathymetry)
names(env_2050_45)
names(env_2050_45) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2050_45)


#2050_60 ENVIRONMENT
r2050_60CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2050R60.tif")
r2050_60CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2050R60.tif")
r2050_60CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2050R60.tif")
r2050_60CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2050R60.tif")
r2050_60SalinityMax <- raster("reproj_SalinityMax2050R60.tif")
r2050_60SalinityMean <- raster("reproj_SalinityMean2050R60.tif")
r2050_60SalinityMin <- raster("reproj_SalinityMin2050R60.tif")
r2050_60SalinityRange <- raster("reproj_SalinityRange2050R60.tif")
r2050_60TemperatureMax <- raster("reproj_TemperatureMax2050R60.tif")
r2050_60TemperatureMean <- raster("reproj_TemperatureMean2050R60.tif")
r2050_60TemperatureMin <- raster("reproj_TemperatureMin2050R60.tif")
r2050_60TemperatureRange <- raster("reproj_TemperatureRange2050R60.tif")

projection(r2050_60CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_60TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2050_60 <- stack(r2050_60CurrentVelocityMax,r2050_60CurrentVelocityMean,r2050_60CurrentVelocityMin,
                     r2050_60CurrentVelocityRange,r2050_60SalinityMax,r2050_60SalinityMean,
                     r2050_60SalinityMin,r2050_60SalinityRange,r2050_60TemperatureMax,r2050_60TemperatureMean,
                     r2050_60TemperatureMin,r2050_60TemperatureRange,rjBathymetry)
names(env_2050_60)
names(env_2050_60) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2050_60)

#2050_85 ENVIRONMENT
r2050_85CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2050R85.tif")
r2050_85CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2050R85.tif")
r2050_85CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2050R85.tif")
r2050_85CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2050R85.tif")
r2050_85SalinityMax <- raster("reproj_SalinityMax2050R85.tif")
r2050_85SalinityMean <- raster("reproj_SalinityMean2050R85.tif")
r2050_85SalinityMin <- raster("reproj_SalinityMin2050R85.tif")
r2050_85SalinityRange <- raster("reproj_SalinityRange2050R85.tif")
r2050_85TemperatureMax <- raster("reproj_TemperatureMax2050R85.tif")
r2050_85TemperatureMean <- raster("reproj_TemperatureMean2050R85.tif")
r2050_85TemperatureMin <- raster("reproj_TemperatureMin2050R85.tif")
r2050_85TemperatureRange <- raster("reproj_TemperatureRange2050R85.tif")

projection(r2050_85CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2050_85TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2050_85 <- stack(r2050_85CurrentVelocityMax,r2050_85CurrentVelocityMean,r2050_85CurrentVelocityMin,
                     r2050_85CurrentVelocityRange,r2050_85SalinityMax,r2050_85SalinityMean,
                     r2050_85SalinityMin,r2050_85SalinityRange,r2050_85TemperatureMax,r2050_85TemperatureMean,
                     r2050_85TemperatureMin,r2050_85TemperatureRange,rjBathymetry)
names(env_2050_85)
names(env_2050_85) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2050_85)


#2100_26 ENVIRONMENT
r2100_26CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2100R26.tif")
r2100_26CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2100R26.tif")
r2100_26CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2100R26.tif")
r2100_26CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2100R26.tif")
r2100_26SalinityMax <- raster("reproj_SalinityMax2100R26.tif")
r2100_26SalinityMean <- raster("reproj_SalinityMean2100R26.tif")
r2100_26SalinityMin <- raster("reproj_SalinityMin2100R26.tif")
r2100_26SalinityRange <- raster("reproj_SalinityRange2100R26.tif")
r2100_26TemperatureMax <- raster("reproj_TemperatureMax2100R26.tif")
r2100_26TemperatureMean <- raster("reproj_TemperatureMean2100R26.tif")
r2100_26TemperatureMin <- raster("reproj_TemperatureMin2100R26.tif")
r2100_26TemperatureRange <- raster("reproj_TemperatureRange2100R26.tif")

projection(r2100_26CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_26TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2100_26 <- stack(r2100_26CurrentVelocityMax,r2100_26CurrentVelocityMean,r2100_26CurrentVelocityMin,
                     r2100_26CurrentVelocityRange,r2100_26SalinityMax,r2100_26SalinityMean,
                     r2100_26SalinityMin,r2100_26SalinityRange,r2100_26TemperatureMax,r2100_26TemperatureMean,
                     r2100_26TemperatureMin,r2100_26TemperatureRange,rjBathymetry)
names(env_2100_26)
names(env_2100_26) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2100_26)


#2100_45 ENVIRONMENT
r2100_45CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2100R45.tif")
r2100_45CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2100R45.tif")
r2100_45CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2100R45.tif")
r2100_45CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2100R45.tif")
r2100_45SalinityMax <- raster("reproj_SalinityMax2100R45.tif")
r2100_45SalinityMean <- raster("reproj_SalinityMean2100R45.tif")
r2100_45SalinityMin <- raster("reproj_SalinityMin2100R45.tif")
r2100_45SalinityRange <- raster("reproj_SalinityRange2100R45.tif")
r2100_45TemperatureMax <- raster("reproj_TemperatureMax2100R45.tif")
r2100_45TemperatureMean <- raster("reproj_TemperatureMean2100R45.tif")
r2100_45TemperatureMin <- raster("reproj_TemperatureMin2100R45.tif")
r2100_45TemperatureRange <- raster("reproj_TemperatureRange2100R45.tif")

projection(r2100_45CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_45TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2100_45 <- stack(r2100_45CurrentVelocityMax,r2100_45CurrentVelocityMean,r2100_45CurrentVelocityMin,
                     r2100_45CurrentVelocityRange,r2100_45SalinityMax,r2100_45SalinityMean,
                     r2100_45SalinityMin,r2100_45SalinityRange,r2100_45TemperatureMax,r2100_45TemperatureMean,
                     r2100_45TemperatureMin,r2100_45TemperatureRange,rjBathymetry)
names(env_2100_45)
names(env_2100_45) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2100_45)


#2100_60 ENVIRONMENT
r2100_60CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2100R60.tif")
r2100_60CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2100R60.tif")
r2100_60CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2100R60.tif")
r2100_60CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2100R60.tif")
r2100_60SalinityMax <- raster("reproj_SalinityMax2100R60.tif")
r2100_60SalinityMean <- raster("reproj_SalinityMean2100R60.tif")
r2100_60SalinityMin <- raster("reproj_SalinityMin2100R60.tif")
r2100_60SalinityRange <- raster("reproj_SalinityRange2100R60.tif")
r2100_60TemperatureMax <- raster("reproj_TemperatureMax2100R60.tif")
r2100_60TemperatureMean <- raster("reproj_TemperatureMean2100R60.tif")
r2100_60TemperatureMin <- raster("reproj_TemperatureMin2100R60.tif")
r2100_60TemperatureRange <- raster("reproj_TemperatureRange2100R60.tif")

projection(r2100_60CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_60TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2100_60 <- stack(r2100_60CurrentVelocityMax,r2100_60CurrentVelocityMean,r2100_60CurrentVelocityMin,
                     r2100_60CurrentVelocityRange,r2100_60SalinityMax,r2100_60SalinityMean,
                     r2100_60SalinityMin,r2100_60SalinityRange,r2100_60TemperatureMax,r2100_60TemperatureMean,
                     r2100_60TemperatureMin,r2100_60TemperatureRange,rjBathymetry)
names(env_2100_60)
names(env_2100_60) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2100_60)


#2100_85 ENVIRONMENT
r2100_85CurrentVelocityMax <- raster("reproj_CurrentVelocityMax2100R85.tif")
r2100_85CurrentVelocityMean <- raster("reproj_CurrentVelocityMean2100R85.tif")
r2100_85CurrentVelocityMin <- raster("reproj_CurrentVelocityMin2100R85.tif")
r2100_85CurrentVelocityRange <- raster("reproj_CurrentVelocityRange2100R85.tif")
r2100_85SalinityMax <- raster("reproj_SalinityMax2100R85.tif")
r2100_85SalinityMean <- raster("reproj_SalinityMean2100R85.tif")
r2100_85SalinityMin <- raster("reproj_SalinityMin2100R85.tif")
r2100_85SalinityRange <- raster("reproj_SalinityRange2100R85.tif")
r2100_85TemperatureMax <- raster("reproj_TemperatureMax2100R85.tif")
r2100_85TemperatureMean <- raster("reproj_TemperatureMean2100R85.tif")
r2100_85TemperatureMin <- raster("reproj_TemperatureMin2100R85.tif")
r2100_85TemperatureRange <- raster("reproj_TemperatureRange2100R85.tif")

projection(r2100_85CurrentVelocityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85CurrentVelocityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85CurrentVelocityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85CurrentVelocityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85SalinityMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85SalinityMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85SalinityMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85SalinityRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85TemperatureMax) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85TemperatureMean) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85TemperatureMin) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"
projection(r2100_85TemperatureRange) <- "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs"


env_2100_85 <- stack(r2100_85CurrentVelocityMax,r2100_85CurrentVelocityMean,r2100_85CurrentVelocityMin,
                     r2100_85CurrentVelocityRange,r2100_85SalinityMax,r2100_85SalinityMean,
                     r2100_85SalinityMin,r2100_85SalinityRange,r2100_85TemperatureMax,r2100_85TemperatureMean,
                     r2100_85TemperatureMin,r2100_85TemperatureRange,rjBathymetry)
names(env_2100_85)
names(env_2100_85) <- c("CVMax", "CVMean", "CVMin", "CVRange", "SalMax", "SalMean", "SalMin", "SalRange",
                        "TempMax", "TempMean","TempMin", "TempRange", "Bathy")
names(env_2100_85)

library(megaSDM)

#Setup a destination folder for the environmental outputs
envoutput <- "EnvironmentalOutput"

#Create the Training / Study environment
TSEnv <- TrainStudyEnv(input_TA = env_present,
                       output = envoutput,
                       clipTrain = c(-180,180,-90,90),
                       clipStudy = c(-180,180,-90,90),
                       desiredCRS = "+proj=cea +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +no_defs",
                       resolution = 5000)

#PredictEnv clips, resamples, and reprojects the future/past environmental layers to match the study region rasters created by the TrainStudyEnv function.

#PredictEnv requires rasters to be organized according to RCP scenario. 
#In the next section, we need to stack the rasterStacks of both 2050 and 2100 according to RCP

#«PredictEnv() function requires the list of time periods to include the time period you're training the data on (in this case, whatever time period the TSEnv$study data is in), in addition to the future data you would like to project the distributions onto. »

#On to PredictEnv, run once per RCP scenario
PredictEnv(studylayers = TSEnv$study,
           futurelayers = list(env_2050_26, env_2100_26),
           time_periods = c("present","2050","2100"),
           output = envoutput,
           scenario_name = "RCP2.6")

PredictEnv(studylayers = TSEnv$study,
           futurelayers = list(env_2050_45, env_2100_45),
           time_periods = c("present","2050","2100"),
           output = envoutput,
           scenario_name = "RCP4.5")

PredictEnv(studylayers = TSEnv$study,
           futurelayers = list(env_2050_60, env_2100_60),
           time_periods = c("present","2050","2100"),
           output = envoutput,
           scenario_name = "RCP6.0")

PredictEnv(studylayers = TSEnv$study,
           futurelayers = list(env_2050_85, env_2100_85),
           time_periods = c("present","2050","2100"),
           output = envoutput,
           scenario_name = "RCP8.5")


####Occurrence data import 

#Occurrence management

occlist <- list.files(pattern = ".csv", full.names = TRUE)
#View(occlist)

occ_output <- "occurrences"

OccurrenceManagement(occlist = occlist,
                     output = occ_output,
                     envextract = TRUE,
                     envsample = TRUE,
                     nbins = 25,
                     envdata = TSEnv$training)

#Background point generation and management

# Get the list of occurrence files again, even if they were written out 
# in the same folder as before. This ensures that the occurrence files 
# are properly formatted.  
occlist <- list.files(occ_output, pattern = ".csv", full.names = TRUE)

# The location to print out the background buffers (.shp) (will be created if it doesn't exist)
buff_output <- "buffers"

# Generates buffers for each species.
BackgroundBuffers(occlist = occlist,
                  envdata = TSEnv$training,
                  buff_output,
                  ncores = 8)

# Set the parameters for the background point generation 
# (how many points, and how spatially-constrained)

# How many background points should be generated per species?
nbg <- 1000

# What proportion of the background points should be sampled from within the buffers?
spatial_weights <- 0.5 

# Should the background points be environmentally subsampled (Varela) or 
# randomly distributed (random)?
sampleMethod <- "Varela" 

# Because we want a partial spatial constraint (50% of points within the buffer), we must make a 
# list of the buffer files to use in the creation of the background points. In the example, 
# these files are created from the BackgroundBuffers function, but they can also be generated
# outside of megaSDM and brought in here.

bufflist <- list.files(buff_output, pattern = ".shp$", full.names = TRUE)

# Define the location where the background points will be printed out to (.shp) 
# (This directory will be created if it doesn't already exist)
bg_output <- "backgrounds"

spplist <- c("Gymnodinium_catenatum",
             "Alexandrium_catenella",
             "Alexandrium_minutum")

BackgroundPoints(spplist = spplist,
                 envdata = TSEnv$training,
                 output = bg_output,
                 nbg = nbg,
                 spatial_weights = spatial_weights,
                 buffers = bufflist,
                 method = sampleMethod,
                 ncores = 8)

bglist <- list.files(bg_output, pattern = ".csv", full.names = TRUE)

#Distribution Modelling and Projection

# Define where the results of the MaxEnt model runs will be printed out to (as .lambdas files)
model_output <- "models"

# "nrep" is set to 4, meaning that the MaxEnt algorithm will run 4 times with different
# subsets of occurrence points for a better representation of the habitat suitability.
MaxEntModel(occlist = occlist,
            bglist = bglist,
            model_output = model_output,
            ncores = 8,
            nrep = 5,
            alloutputs = TRUE)

#MAXENT PROJECTION

# First, create a list of the time periods and climate scenarios used in the analysis
# (starting with the year the model is trained on)
time_periods <- c("present","2050","2100")
scenarios <- c("RCP2.6","RCP4.5","RCP6.0","RCP8.5")
# scenarios <- c("RCP2.6", "RCP4.5","RCP6.0","RCP8.5")
# Define the directory where the current study area rasters are located 
#    (generated from the TrainStudyEnv function or brought in from a separate location)
study_dir <- "EnvironmentalOutput/studyarea"

# Define the directories where the future study area rasters are located 
#    (generated from the PredictEnv function or brought in from a separate location)

# Define a list of directories for the projected climate layers, 
# separated into the different climate scenarios and years:
#    list(c(Scenario1Year1, Scenario1Year2), 
#         c(Scenario2Year1, Scenario2Year2))

predictdir <- list(c("EnvironmentalOutput/RCP2.6/2050",
                     "EnvironmentalOutput/RCP2.6/2100"),
                   c("EnvironmentalOutput/RCP4.5/2050",
                     "EnvironmentalOutput/RCP4.5/2100"),
                   c("EnvironmentalOutput/RCP6.0/2050",
                     "EnvironmentalOutput/RCP6.0/2100"),
                   c("EnvironmentalOutput/RCP8.5/2050",
                     "EnvironmentalOutput/RCP8.5/2100"))

# Define Where the results will be printed out.
# For this example, We'll define a new folder within
# the working directory that is specifically for the model
# projecions and analysis.
library(megaSDM)
result_dir <- "Results"

# Other options are also available (check the documentation page)
MaxEntProj(input = model_output,
           time_periods = time_periods,
           scenarios = scenarios,
           study_dir = study_dir,
           predict_dirs = predictdir,
           output = result_dir,
           aucval = 0.7,
           ncores = 8)

#ANALYSIS SECTION

#TIME MAPS
# The createTimeMaps function uses the binary maps at each time period to create "time maps":These maps detail the step -
#   wise expansions and contractions of the species distribution through those time -
#   steps, allowing for the visualization of both unidirectional range shifts (e.g., a range expansion across all time -
#  steps) and more complex dynamics (e.g.,a range expansion from 2010 - 2050 followed by a range contraction from 2050 
#                                                                                                                  

createTimeMaps(result_dir = result_dir,
               time_periods = time_periods,
               scenarios = scenarios,
               dispersal = FALSE,
               ncores = 8)

#ADDITIONAL STATISTICS

# The additionalStats function generates statistics on species range sizes and changes 
# through the multiple time steps and different scenarios. It also creates several graphs 
# (written out to the directory provided by "result_dir") to visualize these changes.

additionalStats(result_dir = result_dir,
                time_periods = time_periods,
                scenarios = scenarios,
                dispersal = FALSE,
                ncores = 8)





