options(java.parameters = "-Xmx50000m")
library(rgdal)
library(raster)
library(dismo)
library(rJava)

### MODELLING Syagrus romanzoffiana DISTRIBUTION

#### Sciurus aestuans occurence points ####
pontos<- read.csv('syagrus_romanzoffiana.csv', h=T)

pts.geo<-SpatialPoints(pontos[,1:2])
pts.geo.df<-SpatialPointsDataFrame(pontos[,1:2], as.data.frame(pontos['Especie']))

##### Before beginning models ans stacks, we excluded those which were correlated
#### to avoid colinearity, using a threshold of 0.6.
#### After exclusion, we stayed with six:
#### climate and biotic predictors

bio2<-raster('variaveis_atuais/bio2_var.tif')
bio3<-raster('variaveis_atuais/bio3_var.tif')
bio8<-raster('variaveis_atuais/bio8_var.tif')
bio15<-raster('variaveis_atuais/bio15_var.tif')
bio18<-raster('variaveis_atuais/bio18_var.tif')
bio19<-raster('variaveis_atuais/bio19_var.tif')

#### making a stack with all variables

clima <- stack(bio2,bio3,bio8,bio15,bio18,bio19)
clima.df <- as.data.frame(clima, xy=T)
clima.df<- na.omit(clima.df)
gridded(clima.df) <- ~ x+y
clima.st<-stack(clima.df)


#### Making a raster with all variables

clima.mod<-clima.st
clim.pts<- extract(clima.mod, pts.geo, cellnumbers=T)
syagrus.var<- cbind(pontos[,1:2], clim.pts)

### excluding points in the same pixel

# duplicated(syagrus.var[,'cells'])
# dup<- which(duplicated(syagrus.var[,'cells']) == TRUE)
# syagrus.var <- syagrus.var [-dup,]
# syagrus.var.limpo <- na.omit(syagrus.var)
write.table(syagrus.var, 'syagrus_clim.txt')

##### making points of background/pseudo-ausence####

syagrus.oc<-read.table('syagrus_clim.txt', h=T)
write.table(syagrus.oc, 'syagrus_clim.txt')
syagrus.oc<-read.table('syagrus_clim.txt')

clima.mod.df<-na.omit(as.data.frame(clima.mod, xy=T))

back.id<-sample(1:nrow(clima.mod.df), nrow(syagrus.oc))

back<-clima.mod.df[back.id,]

write.table(back, 'background_syagrus.txt')
back<-read.table('background_syagrus.txt')

#### cross validation = 10 
#### each run is a simulation which will compose the ensemble
Sys.setenv(NOAWT=T)
library(rJava)

### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model1<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict1<-predict(maxent.model1,clima.mod)
writeRaster(maxent.predict1,'maxent_atual_syag1.tif')



### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model2<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict2<-predict(maxent.model2,clima.mod)
writeRaster(maxent.predict2,'maxent_atual_syag2.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model3<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict3<-predict(maxent.model3,clima.mod)
writeRaster(maxent.predict3,'maxent_atual_syag3.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model4<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict4<-predict(maxent.model4,clima.mod)
writeRaster(maxent.predict4,'maxent_atual_syag4.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model5<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict5<-predict(maxent.model5,clima.mod)
writeRaster(maxent.predict5,'maxent_atual_syag5.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model6<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict6<-predict(maxent.model6,clima.mod)
writeRaster(maxent.predict6,'maxent_atual_syag6.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model7<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict7<-predict(maxent.model7,clima.mod)
writeRaster(maxent.predict7,'maxent_atual_syag7.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model8<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict8<-predict(maxent.model8,clima.mod)
writeRaster(maxent.predict8,'maxent_atual_syag8.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model9<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict9<-predict(maxent.model9,clima.mod)
writeRaster(maxent.predict9,'maxent_atual_syag9.tif')


### Training and test datasets ###

sample.ocor <- sample(1:nrow(syagrus.oc), round(0.75*nrow(syagrus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=syagrus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=syagrus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model10<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.predict10<-predict(maxent.model10,clima.mod)
writeRaster(maxent.predict10,'maxent_atual_syag10.tif')


### Making ensemble of Syagrus romanzoffiana Distribution

ens.mean.syagrus<-sum(maxent.predict1,maxent.predict2,maxent.predict3,maxent.predict4,maxent.predict5,maxent.predict6,maxent.predict7,maxent.predict8,maxent.predict9,maxent.predict10)
ens.mean.syagrus<-(ens.mean.syagrus/10)

writeRaster(ens.mean.syagrus,'syagrus_var.tif')


### MODELLING WITH BIOTIC RESOURCE

#### Sciurus aestuans occurence points ####
pontos<- read.csv('Sciurus.csv', h=T)

pts.geo<-SpatialPoints(pontos[,1:2])
pts.geo.df<-SpatialPointsDataFrame(pontos[,1:2], as.data.frame(pontos['Especie']))

##### Before beginning models ans stacks, we excluded those which were correlated
#### to avoid colinearity, using a threshold of 0.6.
#### After exclusion, we stayed with six:
#### climate and biotic predictors

bio2<-raster('variaveis_atuais/bio2_var.tif')
bio3<-raster('variaveis_atuais/bio3_var.tif')
bio8<-raster('variaveis_atuais/bio8_var.tif')
bio15<-raster('variaveis_atuais/bio15_var.tif')
bio18<-raster('variaveis_atuais/bio18_var.tif')
bio19<-raster('variaveis_atuais/bio19_var.tif')
syagrus<-raster('variaveis_atuais/syagrus_var.tif')

#### making a stack with all variables

clima <- stack(syagrus,bio2,bio3,bio8,bio15,bio18,bio19)
clima.df <- as.data.frame(clima, xy=T)
clima.df<- na.omit(clima.df)
gridded(clima.df) <- ~ x+y
clima.st<-stack(clima.df)


#### Making a raster with all variables

clima.mod<-clima.st
clim.pts<- extract(clima.mod, pts.geo, cellnumbers=T)
sciurus.var<- cbind(pontos[,1:2], clim.pts)

### excluding points in the same pixel

duplicated(sciurus.var[,'cells'])
dup<- which(duplicated(sciurus.var[,'cells']) == TRUE)
dup
sciurus.var <- sciurus.var [-dup,]
dim(sciurus.var)
sciurus.var.limpo <- na.omit(sciurus.var)
write.table(sciurus.var.limpo, 'sciurus_clim.txt')

##### making points of background/pseudo-ausence####

sciurus.oc<-read.table('sciurus_clim.txt', h=T)
write.table(sciurus.oc, 'sciurus_clim_withrec.txt')
sciurus.oc<-read.table('sciurus_clim_withrec.txt')

clima.mod.df<-na.omit(as.data.frame(clima.mod, xy=T))

back.id<-sample(1:nrow(clima.mod.df), nrow(sciurus.oc))

back<-clima.mod.df[back.id,]

write.table(back, 'background_sciurus_comrec.txt')
back<-read.table('background_sciurus_comrec.txt')

#### cross validation = 10 
#### each run is a simulation which will compose the ensemble
Sys.setenv(NOAWT=T)
library(rJava)

sample(1:1000, 20, replace=FALSE)

### Training and test datasets ###
set.seed(953)

sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### modelling with maxent ####

maxent.model1<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model1@results
maxent.predict<-predict(maxent.model1,clima.mod)
writeRaster(maxent.predict,'maxent_atual_sci_rec_red1.tif')
maxent.eval<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model1)
maxent.eval
maxent.th<-threshold(maxent.eval)
maxent.th$spec_sens
maxent.th$kappa
maxent.th$prevalence
maxent.th$sensitivity
maxent.pred.bin<-maxent.predict
maxent.pred.bin[maxent.pred.bin>=maxent.th$spec_sens]=1
maxent.pred.bin[maxent.pred.bin<maxent.th$spec_sens]=0
writeRaster(maxent.pred.bin, 'maxent_atual_bin_recurso1.tif')
plot(maxent.pred.bin)
table(getValues(maxent.pred.bin))

### Training and test datasets ###
set.seed(179)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model2<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model2@results
maxent.predict2<-predict(maxent.model2,clima.mod)
writeRaster(maxent.predict2,'maxent_atual_sci_rec_red2.tif')
maxent.eval2<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model2)
maxent.eval2
maxent.th2<-threshold(maxent.eval2)
maxent.th2$spec_sens
maxent.th2$kappa
maxent.th2$prevalence
maxent.th2$sensitivity
maxent.pred.bin2<-maxent.predict2
maxent.pred.bin2[maxent.pred.bin2>=maxent.th2$spec_sens]=1
maxent.pred.bin2[maxent.pred.bin2<maxent.th2$spec_sens]=0
writeRaster(maxent.pred.bin2, 'maxent_atual_bin_recurso2.tif')
table(getValues(maxent.pred.bin2))

#### Training and test datasets 
set.seed(272)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model3<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model3@results
maxent.predict3<-predict(maxent.model3,clima.mod)
writeRaster(maxent.predict3,'maxent_atual_sci_rec_red3.tif')
maxent.eval3<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model3)
maxent.eval3
maxent.th3<-threshold(maxent.eval3)
maxent.th3$spec_sens
maxent.th3$kappa
maxent.th3$prevalence
maxent.th3$sensitivity
maxent.pred.bin3<-maxent.predict3
maxent.pred.bin3[maxent.pred.bin3>=maxent.th3$spec_sens]=1
maxent.pred.bin3[maxent.pred.bin3<maxent.th3$spec_sens]=0
writeRaster(maxent.pred.bin3, 'maxent_atual_bin_recurso3.tif')
table(getValues(maxent.pred.bin3))

#### Training and test datasets
set.seed(73)

sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model4<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model4@results
maxent.predict4<-predict(maxent.model4,clima.mod)
writeRaster(maxent.predict4,'maxent_atual_sci_rec_red4.tif')
maxent.eval4<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model4)
maxent.eval4
maxent.th4<-threshold(maxent.eval4)
maxent.th4$spec_sens
maxent.th4$kappa
maxent.th4$prevalence
maxent.th4$sensitivity
maxent.pred.bin4<-maxent.predict4
maxent.pred.bin4[maxent.pred.bin4>=maxent.th4$spec_sens]=1
maxent.pred.bin4[maxent.pred.bin4<maxent.th4$spec_sens]=0
writeRaster(maxent.pred.bin4, 'maxent_atual_bin_recurso4.tif')
table(getValues(maxent.pred.bin4))

### Training and test datasets
set.seed(712)

sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model5<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model5@results
maxent.predict5<-predict(maxent.model5,clima.mod)
writeRaster(maxent.predict5,'maxent_atual_sci_rec_red5.tif')
maxent.eval5<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model5)
maxent.eval5
maxent.th5<-threshold(maxent.eval5)
maxent.th5$spec_sens
maxent.th5$kappa
maxent.th5$prevalence
maxent.th5$sensitivity
maxent.pred.bin5<-maxent.predict5
maxent.pred.bin5[maxent.pred.bin5>=maxent.th5$spec_sens]=1
maxent.pred.bin5[maxent.pred.bin5<maxent.th5$spec_sens]=0
writeRaster(maxent.pred.bin5, 'maxent_atual_bin_recurso5.tif')
table(getValues(maxent.pred.bin5))

### Training and test datasets
set.seed(658)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model6<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model6@results
maxent.predict6<-predict(maxent.model6,clima.mod)
writeRaster(maxent.predict6,'maxent_atual_sci_rec_red6.tif')
maxent.eval6<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model6)
maxent.eval6
maxent.th6<-threshold(maxent.eval6)
maxent.th6$spec_sens
maxent.th6$kappa
maxent.th6$prevalence
maxent.th6$sensitivity
maxent.pred.bin6<-maxent.predict6
maxent.pred.bin6[maxent.pred.bin6>=maxent.th6$spec_sens]=1
maxent.pred.bin6[maxent.pred.bin6<maxent.th6$spec_sens]=0
writeRaster(maxent.pred.bin6, 'maxent_atual_bin_recurso6.tif')
table(getValues(maxent.pred.bin6))

### Training and test datasets

set.seed(417)

sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model7<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model7@results
maxent.predict7<-predict(maxent.model7,clima.mod)
writeRaster(maxent.predict7,'maxent_atual_sci_rec_red7.tif')
maxent.eval7<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model7)
maxent.eval7
maxent.th7<-threshold(maxent.eval7)
maxent.th7$spec_sens
maxent.th7$kappa
maxent.th7$prevalence
maxent.th7$sensitivity
maxent.pred.bin7<-maxent.predict7
maxent.pred.bin7[maxent.pred.bin7>=maxent.th7$spec_sens]=1
maxent.pred.bin7[maxent.pred.bin7<maxent.th7$spec_sens]=0
writeRaster(maxent.pred.bin7, 'maxent_atual_bin_recurso7.tif')
table(getValues(maxent.pred.bin7))

### Training and test datasets 
set.seed(625)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model8<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model8@results
maxent.predict8<-predict(maxent.model8,clima.mod)
writeRaster(maxent.predict8,'maxent_atual_sci_rec_red8.tif')
maxent.eval8<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model8)
maxent.eval8
maxent.th8<-threshold(maxent.eval8)
maxent.th8$spec_sens
maxent.th8$kappa
maxent.th8$prevalence
maxent.th8$sensitivity
maxent.pred.bin8<-maxent.predict8
maxent.pred.bin8[maxent.pred.bin8>=maxent.th8$spec_sens]=1
maxent.pred.bin8[maxent.pred.bin8<maxent.th8$spec_sens]=0
writeRaster(maxent.pred.bin8, 'maxent_atual_bin_recurso8.tif')
table(getValues(maxent.pred.bin8))

### Training and test datasets
set.seed(867)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####
maxent.model9<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model9@results
maxent.predict9<-predict(maxent.model9,clima.mod)
writeRaster(maxent.predict9,'maxent_atual_sci_rec_red9.tif')
maxent.eval9<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model9)
maxent.eval9
maxent.th9<-threshold(maxent.eval9)
maxent.th9$spec_sens
maxent.th9$kappa
maxent.th9$prevalence
maxent.th9$sensitivity
maxent.pred.bin9<-maxent.predict9
maxent.pred.bin9[maxent.pred.bin9>=maxent.th9$spec_sens]=1
maxent.pred.bin9[maxent.pred.bin9<maxent.th9$spec_sens]=0
writeRaster(maxent.pred.bin9, 'maxent_atual_bin_recurso9.tif')
table(getValues(maxent.pred.bin9))

### Training and test datasets
set.seed(941)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model10<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model10@results
maxent.predict10<-predict(maxent.model10,clima.mod)
writeRaster(maxent.predict10,'maxent_atual_sci_rec_red10.tif')
maxent.eval10<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model10)
maxent.eval10
maxent.th10<-threshold(maxent.eval10)
maxent.th10$spec_sens
maxent.th10$kappa
maxent.th10$prevalence
maxent.th10$sensitivity
maxent.pred.bin10<-maxent.predict10
maxent.pred.bin10[maxent.pred.bin10>=maxent.th10$spec_sens]=1
maxent.pred.bin10[maxent.pred.bin10<maxent.th10$spec_sens]=0
writeRaster(maxent.pred.bin10, 'maxent_atual_bin_recurso10.tif')
table(getValues(maxent.pred.bin10))



