##installation of packages
library(raster)
library(pacman)

### Selection of the folder where the files are located
setwd("D:/Artículos/PAPER CEDRELA/GRAFICOS/areas generadas")
mps <- shapefile("F:/Capas ambientales/Presente (1980-2010)/recortadas_cedrela/shape corte1.shp")
#mps
plot(mps)
####load all raster layers
Bio1<-raster("BIO1.tif")
Bio2<-raster("BIO2.tif")
Bio3<-raster("BIO3.tif")
Bio4<-raster("BIO4.tif")
Bio5<-raster("BIO5.tif")
Bio6<-raster("BIO6.tif")
Bio7<-raster("BIO7.tif")
Bio8<-raster("BIO8.tif")
Bio9<-raster("BIO9.tif")
Bio10<-raster("BIO10.tif")
Bio11<-raster("BIO11.tif")
Bio12<-raster("BIO12.tif")
Bio13<-raster("BIO13.tif")
Bio14<-raster("BIO14.tif")
Bio15<-raster("BIO15.tif")
Bio16<-raster("BIO16.tif")
Bio17<-raster("BIO17.tif")
Bio18<-raster("BIO18.tif")
Bio19<-raster("BIO19.tif")

alti<-raster("alt.tif")
ai_<-raster("ai_et.tif")
et_<-raster("et_yr.tif")
def<-raster("def_total.tif")
defl<-raster("loss_total.tif")
defg<-raster("gain_total.tif")
BD<-raster("BD.tif")
CC<-raster("CC.tif")
CEC<-raster("CEC.tif")
CF<-raster("CF.tif")
NIT<-raster("nitrogen.tif")
OCD<-raster("OCD.tif")
PH<-raster("PH.tif")
Sand<-raster("sand.tif")
Slit<-raster("slit.tif")
SOC<-raster("SOC.tif")
SOCS<-raster("SOCS.tif")
PET<-raster("PET.tif")
CMI<-raster("CMI.tif")
FGD<-raster("FGD.tif")
GSL<-raster("GSL.tif")
GSP<-raster("GSP.tif")
GST<-raster("GST.tif")
HURS<-raster("HURS.tif")
NDVI<-raster("NDVI-2018.tif")

####cut back on the study area
biop=crop(BIOp,mps)
biop=mask(biop,mps)

bio1=crop(Bio1, mps)
bio1=mask(bio1, mps)

bio2=crop(Bio2, mps)
bio2=mask(bio2, mps)

bio3=crop(Bio3, mps)
bio3=mask(bio3, mps)

bio4=crop(Bio4, mps)
bio4=mask(bio4, mps)

bio5=crop(Bio5, mps)
bio5=mask(bio5, mps)

bio6=crop(Bio6, mps)
bio6=mask(bio6, mps)

bio7=crop(Bio7, mps)
bio7=mask(bio7, mps)

bio8=crop(Bio8, mps)
bio8=mask(bio8, mps)

bio9=crop(Bio9, mps)
bio9=mask(bio9, mps)

bio10=crop(Bio10, mps)
bio10=mask(bio10, mps)

bio11=crop(Bio11, mps)
bio11=mask(bio11, mps)

bio12=crop(Bio12, mps)
bio12=mask(bio12, mps)

bio13=crop(Bio13, mps)
bio13=mask(bio13, mps)

bio14=crop(Bio14, mps)
bio14=mask(bio14, mps)

bio15=crop(Bio15, mps)
bio15=mask(bio15, mps)

bio16=crop(Bio16, mps)
bio16=mask(bio16, mps)

bio17=crop(Bio17, mps)
bio17=mask(bio17, mps)

bio18=crop(Bio18, mps)
bio18=mask(bio18, mps)

bio19=crop(Bio19, mps)
bio19=mask(bio19, mps)

Alti=crop(alti, mps)
Alti=mask(Alti, mps)

ari=crop(ai_, mps)
ari=mask(ari, mps)

evt=crop(et_, mps)
evt=mask(evt, mps)

Def=crop(def, mps)
Def=mask(Def, mps)

bd=crop(BD, mps)
bd=mask(bd, mps)

cc=crop(CC, mps)
cc=mask(cc, mps)

cec=crop(CEC, mps)
cec=mask(cec, mps)

cf=crop(CF, mps)
cf=mask(cf, mps)

nit=crop(NIT, mps)
nit=mask(nit, mps)

