     ##################################################################################
    ##################    Practical session of Species Distribution   ######################
    ##################      Summer school group4                  #####################
    ##################################################################################

# Install and load important packages

library(dismo)

library(raster)

library(rJava)

library(jsonlite)
##Import occurrence data 

occurence <- read.table("occurence.txt", fill=T, sep="\t", header=T, stringsAsFactors=F, quote="")

dim(occurence)

colnames(occurence)

occur1 <- subset(occurence, !is.na(decimalLongitude) & !is.na(decimalLatitude))



# Handling duplicates

dups <- duplicated(occur1[,c("decimalLongitude","decimalLatitude")])

sum(dups)

occur2 <-occur1[!dups,]

# Extract longitude and latitude data

occur3 <- occur2[, c(99,98)]

head(occur3,5)

names(occur3) <- c("lon","lat")



# Visualise your data on a map

library(geodata)

wrld <- world(path=".")

plot(wrld, xlim=c(-12,4), ylim=c(10,25), axes=TRUE,col="white")

points(occur3$lon, occur3$lat, col="red",pch=20, cex=0.75)

################################################################
library(raster)
getwd()
bio1 <- raster("C:/Users/HP/Desktop/Species Distribution/bio1_wc30s.tif")
bio2 <- raster("C:/Users/HP/Desktop/Species Distribution/bio2_wc30s.tif")
bio3 <- raster("C:/Users/HP/Desktop/Species Distribution/bio3_wc30s.tif")
bio4 <- raster("C:/Users/HP/Desktop/Species Distribution/bio4_wc30s.tif")
bio5 <- raster("C:/Users/HP/Desktop/Species Distribution/bio5_wc30s.tif")
bio6<- raster("C:/Users/HP/Desktop/Species Distribution/bio6_wc30s.tif")
bio7 <- raster("C:/Users/HP/Desktop/Species Distribution/bio7_wc30s.tif")
bio10 <- raster("C:/Users/HP/Desktop/Species Distribution/bio10_wc30s.tif")
bio11 <- raster("C:/Users/HP/Desktop/Species Distribution/bio11_wc30s.tif")
bio12 <- raster("C:/Users/HP/Desktop/Species Distribution/bio12_wc30s.tif")
bio13 <- raster("C:/Users/HP/Desktop/Species Distribution/bio13_wc30s.tif")
bio14 <- raster("C:/Users/HP/Desktop/Species Distribution/bio14_wc30s.tif")
bio15 <- raster("C:/Users/HP/Desktop/Species Distribution/bio15_wc30s.tif")
bio16 <- raster("C:/Users/HP/Desktop/Species Distribution/bio16_wc30s.tif")
bio17 <- raster("C:/Users/HP/Desktop/Species Distribution/bio17_wc30s.tif")
pet <- raster("C:/Users/HP/Desktop/Species Distribution/pet_wc30s.tif")
mimq <- raster("C:/Users/HP/Desktop/Species Distribution/mimq_wc30s.tif")
llds <- raster("C:/Users/HP/Desktop/Species Distribution/llds_wc30s.tif")
miaq <- raster("C:/Users/HP/Desktop/Species Distribution/miaq_wc30s.tif")
mi <- raster("C:/Users/HP/Desktop/Species Distribution/mi_wc30s.tif")
dm<- raster("C:/Users/HP/Desktop/Species Distribution/dm_wc30s.tif")

################### COrrelation analysis  #############################

#we have made a correlation analysis to see the correlation between 21 bioclimatic 
#variables.
library(ENMTools)
corr <- stack(bio1, bio2, bio3, bio4, bio5, bio6, bio7, bio10, bio11, bio12, bio13,
              bio14, bio15, bio16, bio17, pet, mimq, llds, miaq, mi, dm)
env<- extract(corr, cbind(occur3$lon, occur3$lat))
env <-as.data.frame(env)
cor_matrix <- cor(env, method = "pearson")
print(cor_matrix)
#######################  Modelling with chosen variables ###############################

# We have selected e.g 10 bioclimatic variables for the initial model based on their correlation 
#and ecological  consideration.

#initial <- c(llds, bio2, bio3, bio4, bio10, bio11, bio12, bio17, pet, mimq)

#In the final model we have used 5 variables selected based on their contribution 
#in the initial model