#### MODELLING WITHOUT BIOTIC PREDICTOR

#### Repeating the proccess of model with biotic predictor, but without it ####

clima <- stack(bio2,bio3,bio8,bio15,bio18,bio19)
clima.df <- as.data.frame(clima, xy=T)
clima.df<- na.omit(clima.df)
gridded(clima.df) <- ~ x+y
clima.st<-stack(clima.df)

## cleaning data for next background/pseudo-ausence ##

clima.mod<-clima.st
clim.pts<- extract(clima.mod, pts.geo, cellnumbers=T)
sciurus.var<- cbind(pontos[,1:2], clim.pts)

### excluding points in the same pixel ###

duplicated(sciurus.var[,'cells'])
dup<- which(duplicated(sciurus.var[,'cells']) == TRUE)
dup
sciurus.var <- sciurus.var [-dup,]
dim(sciurus.var)
sciurus.var.limpo <- na.omit(sciurus.var)
write.table(sciurus.var.limpo, 'sciurus_clim.txt')

##### making background/pseudo-ausence ####

sciurus.oc<-read.table('sciurus_clim.txt', h=T)

clima.mod.df<-na.omit(as.data.frame(clima.mod, xy=T))

back.id<-sample(1:nrow(clima.mod.df), nrow(sciurus.oc))