ocd=crop(OCD, mps)
ocd=mask(ocd, mps)

ph=crop(PH, mps)
ph=mask(ph, mps)

sand=crop(Sand, mps)
sand=mask(sand, mps)

slit=crop(Slit, mps)
slit=mask(slit, mps)

soc=crop(SOC, mps)
soc=mask(soc, mps)

socs=crop(SOCS, mps)
socs=mask(socs, mps)

pet=crop(PET,mps)
pet=mask(pet, mps)

cmi=crop(CMI, mps)
cmi=mask(cmi, mps)

fgd=crop(FGD, mps)
fgd=mask(fgd, mps)

gsl=crop(GSL, mps)
gsl=mask(gsl, mps)

gsp=crop(GSP, mps)
gsp=mask(gsp, mps)

gst=crop(GST, mps)
gst=mask(gst, mps)

hurs=crop(HURS, mps)
hurs=mask(hurs, mps)

ndvi=crop(NDVI, mps)
ndvi=mask(ndvi, mps)

deflo=crop(defl, mps)
deflo=mask(deflo, mps)

defga=crop(defg, mps)
defga=mask(defga, mps)

####image resolution verification
res(biop)
res(bio1)
res(bio2)
res(bio3)
res(ari)
res(evt)
res(def)
res(NIT)
res(hurs)

#############convert all layers to the same pixel resolution
rbio1 <- resample(bio1,bio1,resample='bilinear')
rbio2 <- resample(bio2,bio1,resample='bilinear')
rbio3 <- resample(bio3,bio1,resample='bilinear')
rbio4 <- resample(bio4,bio1,resample='bilinear')
rbio5 <- resample(bio5,bio1,resample='bilinear')
rbio6 <- resample(bio6,bio1,resample='bilinear')
rbio7 <- resample(bio7,bio1,resample='bilinear')
rbio8 <- resample(bio8,bio1,resample='bilinear')
rbio9 <- resample(bio9,bio1,resample='bilinear')
rbio10 <- resample(bio10,bio1,resample='bilinear')
rbio11 <- resample(bio11,bio1,resample='bilinear')
rbio12 <- resample(bio12,bio1,resample='bilinear')
rbio13 <- resample(bio13,bio1,resample='bilinear')
rbio14 <- resample(bio14,bio1,resample='bilinear')
rbio15 <- resample(bio15,bio1,resample='bilinear')
rbio16 <- resample(bio16,bio1,resample='bilinear')
rbio17 <- resample(bio17,bio1,resample='bilinear')
rbio18 <- resample(bio18,bio1,resample='bilinear')
rbio19 <- resample(bio19,bio1,resample='bilinear')
rdef <- resample(Def,bio1,resample='bilinear')
rbd <- resample(bd,bio1,resample='bilinear')
rcc <- resample(cc,bio1,resample='bilinear')
rcec <- resample(cec,bio1,resample='bilinear')
rcf <- resample(cf,bio1, resample='bilinear')
rnit <- resample(nit,bio1,resample='ngb')
rocd <- resample(ocd,bio1,resample='bilinear')
rph <- resample(ph,bio1,resample='bilinear')
rsand <- resample(sand,bio1,resample='bilinear')
rslit <- resample(slit,bio1,resample='bilinear')
rsoc <- resample(soc,bio1,resample='bilinear')
rsocs <- resample(socs,bio1,resample='bilinear')
rari <- resample(ari,bio1,resample='bilinear')
revt <- resample(evt,bio1,resample='bilinear')
ralt<-resample(Alti,bio1,resample='bilinear')
rpet<-resample(pet,bio1,resample='bilinear')
rcmi<-resample(cmi,bio1,resample='bilinear')
rfgd<-resample(fgd,bio1,resample='bilinear')
rgsl<-resample(gsl,bio1,resample='bilinear')
rgsp<-resample(gsp,bio1,resample='bilinear')
rgst<-resample(gst,bio1,resample='bilinear')
rhurs<-resample(hurs,bio1,resample='bilinear')
rndvi<-resample(ndvi,bio1,resample='bilinear')
rdefl<-resample(deflo,bio1,resample='bilinear')
rdefg<-resample(defga,bio1,resample='bilinear') 

