Cervus_elaphus

library(dismo)

set.seed(10)

group_pres<-kfold(spdf, k)

if (!file.exists(paste(species_directory, "k_folds_presence_2006.csv", sep="/"))){
  write.csv(group_pres, paste(species_directory, "k_folds_presence_2006.csv", sep="/"))
}

group_pres<-read.csv(paste(species_directory, "k_folds_presence_2006.csv", sep="/"))
group_pres<-group_pres[,-1]  #all presence points

sp_train<-sp[group_pres!=1,] #presence points split into training and test
sp_test<-sp[group_pres ==1,]

targetRecords <- occ_search(
  scientificName=c(
    'Dama_dama',
    'Alces_alces',
    'Capreolus capreolus',
     'Sus scrofa',
    'Martes_martes',
    'Meles_meles',
    "Vulpes_vulpes",
    "Mustela_erminea",
    "Lutra_lutra"
 
  ),
  hasCoordinate=TRUE,
  year='2000,2016',
  hasGeospatialIssue=FALSE,
  fields='minimal',
  decimalLongitude='-11,33',
  decimalLatitude='35,68',
  limit=11000
)

if (exists('targetSites')) rm(targetSites)

for (i in 1:length(targetRecords)) {

  thisGenusCoords <- data.frame(
    longitude=targetRecords[[i]]$data$decimalLongitude,
    latitude=targetRecords[[i]]$data$decimalLatitude
  )

  targetSites <- if (exists('targetSites')) {
      
    rbind(targetSites, thisGenusCoords)
    
  } else {
    
    thisGenusCoords
    
  }
  
}

targetEnv <- as.data.frame(extract(pred_crop, targetSites))

outside <- which(is.na(rowSums(targetEnv))) 

targetSites <- targetSites[-outside, ]
targetEnv <- targetEnv[-outside, ]

targetBg <- cbind(targetSites, targetEnv)
names(targetBg)[1:2] <- c('longitude', 'latitude')

nrow(targetBg)

targetBg <- targetBg[sample(1:nrow(targetBg), 10000), ]
nrow(targetBg)

#write.csv(targetBg, paste(species_directory,"background_data_et_ungulates.csv", sep="/"))

bg_pse<-read.csv(paste(species_directory,"background_data_et_ungulates.csv", sep="/"))

bg_pse<-bg_pse[,-1]

k<-4

targetBg_k<-kfold(bg_pse, k)

#write.csv(targetBg_k, paste(species_directory, "kfolds_background_data_et_ungulates.csv", sep="/"))

group_back<-read.csv(paste(species_directory,"kfolds_background_data_et_ungulates.csv", sep="/"))

group_back<-group_back[,-1]

sp_backg_train <- bg_pse[group_back != 1,1:2] #background points split into training and test
sp_backg_test <- bg_pse[group_back == 1,1:2]
colnames(sp_backg_train)<-c("lon", "lat")
colnames(sp_backg_test)<-c("lon", "lat")

sp_train<-sp_train@coords #presence training points
colnames(sp_train)<-c("lon", "lat")
train <- rbind(sp_train, sp_backg_train)  #presence and background training points combined

sp_test<-sp_test@coords #presence test points
colnames(sp_test)<-c("lon", "lat")
test<-rbind(sp_test, sp_backg_test)    #presence and background test points combined

pb_train <- c(rep(1, nrow(sp_train)), rep(0, nrow(sp_backg_train)))
envtrain <- extract(pred_nf, train)
envtrain <- data.frame( cbind(pa=pb_train, envtrain) ) #training points with environmental  variables

envtrain<-na.omit(envtrain)

pb_test<-c(rep(1, nrow(sp_test)), rep(0, nrow(sp_backg_test)))
envtest <- extract(pred_nf, test)
envtest <- data.frame( cbind(pa=pb_test, envtest) ) #test points with environmental  variables
envtest<-na.omit(envtest)

env_all<-rbind(envtrain, envtest) #all of the training and test data with environmental variables

env_pres<-extract(pred_nf, sp)
pa<-rep(1, nrow(env_pres))
env_pres<-data.frame(pa, env_pres)  #presence points with environmental variables
colnames(env_pres)<-c("pa", bioclim_names, "Hyde_Year_2005")
env_pres_xy<-data.frame(sp@coords, env_pres)
env_pres_xy<-na.omit(env_pres_xy)

env_back<-extract(pred_nf, bg_pse[,1:2])
ap<-rep(0, nrow(env_back))
env_back<-data.frame(ap, env_back) #background points with environmental variables
colnames(env_back)<-c("pa",bioclim_names, "Hyde_Year_2005")
env_back_xy<-data.frame(bg_pse[,1:2], env_back)
env_back_xy<-na.omit(env_back_xy)

group_pres<-group_pres[1:nrow(env_pres_xy)]
group_back<-group_back[1:nrow(env_back_xy)]

evl_bc<- list()
tss_bc<-list()
for (i in 1:k){
  pres_train<-sp[group_pres!=i,]
  pres_test<-sp[group_pres==i,]
  backg_test<-bg_pse[group_back==i,1:2]
  bc <- bioclim(pred_nf,pres_train)
  evl_bc[[i]] <- dismo:::evaluate(pres_test, backg_test, bc,pred_nf,type="response")#test presence, test absence, model, predictor variables
     #print(i)
}

auc_bc <- sapply( evl_bc, function(x){slot(x, "auc")} )
print(auc_bc)

bc_auc<-mean(auc_bc)

#library(mgcv)