back<-clima.mod.df[back.id,]

write.table(back, 'background_sciurus.txt')

points(back[,'x'], back[,'y'], col='red')
points(sciurus.oc[,'Longitude'], sciurus.oc[,'Latitude'])

back<-read.table('background_sciurus.txt')

## cross validation = 10 
Sys.setenv(NOAWT=T)
library(rJava)

### Training and test datasets
set.seed(774)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model1<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model1@results
maxent.predict<-predict(maxent.model1,clima.mod)
writeRaster(maxent.predict,'maxent_atual_sci_semrec_red1.tif', overwrite=T)
maxent.eval<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model1)
maxent.eval
maxent.th<-threshold(maxent.eval)
maxent.th$spec_sens
maxent.th$kappa
maxent.th$prevalence
maxent.th$sensitivity
maxent.pred.bin<-maxent.predict
maxent.pred.bin[maxent.pred.bin>=maxent.th$spec_sens]=1
maxent.pred.bin[maxent.pred.bin<maxent.th$spec_sens]=0
writeRaster(maxent.pred.bin, 'maxent_atual_bin_semrec.tif')
table(getValues(maxent.pred.bin))

### Training and test datasets
set.seed(516)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model2<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model2@results
maxent.predict2<-predict(maxent.model2,clima.mod)
writeRaster(maxent.predict2,'maxent_atual_sci_semrec_red2.tif')
maxent.eval2<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model2)
maxent.eval2
maxent.th2<-threshold(maxent.eval2)
maxent.th2$spec_sens
maxent.th2$kappa
maxent.th2$prevalence
maxent.th2$sensitivity
maxent.pred.bin2<-maxent.predict2
maxent.pred.bin2[maxent.pred.bin2>=maxent.th2$spec_sens]=1
maxent.pred.bin2[maxent.pred.bin2<maxent.th2$spec_sens]=0
writeRaster(maxent.pred.bin, 'maxent_atual_bin_semrec2.tif')
table(getValues(maxent.pred.bin2))