########## Modelling with the chosen variables #######################################
historical <- stack(bio12,bio4,mimq, bio11, bio2)
plot(historical)
env_occur <- extract(historical, cbind(Bombax3$lon, Bombax3$lat))
env_occur <-as.data.frame(env_occur)
env_occur <- cbind.data.frame(Bombax3$lon, Bombax3$lat, env_occur)
env_occur1 <- as.data.frame(na.omit(env_occur))
names(env_occur1)
# Rename the columns
names(env_occur1) <- c("lon","lat","bio12","bio4","mimq","bio11", "bio2")
###The map of Mali
library(sf)
Ml <- sf::st_read("country_adm0.shp")
plot(Ml)
plot(Ml$geometry)
##To get the inform about env condition on country based on the map of africa#
historical1 <- crop(historical, extent(Ml))
historical1 <- mask(historical1, Ml)
names(historical1)
plot(historical1)
########
names(historical1) <- c("bio12","bio4","mimq", "bio11", "bio2")
#View(historical1)
# # Getting all data from the rasters of study area
myxy <- xyFromCell(historical1[[1]],cell=1:ncell(historical1[[1]]))
myvar <- getValues(historical1)
study_area <-cbind.data.frame(myxy, myvar)

# Delete na from the database
study_area <- as.data.frame(na.omit(study_area))
names(study_area) <- c("lon","lat", "bio12","bio4","mimq", "bio11", "bio2")

# generate a matrix containing 1 to depict presence and 0 to depict background (the study area)
pa <-c(rep(1,nrow(env_occur1)),rep(0,nrow(study_area)))

model_data <- rbind(env_occur1[,3:7],study_area[,3:7]) # Note that we excluded the x and y information.

# Specify where you want to store the results
##setwd("C:/Users/HP/Desktop/MBIOST_Benin/M1/Cours/Master_2/Species Distribution")
wd <- getwd()

# Call packages that are necessary to run MaxEnt in R (rJava)
options(java.parameters = "-Xmx4g")

library(rJava)
library(dismo)
# Model
MSWD <- maxent(x=model_data, p=pa, path=paste(wd, "/thmaxSWD", sep=""), args=c("-P", "noautofeature", "nothreshold", "noproduct", "maximumbackground=1534459", "noaddsamplestobackground"))

# visualization
windows()
plot(MSWD)
response(MSWD)

MSWD.pred1 <- predict(MSWD, historical1, args="outputformat=raw", sep="", progress="text")                                  
windows()
plot(MSWD.pred1)
getwd()
## Projecting to the future bioclimatic environment

# importing future bioclimate data
bio2f<- raster("C:/Users/HP/Desktop/Species Distribution/bio2_rcp45_2055_wc30s.tif")
bio12f <- raster("C:/Users/HP/Desktop/Species Distribution/bio12_rcp45_2055_wc30s.tif")
bio4f <- raster("C:/Users/HP/Desktop/Species Distribution/bio4_rcp45_2055_wc30s.tif")
mimqf <- raster("C:/Users/HP/Desktop/Species Distribution/mimq_rcp45_2055_wc30s.tif")
bio11f <- raster("C:/Users/HP/Desktop/Species Distribution/bio11_rcp45_2055_wc30s.tif")

#plot(lldsf)

future <- stack(bio2f, bio12f,bio4f,mimqf,bio11f)
plot(future)

# Crop the African maps of the future bioclimate to Benin
library(sf)
Ml <- sf::st_read("MLI_adm0.shp")
plot(Ml)
plot(Ml$geometry)

future1 <- crop(future,extent(Ml))
plot(future1)
future1<- mask(future1,Ml)
plot(future1)
names(future1)
names(future1) <- c("bio2","bio12","bio4","mimq","bio11")

# predicting the future
MSWD.pred2 <- predict(MSWD, future1, args="outputformat=raw", sep="", progress="text")
plot(MSWD.pred2)

# visualiseles the modeled habitat suitability maps

MSWD.pred2 <- crop(MSWD.pred2,extent(MSWD.pred1))
MSWD.pred2 <- resample(MSWD.pred2,MSWD.pred1, method="bilinear")
MSWD.pred2 <- mask(MSWD.pred2,MSWD.pred1)
results <- stack(MSWD.pred1,MSWD.pred2)
windows()
plot(results)

names(results) <- c("mdesp historical","mdesp future")
windows()
plot(results)

# Save the results for further handling in ArcMap
#writeRaster(project, filename = "projection.tif")

writeRaster(MSWD.pred1,filename="Bombax.tif")
writeRaster(MSWD.pred2,filename="Bombaxf.tif")
writeRaster(results,filename="resultsf.tif")

##END