evl_gam<- list()
tss_gam<-list()
for (i in 1:k){
  pres_train<-env_pres_xy[group_pres!=i ,-c(1,2)]
  pres_test<-env_pres_xy[(group_pres==i) ,-c(1,2)]
  back_test<-env_back_xy[(group_back==i),-c(1,2)]
  back_train<-env_back_xy[(group_back!=i),-c(1,2)]
  envtrain<-rbind(pres_train, back_train)
  gm1<-mgcv:::gam(pa~ s(Bio_2_2006_2016_average)+s(Bio_4_2006_2016_average)+ s(Bio_5_2006_2016_average)+ 
           s(Bio_6_2006_2016_average)+ s(Bio_7_2006_2016_average)+ s(Bio_10_2006_2016_average)+ 
           s(Bio_11_2006_2016_average) + s(Hyde_Year_2005), family = binomial(link = "logit"),data=envtrain)


  evl_gam[[i]] <- dismo:::evaluate(p = pres_test, a = back_test,model= gm1,type="response")
    }

auc_gam <- sapply( evl_gam, function(x){slot(x, "auc")} )
print(auc_gam)

gam_auc<-mean(auc_gam)

model<-pa ~  Bio_2_2006_2016_average+Bio_4_2006_2016_average+ 
  Bio_5_2006_2016_average+ Bio_6_2006_2016_average+ Bio_7_2006_2016_average+ 
  Bio_10_2006_2016_average+ Bio_11_2006_2016_average + Hyde_Year_2005


evl_rf<- list()
for (i in 1:k){
  pres_train<-env_pres_xy[group_pres!=i,-c(1,2)]
  pres_test<-env_pres_xy[(group_pres==i) ,-c(1,2)]
  back_test<-env_back_xy[(group_back==i),-c(1,2)]
  back_train<-env_back_xy[group_back !=i,-c(1,2)]
  envtrain<-rbind(pres_train, back_train)
  rf1 <- randomForest(model, data=envtrain)
  evl_rf[[i]] <- dismo:::evaluate(pres_test, back_test, rf1,type="response")
  }

auc_rf <- sapply( evl_rf, function(x){slot(x, "auc")} )
print(auc_rf)

rf_auc<-mean(auc_rf)

#total models
library(randomForest)
library(boot)
bc <- bioclim(pred_crop, sp)
gm1<-gam(pa~ s(Bio_2_2006_2016_average)+s(Bio_4_2006_2016_average)+ s(Bio_5_2006_2016_average)+ 
           s(Bio_6_2006_2016_average)+ s(Bio_7_2006_2016_average)+ s(Bio_10_2006_2016_average)+ 
           s(Bio_11_2006_2016_average) + s(Hyde_Year_2005), family = binomial(link = "logit"),data=env_all)
rf1 <- randomForest(model, data=env_all)

pb <- predict(pred_crop, bc, ext=e_iucn, progress='') #bioclim predict
pg <- predict(pred_crop, gm1, ext=e_iucn) #gam predict
pgl<-raster:::calc(pg, fun=function(x){ exp(x)/(1+exp(x))}) #backtransforming from logit space
pr <- predict(pred_crop, rf1, ext=e_iucn) #random forest predict

models <- stack(pb, pgl, pr)
names(models) <- c("bioclim", "gam", "random forest")
plot(models)

auc<-c(bc_auc, gam_auc, rf_auc)
w <- (auc-0.5)^2
wm <- weighted.mean( models[[c("bioclim", "gam", "random.forest")]], w)
plot(wm, main="Ensemble model")
points(sp)

abs_xy<-env_back_xy[,c(1,2)]

ens<-ensemble_evaluate(sp,abs_xy , wm)
thresh<-threshold(ens, stat="spec_sens")

thresh

bc <- bioclim(pred_nf, sp)
gm1<-mgcv:::gam(pa~ s(Bio_2_2006_2016_average)+ s(Bio_4_2006_2016_average)+ s(Bio_5_2006_2016_average)+ s(Bio_6_2006_2016_average)+ s(Bio_7_2006_2016_average)+s(Bio_10_2006_2016_average)+ s(Bio_11_2006_2016_average) + s(Hyde_Year_2005), data=env_all)
rf1 <- randomForest:::randomForest(model, data=env_all)

for (i in 1:length(years)){
 
  if (Europe_only){
      pred_nf<-stack(paste(wd, "/Bioclim/", years[i], "_bioclim_variable_stack.tif", sep="" ))
  }else{
     pred_nf<-stack(paste(wd, "/Bioclim/Global/", years[i], "_bioclim_variable_stack.tif", sep="" ))
  }

  pred_nf<-pred_nf[[bioclim_layers]]
  names(pred_nf)<-bioclim_names
  hyde_yr<-raster(paste(wd, "/hyde/Hyde_Year_", years[i], ".tif",sep=""))
  names(hyde_yr)<-"Hyde_Year_2005"
  pred_nf<-stack(pred_nf, hyde_yr)
  
  pb <- predict(pred_nf, bc, ext=e_big, progress='')
  pg <- predict(pred_nf, gm1, ext=e_big) #gam predict
  pgl<-raster:::calc(pg, fun=function(x){ exp(x)/(1+exp(x))})
  pr <- predict(pred_nf, rf1, ext=e_big)

  models <- stack(pb, pgl, pr)
  names(models) <- c("bioclim", "gam", "random forest")
  wm <- weighted.mean( models[[c("bioclim", "gam", "random.forest")]], w)

  pa_sss<-wm>thresh

 writeRaster(wm , paste(sdm_folder, "/hyde_weighted_ensemble_sdm_", years[i], ".tif", sep=""), overwrite=TRUE)
 writeRaster(pa_sss , paste(sdm_folder, "/hyde_pres_abs_sss_weighted_ensemble_sdm_", years[i], ".tif", sep=""), overwrite=TRUE)

  plot(pa_sss, main=years[i])
}