### Training and test datasets
set.seed(368)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model3<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model3@results
maxent.predict3<-predict(maxent.model3,clima.mod)
writeRaster(maxent.predict3,'maxent_atual_sci_semrec_red3.tif')
maxent.eval3<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model3)
maxent.eval3
maxent.th3<-threshold(maxent.eval3)
maxent.th3$spec_sens
maxent.th3$kappa
maxent.th3$prevalence
maxent.th3$sensitivity
maxent.pred.bin3<-maxent.predict3
maxent.pred.bin3[maxent.pred.bin3>=maxent.th3$spec_sens]=1
maxent.pred.bin3[maxent.pred.bin3<maxent.th3$spec_sens]=0
writeRaster(maxent.pred.bin3, 'maxent_atual_bin_semrec3.tif')
table(getValues(maxent.pred.bin3))

### Training and test datasets
set.seed(290)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model4<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model4@results
maxent.predict4<-predict(maxent.model4,clima.mod)
writeRaster(maxent.predict4,'maxent_atual_sci_semrec_red4.tif')
maxent.eval4<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model4)
maxent.eval4
maxent.th4<-threshold(maxent.eval4)
maxent.th4$spec_sens
maxent.th4$kappa
maxent.th4$prevalence
maxent.th4$sensitivity
maxent.pred.bin4<-maxent.predict4
maxent.pred.bin4[maxent.pred.bin4>=maxent.th4$spec_sens]=1
maxent.pred.bin4[maxent.pred.bin4<maxent.th4$spec_sens]=0
writeRaster(maxent.pred.bin4, 'maxent_atual_bin_semrec4.tif')
table(getValues(maxent.pred.bin4))

### Training and test datasets
set.seed(796)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent

maxent.model5<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model5@results
maxent.predict5<-predict(maxent.model5,clima.mod)
writeRaster(maxent.predict5,'maxent_atual_sci_semrec_red5.tif')
maxent.eval5<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model5)
maxent.eval5
maxent.th5<-threshold(maxent.eval5)
maxent.th5$spec_sens
maxent.th5$kappa
maxent.th5$prevalence
maxent.th5$sensitivity
maxent.pred.bin5<-maxent.predict5
maxent.pred.bin5[maxent.pred.bin5>=maxent.th5$spec_sens]=1
maxent.pred.bin5[maxent.pred.bin5<maxent.th5$spec_sens]=0
writeRaster(maxent.pred.bin5, 'maxent_atual_bin_semrec5.tif')
table(getValues(maxent.pred.bin5))

### Training and test datasets
set.seed(994)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent

maxent.model6<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model6@results
maxent.predict6<-predict(maxent.model6,clima.mod)
writeRaster(maxent.predict6,'maxent_atual_sci_semrec_red6.tif')
maxent.eval6<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model6)
maxent.eval6
maxent.th6<-threshold(maxent.eval6)
maxent.th6$spec_sens
maxent.th6$kappa
maxent.th6$prevalence
maxent.th6$sensitivity
maxent.pred.bin6<-maxent.predict6
maxent.pred.bin6[maxent.pred.bin6>=maxent.th6$spec_sens]=1
maxent.pred.bin6[maxent.pred.bin6<maxent.th6$spec_sens]=0
writeRaster(maxent.pred.bin6, 'maxent_atual_bin_semrec6.tif')
table(getValues(maxent.pred.bin6))

### Training and test datasets
set.seed(705)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent

maxent.model7<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model7@results
maxent.predict7<-predict(maxent.model7,clima.mod)
writeRaster(maxent.predict7,'maxent_atual_sci_semrec_red7.tif', overwrite=T)
maxent.eval7<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model7)
maxent.eval7
maxent.th7<-threshold(maxent.eval7)
maxent.th7$spec_sens
maxent.th7$kappa
maxent.th7$prevalence
maxent.th7$sensitivity
maxent.pred.bin7<-maxent.predict7
maxent.pred.bin7[maxent.pred.bin7>=maxent.th7$spec_sens]=1
maxent.pred.bin7[maxent.pred.bin7<maxent.th7$spec_sens]=0
writeRaster(maxent.pred.bin7, 'maxent_atual_bin_semrec7.tif')
table(getValues(maxent.pred.bin7))

