
                                                ###################################################
                                                #                                                 #
                                                #      READING IN AND CLEANING OCCURRENCE DATA    #
                                                #                                                 # 
                                                ###################################################


# Importing data, view and edit data file as needed
SABAP2raw <- read.delim("FILEPATH")
recs_SABAP2 <- SABAP2raw[,c( "decimalLongitude","decimalLatitude","species", "countryCode")]

# Cleaning and cross-checking
library(CoordinateCleaner)
subset(recs_SABAP2, !is.na(decimalLatitude))
cl_recsSABAP2 <- clean_coordinates(recs_SABAP2, lon="decimalLongitude",vlat="decimalLatitude", countries="countryCode",tests=c("centroids","outliers"))
recs_SABAP2 <- recs_SABAP2[cl_recsSABAP2$.summary,]
head(recs_SABAP2)



                                               #####################################################
                                               #                                                   #
                                               #          VIEWING/MAPPING OCCURENCE DATA           #
                                               #                                                   #
                                               #####################################################


#View these data on a map
library(maptools)
data(wrld_simpl) #Get the world map
sa = wrld_simpl[wrld_simpl$ISO2%in%c('ZA','NA','BW','ZW','MZ','LS','SZ'),] 
plot(sa) #Plot southern African countries
points(recs_SABAP2$decimalLongitude, recs_SABAP2$decimalLatitude, col='red')



                                              #######################################################
                                              #                                                     #
                                              #    IMPORTING BIOCLIMATIC DATA/PREDICTOR VARIABLES   # 
                                              #                                                     #
                                              #######################################################


# Downloading data
library(raster)
ext <- extent( c(0,60,-40,0))    #defining extent
clim_curr <- getData("worldclim", var="bio", res=5, download=T)
clim_curr <- crop(clim_curr,ext)

# Testing for collinearity		
library(usdm)
cortest2 <-vifcor(clim_curr,0.7)
cortest2

# Dropping collinear predictors as identified by vifcor() 	
climSABAP2 <- exclude(clim_curr,cortest2) 
names(climSABAP2) #Final set of predictors



                                             ########################################################
                                             #                                                      #
                                             #       JOINING OCCURENCE AND (BIO)CLIMATIC DATA       #
                                             #                                                      #
                                             ########################################################


# Joining species and climate data
data_SABAP2 <- cbind(recs_SABAP2, extract(x = climSABAP2, y = data.frame(recs_SABAP2[,c('decimalLongitude','decimalLatitude')]))) 
final_SABAP2 <- data_SABAP2[,c("decimalLongitude","decimalLatitude", "species", "bio2", "bio8", "bio9", "bio15", "bio18", "bio19")]
final_SABAP2$species <-1
head(final_SABAP2)

# Making sure objects are spatial objects
class(final_SABAP2)
coordinates(final_SABAP2) <- ~decimalLongitude + decimalLatitude
class(final_SABAP2)

# Checking and fixing projection issues
projection(final_SABAP2)
proj4string(final_SABAP2) <- projection(climSABAP2)
projection(final_SABAP2)




                                            #########################################################
                                            #                                                       #
                                            #               MODEL FITTING AND ASSESSMENT            #
                                            #                                                       #
                                            #########################################################


library(sdm)
installAll() # use only if running sdm for the first time
  
SABAP2.sdm <- sdmData(species~., train = final_SABAP2, predictors = climSABAP2, bg = list(n=1000))

# Fitting models via four algorithms;
# getmethodNames() displays all available algorithms
fit_SABAP2 <- sdm(species~., SABAP2.sdm, methods=c("gam","svm","rf","mars"), replication=c("boot"), n=10)
fit_SABAP2
roc(fit_SABAP2) 



                                            #########################################################
                                            #                                                       #
                                            #              MAKING AND USING ENSEMBLES               #
                                            #                                                       #
                                            #########################################################
 

# Current species distribution
SABAP2curr <- ensemble(fit_SABAP2, climSABAP2, filename = 'ens.img', setting = list(method='weighted', stat='AUC')) #creating ensemble
plot(SABAP2curr)
points(final_SABAP2)

# Future species distribution
clim_fut <-raster::getData("CMIP5", var='bio',res=5,model='AC',rcp=85,year=70)
names(clim_fut) <- names(clim_curr)
clim_fut <- exclude(clim_fut,cortest2)
clim_fut <- crop(clim_fut,ext)

# Ensemble prediction
SABAP2fut <- ensemble(fit_SABAP2, clim_fut, 'ensf.img', setting = list(method='weighted', stat='AUC')) 
plot(SABAP2fut)
points(final_SABAP2)



                                           ##########################################################
                                           #                                                        #
                                           #                INFERENCE ON RANGE CHANGES              #
                                           #                                                        #
                                          ###########################################################        


# Compare current and future distributions
rangeChange <- SABAP2fut - SABAP2curr
clz <- colorRampPalette(c('red', 'orange', 'gray', 'green', 'blue'))
plot(rangeChange, col=clz(200))