junt1 <- stack(rbio1,rbio2,rbio3,rbio4,rbio5,rbio6,rbio7,rbio8,rbio9,rbio10,rbio11,
                rbio12,rbio13,rbio14,rbio15,rbio16,rbio17,rbio18,rbio19,rdef,rbd,rcc,
                rcec,rcf,rnit,rocd,rph,rsand,rslit,rsoc,rsocs,ralt,rpet,
                rcmi,rhurs,rndvi,rdefl,rdefg)

#plot(junt1)
writeRaster(junt1, filename = "recort_present2.grd", overwrite = T)

#####PREDICTIONN##############
library(dismo)
library(rJava)
library(raster)
library(caTools)
library(sp)
library(terra)

#--------------
predictors <- read.csv("Cedrela_angus.csv", ";", header = TRUE)
predictors
samp_train <- predictors[,(3:4)]
head(samp_train)
plot(samp_train)

layers <- stack("recort_present2.grd")
plot(layers,1)
points(samp_train$lon, samp_train$lat, col="red", pch=20, cex=0.75)
class(layers)

#-----Variables selection------
presvals <-raster::extract(layers, samp_train)
presvals
write.csv(presvals, file ="present.csv")
head(presvals)
summary(presvals)
hist(presvals[,20], main="Bioclimatic Profile BIO1",xlab="MAT",
     ylab="Presence")
dat<-presvals
names(dat)
head(dat)

library(adehabitatHR)
library(virtualspecies)
library(corrplot)
cor_specie<-cor(dato, method = "pearson")
cor_specie
write.csv(cor_specie, "Correlacion bio cedrela.csv")

corrplot(cor(cor_specie))
corrplot(cor(cor_specie,method = "pearson"), tl.col = "black", diag=F)
corrplot.mixed(cor(cor_specie),number.cex = .8, tl.col ="black")

#-----Selection of environmental variables r>0.8-----
predictors <- stack(layers)
predictorss <- removeCollinearity(predictors,multicollinearity.cutoff = 0.85,
                                  select.variables =TRUE, sample.points = FALSE)
predictorss

predictors1<-layers
predictors2<-stack(layers$BIO5, layers$BIO2, layers$BIO3, layers$BIO7,layers$BIO12,
                   layers$BIO15,layers$BIO19,layers$BD,layers$CC, layers$CEC,layers$CF,
                   layers$nitrogen,layers$OCD,layers$PH,layers$sand, layers$slit,layers$SOC,
                   layers$SOCS,layers$PET,layers$HURS,layers$NDVI.2018)#seleccion de variables que entrar? al modelo

predictors3<-stack( layers$BIO3, layers$BIO4,layers$BIO6, layers$BIO10, 
                    layers$BIO15, layers$BIO18, layers$BIO19, layers$CMI,
                    layers$alt,layers$SOCS,layers$OCD,layers$slit,layers$CC,
                    layers$NDVI.2018)
#layers$def_total,layers$loss_total

#-----------Entrenamiento de modelos--------
library(dismo)
pred_nf <- dropLayer(predictors3, 'biome')
#pred_nf
set.seed(7)
group <- kfold(samp_train, 5)
pres_train <- samp_train[group !=1, ]
pres_train
pres_test <- samp_train[group ==1, ]
pres_test
#obs<-rbind(pres_train,pres_test)
#obs
#-----------Validation set--------------
set.seed(1)
backg <- randomPoints(pred_nf, n=1000)
colnames(backg) = c('lon','lat')
group <- kfold(backg, 5)
backg_train <- backg[group !=1, ]
backg_train
write.csv(backg_train,"bg_train.csv")
backg_test <- backg[group ==1, ]
backg_test
obs<-rbind(backg_train,backg_test)
write.csv(obs,"obs_eval.csv")

#######RANDOM FOREST#######################
library(rgeos)
library(knitr)
library(randomForest)

pres_train <- as.data.frame(pres_train[,(1:2)])
pres_test
colnames(pres_train) <- colnames(backg_train[,(1:2)])
train <- rbind(pres_train,backg_train)
test<-rbind(pres_test,backg_test)
pb_train <- c(rep(1, nrow(pres_train)), rep(0,nrow(backg_train[,(1:2)])))
envtrain <- raster::extract(predictors3, train)
envtrain <- data.frame(cbind(pa=pb_train, envtrain))
head(envtrain)

testpres <- data.frame(raster::extract(predictors3,pres_test))
testbackg <- data.frame(raster::extract(predictors3, backg_test))