### Training and test datasets
set.seed(533)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent ####

maxent.model8<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model8@results
maxent.predict8<-predict(maxent.model8,clima.mod)
writeRaster(maxent.predict8,'maxent_atual_sci_semrec_red8.tif')
maxent.eval8<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model8)
maxent.eval8
maxent.th8<-threshold(maxent.eval8)
maxent.th8$spec_sens
maxent.th8$kappa
maxent.th8$prevalence
maxent.th8$sensitivity
maxent.pred.bin8<-maxent.predict8
maxent.pred.bin8[maxent.pred.bin8>=maxent.th8$spec_sens]=1
maxent.pred.bin8[maxent.pred.bin8<maxent.th8$spec_sens]=0
writeRaster(maxent.pred.bin8, 'maxent_atual_bin_semrec8.tif')
table(getValues(maxent.pred.bin8))

### Training and test datasets
set.seed(393)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent

maxent.model9<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model9@results
maxent.predict9<-predict(maxent.model9,clima.mod)
writeRaster(maxent.predict9,'maxent_atual_sci_semrec_red9.tif')
maxent.eval9<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model9)
maxent.eval9
maxent.th9<-threshold(maxent.eval9)
maxent.th9$spec_sens
maxent.th9$kappa
maxent.th9$prevalence
maxent.th9$sensitivity
maxent.pred.bin9<-maxent.predict9
maxent.pred.bin9[maxent.pred.bin9>=maxent.th9$spec_sens]=1
maxent.pred.bin9[maxent.pred.bin9<maxent.th9$spec_sens]=0
writeRaster(maxent.pred.bin9, 'maxent_atual_bin_semrec9.tif')
table(getValues(maxent.pred.bin9))

### Training and test datasets
set.seed(165)
sample.ocor <- sample(1:nrow(sciurus.oc), round(0.75*nrow(sciurus.oc)))
sample.back <- sample(1:nrow(back), round(0.75*nrow(back)))

treino <- prepareData(x=clima.mod, p=sciurus.oc[sample.ocor,1:2], b=back[sample.back,1:2])

teste<- prepareData (x= clima.mod, p=sciurus.oc[- sample.ocor, 1:2], b=back[-sample.back, 1:2])

#### Modelling with maxent

maxent.model10<-maxent(teste[,2:7],p=teste[,'pb'])
maxent.model10@results
maxent.predict10<-predict(maxent.model10,clima.mod)
writeRaster(maxent.predict10,'maxent_atual_sci_semrec_red10.tif')
maxent.eval10<-evaluate(p=teste[teste[,'pb']==1,-1], a=teste[teste['pb']==0,-1], model=maxent.model10)
maxent.eval10
maxent.th10<-threshold(maxent.eval10)
maxent.th10$spec_sens
maxent.th10$kappa
maxent.th10$prevalence
maxent.th10$sensitivity
maxent.pred.bin10<-maxent.predict10
maxent.pred.bin10[maxent.pred.bin10>=maxent.th10$spec_sens]=1
maxent.pred.bin10[maxent.pred.bin10<maxent.th10$spec_sens]=0
writeRaster(maxent.pred.bin10, 'maxent_atual_bin_semrec10.tif')
table(getValues(maxent.pred.bin10))


#### Making ensembles of suitability rasters ####

max1.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red1.tif')
max2.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red2.tif')
max3.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red3.tif')
max4.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red4.tif')
max5.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red5.tif')
max6.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red6.tif')
max7.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red7.tif')
max8.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red8.tif')
max9.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red9.tif')
max10.comrec<-raster('Modelo_com_recurso/maxent_atual_sci_rec_red10.tif')

beginCluster()
ens.mean.comrec<-sum(max1.comrec,max2.comrec,max3.comrec,max4.comrec,max5.comrec,max6.comrec,max7.comrec,max8.comrec,max9.comrec, max10.comrec)
ens.mean.comrec<-(ens.mean.comrec/10)

plot(ens.mean.comrec, main='Ensemble - with resource')
scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')
writeRaster(ens.mean.comrec,'Ensemble_sciurus_withresource_final.tif',overwrite=T)

max1.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red1.tif')
max2.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red2.tif')
max3.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red3.tif')
max4.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red4.tif')
max5.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red5.tif')
max6.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red6.tif')
max7.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red7.tif')
max8.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red8.tif')
max9.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red9.tif')
max10.semrec<-raster('Modelo_sem_recurso/maxent_atual_sci_semrec_red10.tif')