presvals <- raster::extract(predictors3, samp_train, df = TRUE)
set.seed(0)
backgr <- randomPoints(predictors3, nrow(presvals))

absvals<- raster::extract(predictors3, backgr)
pb <- c(rep(1,nrow(presvals)), rep(0,nrow(absvals)))
sdmdata <- data.frame(cbind(pb, rbind(presvals[-1], absvals)))
sdmdata[,'pb']=as.factor(sdmdata[,'pb'])
head(sdmdata)
tail(sdmdata)

#-------Saving of training data-----------------
delete.na <- function(DF, n=0) { DF[rowSums(is.na(DF)) <=n,]}
sdmdata <- delete.na(sdmdata)

#----modeling----
samp <- sample(nrow(sdmdata), round(0.75 * nrow(sdmdata)))
traindata <-sdmdata[samp,]
testdata <- sdmdata[-samp,]

envtrain <- delete.na(envtrain)
model <- as.factor(pa) ~.
rf1 <- randomForest(pa~., data=envtrain)
rf1

predictions <- as.data.frame(predict(rf1, test, type = "prob"))

plot(rf1,main = "Error decrease RF model")
erf <- evaluate(traindata, testdata, rf1)
erf
varImpPlot(rf1,type=2,main = "variable importance - Random Forest",
           pch=19)
plot(erf, 'ROC')
pr <- predict(predictors3, rf1) #prediction map
pr

plot(pr, main='Random Forest')
tr3 <- threshold(erf, 'spec_sens')
tr3
irf<-importance(rf1)
irf
write.csv(irf, "imp_future100_370_cd.csv")

varImpPlot(rf1, sort = TRUE)

pd_rf_test <- predict(rf1, presvals)
pd_rf_test 
plot(pr >0.25, main ='presence/absence')
points(pres_train, pch='+')
points(backg_train, pch='-', cex= 0.25)
writeRaster(pr, "cedrela_future100-585_sd", format="GTiff",overwrite=TRUE)

#####MAXENT######################
library(dismo)
library(rJava)
fold <- kfold(samp_train, k=5)
occtest <- samp_train[fold == 1, ]
occtrain <- samp_train[fold != 1, ]

me1 <- maxent(predictors3, samp_train, args=c(
  'maximumbackground=10000',
  'defaultprevalence=1.00',
  'betamultiplier=1.0',
  'pictures=true',
  'linear=true',
  'quadratic=true',
  'product=true',
  'threshold=false',
  'hinge=true',
  'threads=1',
  'responsecurves=true',
  'jackknife=true',
  'askoverwrite=true',
  'replicatetype=Crossvalidate',
  'replicates=10',
  'randomseed=TRUE',
  'plots=TRUE',
  'writeplotdata=TRUE'
))

me1 <- maxent(predictors3, samp_train, args=c(
  'maximumbackground=10000',
  'defaultprevalence=1.00',
  'betamultiplier=1.0',
  'pictures=true',
  'linear=true',
  'quadratic=true',
  'product=true',
  'threshold=false',
  'hinge=true',
  'threads=1',
  'responsecurves=true',
  'jackknife=true',
  'askoverwrite=true',
  'replicatetype=Crossvalidate',
  'replicates=10',
  'randomseed=TRUE',
  'plots=TRUE',
  'writeplotdata=TRUE'
))
# plot showing importance of each variable
rem<-me1@results
write.csv(rem, "maxent_results present sd.csv")
plot(me1@models)

# response curves
res1<-response(me1)

# predict to entire dataset
r <- predict(me1, predictors3) 
Summary(r)
plot(r)
points(occtrain)
writeRaster(r$layer.5, "maxent_present4.tif", overwrite = TRUE)

#testing
# background data
bg <- randomPoints(predictors3, 1000)

#simplest way to use 'evaluate'
e1 <- evaluate(me1, p=occtest, a=bg, x=predictors3)
e1

# alternative 1
# extract values
pvtest <- data.frame(extract(predictors3, occtest))
avtest <- data.frame(extract(predictors3, bg))

e2 <- evaluate(me1, p=pvtest, a=avtest)
e2
# alternative 2 
# predict to testing points 
testp <- predict(me1, pvtest) 
head(testp)
testa <- predict(me1, avtest) 

e3 <- evaluate(p=testp, a=testa)
e3
threshold(e3)

plot(e3, 'ROC')