beginCluster()
ens.mean.semrec<-sum(max1.semrec,max2.semrec,max3.semrec,max4.semrec,max5.semrec,max6.semrec,max7.semrec,max8.semrec,max9.semrec, max10.semrec)
ens.mean.semrec<-(ens.mean.semrec/10)

plot(ens.mean.semrec, main='Ensemble - without resource')
scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')

writeRaster(ens.mean.semrec,'Ensemble_sciurus_withoutresource_final.tif',overwrite=T)


#### Making ensembles with presence/ausence rasters ####

max1.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso1.tif')
max2.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso2.tif')
max3.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso3.tif')
max4.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso4.tif')
max5.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso5.tif')
max6.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso6.tif')
max7.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso7.tif')
max8.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso8.tif')
max9.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso9.tif')
max10.comrec.bin<-raster('modelo_com_recurso_binario/maxent_atual_bin_recurso10.tif')

beginCluster()
ens.mean.comrec.bin<-sum(max1.comrec.bin,max2.comrec.bin,max3.comrec.bin,max4.comrec.bin,max5.comrec.bin,max6.comrec.bin,max7.comrec.bin,max8.comrec.bin,max9.comrec.bin,max10.comrec.bin)
ens.mean.comrec.bin<-(ens.mean.comrec.bin/10)

plot(ens.mean.comrec.bin, main='Ensemble - with resources(pres/abs)')
scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')

writeRaster(ens.mean.comrec.bin,'Ensemble_sciurus_comrecurso_bin.tif', overwrite=T)

max1.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec.tif')
max2.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec2.tif')
max3.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec3.tif')
max4.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec4.tif')
max5.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec5.tif')
max6.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec6.tif')
max7.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec7.tif')
max8.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec8.tif')
max9.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec9.tif')
max10.semrec.bin<-raster('modelo_sem_recurso_binario/maxent_atual_bin_semrec10.tif')

beginCluster()
ens.mean.semrec.bin<-sum(max1.semrec.bin,max2.semrec.bin,max3.semrec.bin,max4.semrec.bin,max5.semrec.bin,max6.semrec.bin,max7.semrec.bin,max8.semrec.bin,max9.semrec.bin,max10.semrec.bin)
ens.mean.semrec.bin<-(ens.mean.semrec.bin/10)

plot(ens.mean.semrec.bin, main='Ensemble - without resources (pres/abs)')
scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')

writeRaster(ens.mean.semrec.bin,'Ensemble_sciurus_semrecurso_bin.tif', overwrite=T)


#### Multiplication of ensembles ####
## Suitability rasters
ensemble.com.recurso<-raster('ensemble_raster_final/Ensemble_sciurus_withresource_final.tif')
ensemble.sem.recurso<-raster('ensemble_raster_final/Ensemble_sciurus_withoutresource_final.tif')
ensemble.com.recurso<-ensemble.com.recurso*100
ensemble.sem.recurso<-ensemble.sem.recurso*10
e<-extent(-91.66681, -34.83347, -55.90014, 12.59153)
s<- raster(e, nrows=8219, ncols= 6825, crs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
ensemble.com.recurso<-resample(ensemble.com.recurso,s, method='ngb')
ensemble.sem.recurso<-resample(ensemble.sem.recurso,s, method='ngb')

subtraction.com.sem<-overlay(ensemble.com.recurso, ensemble.sem.recurso, fun=function(x,y){return(x-y)})
writeRaster(subtraction.com.sem, 'Subtraction_with-without_ensemble_continuo', overwrite=T)
plot(subtraction.com.sem, main='Subtraction (with-without)')
scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')


#### plotting

library(raster)
com.recurso<-raster('Ensemble_sciurus_withresource_final.tif')
c.recurso<-com.recurso*100
sem.recurso<-raster('Ensemble_sciurus_withoutresource_final.tif')

ensemble.com.recurso.bin<-raster('ensemble_raster_final/ensemble_sciurus_comrecurso_bin.tif')
ensemble.sem.recurso.bin<-raster('ensemble_raster_final/ensemble_sciurus_semrecurso_bin.tif')
e<-extent(-91.66681, -34.83347, -55.90014, 12.59153)
s<- raster(e, nrows=8219, ncols= 6825, crs='+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0')
c.recurso<-resample(com.recurso,s, method='ngb')
sem.recurso<-resample(sem.recurso,s, method='ngb')

subtraction.com.sem<-overlay(c.recurso, sem.recurso, fun=function(x,y){return(x*y)})
writeRaster(subtraction.com.sem, 'Subtraction_with-without_bin_ensemble')
plot(subtraction.com.sem, main='Subtraction (with-without) - pres/abs')
scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')

#### plotting

library(rgdal)
library(rasterVis)
library(RColorBrewer)
library(latticeExtra)

colr <- colorRampPalette(brewer.pal(11, 'YlGn'))
south.america<-readOGR('South_America.shp')
south.america <- spTransform(south.america, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
semrec.bin<-raster('Ensemble_sciurus_comrecurso_bin.tif')
crs(semrec.bin) <-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
comrec.bin<-raster('Ensemble_sciurus_semrecurso_bin.tif')
com.recurso
sem.recurso

levelplot(semrec.bin, main='Model without resource',
          margin=FALSE,                       
          colorkey=list(
            labels=list(at=c(0,0.2,0.4,0.6,0.8), font=4),
            axis.line=list(col='black')       
          ),    
          par.settings=list(
            axis.line=list(col='transparent') 
          ),
          scales=list(draw=TRUE),            
          col.regions=colr,                   
          at=seq(0,1, by=0.2)) +
  layer(sp.polygons(south.america))

levelplot(c.recurso, main='Model with resource',
          margin=FALSE,                       
          colorkey=list(
            labels=list(at=c(0,0.2,0.4,0.6,0.8), font=4),
            axis.line=list(col='black')       
          ),    
          par.settings=list(
            axis.line=list(col='transparent') 
          ),
          scales=list(draw=TRUE),            
          col.regions=colr,                   
          at=seq(0,1, by=0.2)) +
  layer(sp.polygons(south.america))

levelplot(sem.recurso, main='Model with resource',
          margin=FALSE,                       
          colorkey=list(
            labels=list(at=c(0,0.2,0.4,0.6,0.8), font=4),
            axis.line=list(col='black')       
          ),    
          par.settings=list(
            axis.line=list(col='transparent') 
          ),
          scales=list(draw=TRUE),            
          col.regions=colr,                   
          at=seq(0,1, by=0.2)) +
  layer(sp.polygons(south.america))

plotando<-levelplot(subtraction.com.sem, main='Concordance between maps',
                    margin=FALSE,                       
                    colorkey=list(
                      labels=list(at=c(0,20,40,60,80), font=4),
                      axis.line=list(col='black')       
                    ),    
                    par.settings=list(
                      axis.line=list(col='transparent') 
                    ),
                    scales=list(draw=TRUE),            
                    col.regions=colr,                   
                    at=seq(0,80, by=20)) +
  layer(sp.polygons(south.america))

plotando 

#### plotting a scalebar in final maps 

scalebar(2000, xy=c(-120, -50), type='bar', divs=2, belo = 'km')

### Barplot of contribution of variables, first without biotic predictor

without_syagrus<-read.csv('contribution_without_syagrus.csv')

# Change the colors manually

library(ggplot2)
p <- ggplot(data=without_syagrus, aes(x=reorder(predictor, -contribution), y=contribution, fill=predictor) +
  geom_bar(stat="identity", color="black", position=position_dodge())+
  theme_minimal()
# Use custom colors

p + scale_fill_hue(l=40, c=35)
# Use brewer color palettes
#p + scale_fill_brewer(palette="Greys")

p + labs(tittle = 'Variables contribution', x= 'Predictor variables', y='Contribution', tag = 'A')+
theme(axis.title.x = element_text(size=14))+
  theme(axis.title.y = element_text(size=14))

### Now with biotic predictor

with_syagrus<-read.csv('contribution_with_syagrus.csv')

p <- ggplot(data=with_syagrus, aes(x=reorder(Predictor, -Contribution), y=Contribution, fill=Predictor)) +
geom_bar(stat="identity", color="black", position=position_dodge())+
  theme_minimal() 

# Use custom colors
            
p + scale_fill_hue(l=40, c=35)
            
# Use brewer color palettes
#p + scale_fill_brewer(palette="Greys")
            
p + labs(tittle = 'Variables contribution', x= 'Predictor variables', y='Contribution', tag = 'B')+
theme(axis.title.x = element_text(size=14))+
theme(axis.title.y = element_text(size=14))+
  scale_fill_hue(l=40, c=35)

### Making boxplot of presence pixels
library(ggplot2)
presence<-read.csv('boxplot.csv')

box<-ggplot(presence, aes(x=model, y=presence_pixels, color=model)) +
  geom_boxplot(coef=1e30,fill='#A4A4A4', color="black")+
  theme_classic()
box + labs(x= 'Model', y='Number of Presence Pixels')+
  theme(axis.title.x = element_text(size=14))+
  theme(axis.title.y = element_text(size=14))