# ----- to be used to process  model outptut ------------------------------
# ------ load libraries to be used for processing output from  SSDM -----------------
libs <- c('tidyverse','lubridate','purrr','here','SSDM','raster','parallel',
          'rworldmap','rworldxtra','sf','raster');
lapply(libs, library,character.only=TRUE)

# ------- search and  saved SSDM object for each species ------------------
list_no_thin <- list.files('analysis/out_ssdm_no_thin/',pattern = '.RData')
list_thin <- list.files('analysis/out_ssdm_thin/',pattern = '.RData')

# ------- prepare list of species and season ------------------------------
species_season_no_thin <- str_remove_all(list_no_thin,pattern = '.RData')
species_season_thin  <-str_remove_all(list_thin,pattern = '.RData')

# ------ load actual objects ----------------------------------------------
mod_obj_no_thin <-sapply(here::here('analysis/out_ssdm_no_thin',list_no_thin), 
                        function(x) mget(load(x)),simplify = TRUE)
mod_obj_thin <- sapply(here::here('analysis/out_ssdm_thin',list_thin),
                       function(x) mget(load(x)),simplify = TRUE)
names(mod_obj_no_thin) =  species_season_no_thin
names(mod_obj_thin) = species_season_thin

# --------- list of models considered -------------------------------------
models_sel <- c('CTA', 'SVM','GLM','MARS','GBM','RF','ANN')


# ------- prepare data thining: data with  no thining ---------------------
eval_all_models_noThin = tibble(species_season = species_season_no_thin,
                         evals = vector('list',length(species_season_no_thin)))
imp_all_models_noThin = tibble(species_season = species_season_no_thin,
                         imps = vector('list',length(species_season_no_thin)))
pred_dist_all_models_noThin = tibble(species_season = species_season_no_thin,
                         pred_rasters = vector('list',length(species_season_no_thin)))

#   ---- prepare data storage: data with thining ---------------------------------------------
eval_all_models_Thin = tibble(species_season = species_season_thin,
                                evals = vector('list',length(species_season_thin)))
imp_all_models_Thin = tibble(species_season = species_season_thin,
                               imps = vector('list',length(species_season_thin)))
pred_dist_all_models_Thin = tibble(species_season = species_season_thin,
                                     pred_rasters = vector('list',length(species_season_thin)))


# --------------------- get source file for  data processing --------------
source('R/000_custom_utility_functions.R')

# ----- get all model evaluations: data with no thining ---------------------------------
for(i in 1:nrow(eval_all_models_noThin)){
  sel_obj =  mod_obj_no_thin[str_detect(names(mod_obj_no_thin),eval_all_models_noThin$species_season[i])]
  
  eval_all_models_noThin$evals[[i]] = get_mean_perf_model(model_objects = sel_obj)
  imp_all_models_noThin$imps[[i]] = get_mean_var_imp_model(model_objects = sel_obj)
  pred_dist_all_models_noThin$pred_rasters[[i]] = get_mean_dist_model(model_objects = sel_obj)  
} 

# -------- get all model evaluation impotance: data with thining ----------
for(i in 1:nrow(eval_all_models_Thin)){
  sel_obj =  mod_obj_thin[str_detect(names(mod_obj_thin),eval_all_models_Thin$species_season[i])]
  
  eval_all_models_Thin$evals[[i]] = get_mean_perf_model(model_objects = sel_obj)
  imp_all_models_Thin$imps[[i]] = get_mean_var_imp_model(model_objects = sel_obj)
  pred_dist_all_models_Thin$pred_rasters[[i]] = get_mean_dist_model(model_objects = sel_obj)  
} 


# --------------- process ensemble outputs --------------------------------
# ----- additional utility function ---------------------------------------
flattenSquareMatrix <- function(m) {
  #if( (class(m) != "matrix") | (nrow(m) != ncol(m))) stop("Must be a square matrix.") 
  if(!identical(rownames(m), colnames(m))) stop("Row and column names must be equal.")
  ut <- upper.tri(m)
  data.frame(i = rownames(m)[row(m)[ut]],
             j = rownames(m)[col(m)[ut]],
             cor=t(m)[ut],
             p=m[ut])
}

clean_mat_flatten <- function(x) {
  o_row =  rownames(x)
  o_col = colnames(x)
  n_row = str_extract_all(o_row,pattern = 'CTA|GLM|GBM|RF|ANN|MARS|SVM')%>%unlist()
  n_col= str_extract_all(o_col,pattern = 'CTA|GLM|GBM|RF|ANN|MARS|SVM')%>%unlist()
  
  rownames(x) = n_row
  colnames(x) = n_col
  
  fl_mat = flattenSquareMatrix(x)
  fl_mat
  #x
}

# --------- get all ensemble output:  no thining --------------------------
imp_ens_noThin =  purrr::map_df(mod_obj_no_thin,.f = function(x = .x) x@variable.importance%>%t()%>%
                                as.data.frame()%>%rownames_to_column(var = 'predictors')%>%
                                mutate(species_season = x@name))

eval_ens_noThin =  purrr::map_df(mod_obj_no_thin,.f = function(x = .x) x@evaluation%>%t()%>%
                                as.data.frame()%>%rownames_to_column(var = 'metrics')%>%
                                mutate(species_season = x@name))


pred_dist_ens_noThin = purrr::map_df(mod_obj_no_thin,.f = function(x = .x) tibble(species_season = x@name,
                                                        preds_dist = list(x@projection),
                                                        uncertainty = list(x@uncertainty)))

cor_mon_ens_noThin = purrr::map_df(mod_obj_no_thin,.f = function(x =.x) {
  tibble(species_season = x@name,alg_cor = list(x@algorithm.correlation%>%
                as.data.frame()%>%clean_mat_flatten()))})

system.time(all_resp_curve_noThin <- purrr::map(mod_obj_no_thin,.f = function(x = .x) get_partial_effect(model_object = x)))

# ---------------- get all ensemble output: thining -----------------------
pred_dist_ens_thin = purrr::map_df(mod_obj_thin,.f = function(x = .x) tibble(species_season = x@name,
                                                                             pred_dist = list(x@projection),
                                                                             uncertainty = list(x@uncertainty)))

cor_mon_ens_thin = purrr::map_df(mod_obj_thin,.f = function(x =.x) {
  tibble(species_season = x@name,alg_cor = list(x@algorithm.correlation%>%
                                                  as.data.frame()%>%clean_mat_flatten()))})
imp_ens_thin =  purrr::map_df(mod_obj_thin,.f = function(x = .x) x@variable.importance%>%t()%>%
                                  as.data.frame()%>%rownames_to_column(var = 'predictors')%>%
                                  mutate(species_season = x@name))

eval_ens_thin =  purrr::map_df(mod_obj_thin,.f = function(x = .x) x@evaluation%>%t()%>%
                                   as.data.frame()%>%rownames_to_column(var = 'metrics')%>%
                                   mutate(species_season = x@name))



system.time(all_resp_curve_thin <- purrr::map(mod_obj_thin,.f = function(x = .x) get_partial_effect(model_object = x)))

#purrr::map_df(all_resp_curve_thin[[1]],.f = rbind,.id = 'model')

system.time(all_occ_bkgnd_thin <- purrr::map(mod_obj_thin,.f = function(x = .x) get_occ_background(model_object = x)))
system.time(all_occ_bkgnd_noThin <- purrr::map(mod_obj_no_thin,.f = function(x = .x) get_occ_background(model_object = x)))


# -------- process all ensemble outputs -----------------------------------
# -------------- pair-wise correlation models: ensemble model -------------
cor_mon_ens_noThin <-cor_mon_ens_noThin%>%
  mutate(species_season = species_season%>%
           str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)
cor_mon_ens_thin <- cor_mon_ens_thin%>%
  mutate(species_season = species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)

cor_mon_ens_noThin%>%unnest(cols = alg_cor)%>% 
  ggplot(aes(i, j, alpha = cor)) +
  geom_tile(show.legend = FALSE) +
  geom_text(aes(label = round(cor,1)), colour = "white", alpha = 1, size = 8)+
  #labs(title =  plt_titles%>%filter(type=='all')%>%pull(plt_t))+
  facet_grid(species~ season)+
  theme_bw()

cor_mon_ens_thin%>%unnest(cols = alg_cor)%>% 
  ggplot(aes(i, j, alpha = cor)) +
  geom_tile(show.legend = FALSE) +
  geom_text(aes(label = round(cor,1)), colour = "white", alpha = 1, size = 8)+
  #labs(title =  plt_titles%>%filter(type=='all')%>%pull(plt_t))+
  facet_grid(species~ season)+
  theme_bw()

# ----------------- ensemble  variable importance ---------------------------
proc_ens_imp_thin <- imp_ens_thin%>%as_tibble()%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  mutate(predictors = str_remove_all(predictors,pattern = '_s|_w'))

proc_ens_imp_noThin <- imp_ens_noThin%>%as_tibble()%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  mutate(predictors = str_remove_all(predictors,pattern = '_s|_w'))

proc_ens_imp_noThin%>%ggplot()+
  geom_bar(aes(predictors,Axes.evaluation,fill =season),stat = 'identity',position = 'dodge')+
  labs(y =  'index of importance')+
  coord_flip()+
  facet_wrap(~species)+theme_bw()

proc_ens_imp_thin%>%ggplot()+
  geom_bar(aes(predictors,Axes.evaluation,fill =season),stat = 'identity',position = 'dodge')+
  labs(y = 'index of importance')+
  coord_flip()+
  facet_wrap(~species)+theme_bw()

proc_both_ens_thinUnThin <- rbind(proc_ens_imp_thin%>%mutate(type = 'thinned'),
                                  proc_ens_imp_noThin%>%mutate(type ='un-thinned'))

proc_both_ens_thinUnThin%>%
  ggplot()+
  geom_bar(aes(predictors,Axes.evaluation,fill =season),stat = 'identity',position = 'dodge')+
  labs(y = 'index of importance')+
  coord_flip()+
  facet_grid(type~species)+
  theme_bw()

# ----- ensemble performance --------------------------------------------------
proc_ens_perf_thin <- eval_ens_thin%>%as_tibble()%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)
proc_ens_perf_noThin <- eval_ens_noThin%>%as_tibble()%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)

proc_ens_perf_thin%>%
  filter(metrics=='AUC'|metrics=='sensitivity'|metrics=='specificity'|metrics=='Kappa')%>%
  ggplot()+
  geom_bar(aes(metrics,fp,fill =season),stat = 'identity',position = 'dodge')+
  coord_flip()+
  facet_wrap(~species)+theme_bw()

proc_ens_perf_noThin%>%
  filter(metrics=='AUC'|metrics=='sensitivity'|metrics=='specificity'|metrics=='Kappa')%>%
  ggplot()+
  geom_bar(aes(metrics,fp,fill =season),stat = 'identity',position = 'dodge')+
  coord_flip()+
  facet_wrap(~species)+theme_bw()


# -------  model specific  variable importance ----------------------------
proc_imp_model_thin <- imp_all_models_Thin%>%unnest(cols = imps)%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  mutate(predictors = str_remove_all(predictors,pattern = '_s|_w'))
proc_imp_model_noThin <- imp_all_models_noThin%>%unnest(cols = imps)%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  mutate(predictors = str_remove_all(predictors,pattern = '_s|_w'))

proc_imp_model_thin%>%
  ggplot()+
  geom_bar(aes(predictors,m_imp,fill =season),stat = 'identity',position = 'dodge')+
  labs(y =  'index of importance')+
  coord_flip()+
  facet_grid(model~species)+theme_bw()
proc_imp_model_noThin%>%
  ggplot()+
  geom_bar(aes(x = predictors,y = m_imp,fill=season),stat = 'identity', position = 'dodge')+
  labs(y = 'index of importance')+
  coord_flip()+
  facet_grid(model~species)+
  theme_bw()

# --------  model performance ---------------------------------------------
proc_eval_model_thin <- eval_all_models_Thin%>%as_tibble()%>%
mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  unnest(cols = evals)
proc_eval_model_noThin <- eval_all_models_noThin%>%as_tibble()%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  unnest(cols = evals)

proc_eval_model_thin%>%
  filter(predictors=='AUC'|predictors=='sensitivity'|predictors=='specificity'|predictors=='Kappa')%>%
  ggplot()+
  geom_bar(aes(predictors,m_imp,fill =season),stat = 'identity',position = 'dodge')+
  labs(y = 'performance measure',x ='performance metrics')+
  coord_flip()+
  facet_grid(model~species)+theme_bw()

proc_eval_model_noThin%>%
  filter(predictors=='AUC'|predictors=='sensitivity'|predictors=='specificity'|predictors=='Kappa')%>%
  ggplot()+
  geom_bar(aes(predictors,m_imp,fill =season),stat = 'identity',position = 'dodge')+
  labs(y = 'performance measure',x ='performance metrics')+
  coord_flip()+
  facet_grid(model~species)+theme_bw()



# -----  processing the ensemble performance ------------------------------
proc_ens_imp_thin_sel <- proc_ens_perf_thin%>%
  filter(metrics=='AUC'|metrics=='sensitivity'|metrics=='specificity'|metrics=='Kappa')%>%
  dplyr::select(metrics,species,season,fp)
proc_ens_imp_no_thin_sel <- proc_ens_perf_noThin%>%
  filter(metrics=='AUC'|metrics=='sensitivity'|metrics=='specificity'|metrics=='Kappa')%>%
  dplyr::select(metrics,species,season,fp)

# ------- write model and ensemble performance ----------------------------
save(list = c('proc_eval_model_thin','proc_eval_model_noThin','proc_ens_perf_thin',
              'proc_ens_perf_noThin','proc_ens_imp_thin','proc_ens_imp_noThin',
              'proc_imp_model_noThin','proc_imp_model_thin','proc_both_ens_thinUnThin',
              'proc_ens_imp_no_thin_sel','proc_ens_imp_thin_sel'),
     file = 'analysis/out_tables/ensemble_model_perf_importance.RData')


# -------- combine model and ensemble performance -------------------------
load('analysis/out_tables/ensemble_model_perf_importance.RData')


proc_eval_model_thin%>%
  filter(predictors=='AUC'|predictors=='sensitivity'|predictors=='specificity'|predictors=='Kappa')%>%
  dplyr::rename(metrics = predictors)%>%
  ggplot()+
  geom_bar(aes(metrics,m_imp,fill =season),stat = 'identity',position = 'dodge')+
  geom_hline(data = proc_ens_imp_thin_sel%>%filter(metrics=='AUC'),aes(color=season,yintercept = fp),linetype=2)+
  labs(y = 'performance measure',x ='performance metrics')+
  coord_flip()+
  facet_grid(model~species)+
  theme_bw()
  
proc_eval_model_noThin%>%
  filter(predictors=='AUC'|predictors=='sensitivity'|predictors=='specificity'|predictors=='Kappa')%>%
  dplyr::rename(metrics = predictors)%>%
  ggplot()+
  geom_bar(aes(metrics,m_imp,fill =season),stat = 'identity',position = 'dodge')+
  geom_hline(data = proc_ens_imp_no_thin_sel%>%filter(metrics=='AUC'),aes(color=season,yintercept = fp),linetype=2)+
  labs(y = 'performance measure',x ='performance metrics')+
  coord_flip()+
  facet_grid(model~species)+
  theme_bw()

# ------ process  and present model partial effects -----------------------
resp_all_models_noThin = tibble(species_season = species_season_no_thin,
          resp = lapply(all_resp_curve_noThin,function(x) purrr::map_df(.x = x,.f = rbind,.id = 'Model')))
resp_all_models_thin = tibble(species_season = species_season_thin,
          resp = lapply(all_resp_curve_thin, function(x) purrr::map_df(.x = x,.f = rbind,.id = 'Model')))

# ------- process raw partial effect --------------------------------------
proc_resp_model_thin <- resp_all_models_thin%>%unnest(cols = resp)%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  mutate(vars = str_remove_all(vars,pattern = '_s|_w'))

proc_resp_model_noThin <- resp_all_models_noThin%>%unnest(cols = resp)%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  mutate(vars = str_remove_all(vars,pattern = '_s|_w'))


proc_resp_model_noThin%>%filter(Model=='RF')%>%
  ggplot()+
  geom_line(aes(n_vals,m_pred,color=season))+
  facet_grid(species~vars,scale = 'free')+
  theme_bw()

proc_resp_model_noThin%>%filter(Model=='GBM')%>%
  ggplot()+
  geom_line(aes(n_vals,m_pred,color=season))+
  facet_grid(species~vars,scale = 'free')+
  theme_bw()

proc_resp_model_noThin%>%filter(Model=='CTA')%>%
  ggplot()+
  geom_line(aes(n_vals,m_pred,color=season))+
  facet_grid(species~vars,scale = 'free')+
  theme_bw()

proc_resp_model_thin%>%filter(Model=='RF')%>%
  ggplot()+
  geom_line(aes(n_vals,m_pred,color=season))+
  facet_grid(species~vars,scale = 'free')+
  theme_bw()

proc_resp_model_thin%>%filter(Model=='GBM')%>%
  ggplot()+
  geom_line(aes(n_vals,m_pred,color=season))+
  facet_grid(species~vars,scale = 'free')+
  theme_bw()


# ---------  save processed response curves to file -----------------------
save(list = c('proc_resp_model_noThin','proc_resp_model_thin'),
     file = 'analysis/out_tables/processed_response_curves.RData')

#  -----  prepre occurrence maps with basemap -----------------------------
world <- getMap(resolution = "high")
### convert the spatialPolygonsDataFrame to an sf object
world_sf <- st_as_sf(world)
al_shape <- world_sf%>%
  st_crop(xmin=0.05,xmax=80,ymin=-38.45,ymax=-16.05)%>%lwgeom::st_make_valid()

# --------- utility function plot  predicted distribution -----------------
plot_pred_dist = function(pred_dat,bkgd_plt,which_fill){
  ggplot2::ggplot()+
    geom_tile(data = pred_dat,aes(x,y,fill=.data[[which_fill]]))+
    scale_fill_gradientn(colours = fields::tim.colors(100),name = 'probability of\n occurrence')+
    geom_sf(data = bkgd_plt,fill='gray50')+
    labs(x = 'longitude',y='latitude')+
    facet_wrap(~Model)+
    theme_bw()
}

# -----------------  ensemble and model prediction ------------------------
proc_pred_dist_ens_thin <-  pred_dist_ens_thin%>%
  mutate(pred = purrr::map(pred_dist,.f = function(x =.x) rasterToPoints(x)%>%as_tibble()),
         uncer = purrr::map(uncertainty,.f = function(x =.x) rasterToPoints(x)%>%
          as_tibble()%>%dplyr::rename(uncertainty = uncertainty.map)%>%dplyr::select(-c(x,y))))%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  dplyr::select(species_season,species,season,pred,uncer)%>%
  unnest(cols = c(pred,uncer))

proc_pred_dist_ens_noThin <-  pred_dist_ens_noThin%>%
  mutate(pred = purrr::map(preds_dist,.f = function(x =.x) rasterToPoints(x)%>%as_tibble()),
         uncer = purrr::map(uncertainty,.f = function(x =.x) rasterToPoints(x)%>%
      as_tibble()%>%dplyr::rename(uncertainty = uncertainty.map)%>%dplyr::select(-c(x,y))))%>%
  mutate(species_season =  species_season%>%str_remove_all(pattern = ".Ensemble.SDM"))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  dplyr::select(species_season,species,season,pred,uncer)%>%
  unnest(cols = c(pred,uncer))
  
# ----------  prepare ensemble distribution maps ---------------------------
proc_pred_dist_ens_thin_map <- proc_pred_dist_ens_thin%>%
  nest(data=c(x,y,Probability,uncertainty))%>%
  mutate(dist_maps = purrr::map(.x = data,
    .f = ~plot_pred_dist(pred_dat = .x,bkgd_plt = al_shape,which_fill ='Probability')))

proc_pred_dist_ens_thin%>%
  ggplot()+
  geom_tile(aes(x,y,fill=Probability))+
  scale_fill_gradientn(colours = fields::tim.colors(100),name='probability of\n occurrence')+
  geom_sf(data = al_shape,fill='gray50')+
  labs(x = 'longitude',y='latitude')+
  facet_grid(species~season)+
  theme_bw()

proc_pred_dist_ens_noThin%>%
  ggplot()+
  geom_tile(aes(x,y,fill=Probability))+
  scale_fill_gradientn(colours = fields::tim.colors(100),name='probability of\n occurrence')+
  geom_sf(data = al_shape,fill='gray50')+
  labs(x = 'longitude',y='latitude')+
  facet_grid(species~season)+
  theme_bw()

proc_pred_dist_ens_thin%>%
  ggplot()+
  geom_tile(aes(x,y,fill=uncertainty))+
  scale_fill_gradientn(colours = fields::tim.colors(100),name='uncertainty')+
  geom_sf(data = al_shape,fill='gray50')+
  labs(x = 'longitude',y='latitude')+
  facet_grid(species~season)+
  theme_bw()

proc_pred_dist_ens_noThin%>%
  ggplot()+
  geom_tile(aes(x,y,fill=uncertainty))+
  scale_fill_gradientn(colours = fields::tim.colors(100),name='uncertainty')+
  geom_sf(data = al_shape,fill='gray50')+
  labs(x = 'longitude',y='latitude')+
  facet_grid(species~season)+
  theme_bw()
# --------------------  model prediction  maps ----------------------------

#  -------  process model base predicted distirbution ---------------------
proc_pred_dist_model_thin <-  pred_dist_all_models_Thin%>%
  mutate(pred = purrr::map(pred_rasters,.f = function(x =.x) rasterToPoints(x)%>%
            as_tibble()%>%pivot_longer(cols = CTA:ANN,names_to = 'Model',values_to ='probs')))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  dplyr::select(species_season,species,season,pred)

proc_pred_dist_model_noThin <-  pred_dist_all_models_noThin%>%
  mutate(pred = purrr::map(pred_rasters,.f = function(x =.x) rasterToPoints(x)%>%
                             as_tibble()%>%pivot_longer(cols = CTA:ANN,names_to = 'Model',values_to ='probs')))%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)%>%
  dplyr::select(species_season,species,season,pred)


#  -------  plo distribution maps -----------------------------------------
proc_pred_dist_model_thin <- proc_pred_dist_model_thin%>%
  mutate(dist_maps  =purrr::map(.x = pred,
.f = ~plot_pred_dist(pred_dat = .x,bkgd_plt = al_shape,which_fill = 'probs')))

proc_pred_dist_model_noThin <- proc_pred_dist_model_noThin%>%
  mutate(dist_maps  =purrr::map(.x = pred,
  .f = ~plot_pred_dist(pred_dat = .x,bkgd_plt = al_shape,which_fill = 'probs')))


# ------- save ensemble and model prediction ------------------------------
save(list = c('proc_pred_dist_model_noThin','proc_pred_dist_model_thin','al_shape',
              'proc_pred_dist_ens_thin','proc_pred_dist_ens_noThin'),
     file = 'analysis/out_tables/model_ensemble_prediction.RData')

# ----------------------- processed all background data -------------------
# -------- rasrer file  used in distribution modelling --------------------
all_env_sum <- raster::stack('spatial/input_raster/env_all.grd')
# ------ select part of the data to be used for visualization -------------
tmp_sel_fun  =  function(dat){
  lapply(X = dat, function(x = X) {
    if(nrow(x)==0){
      x
      }
      else{
        x%>%dplyr::rename(lon = X, lat =Y,reps = model_index)%>%
         dplyr::select(c(reps,lon,lat,Presence,train))%>%as_tibble()
        }})
  
}


# --- subset the  occurrence data -----------------------------------------
sub_occ_thin <-  purrr::map(.x = all_occ_bkgnd_thin, 
                            .f = ~tmp_sel_fun(dat = .x)%>%purrr::map_df(.f = rbind,.id = 'Model'))
sub_occ_noThin <-  purrr::map(.x = all_occ_bkgnd_noThin, 
                            .f = ~tmp_sel_fun(dat = .x)%>%purrr::map_df(.f = rbind,.id = 'Model'))




# ------- prepare basemap -------------------------------------------------
basemap_al = prep_base_map(env_rast = all_env_sum,bathy_loc = 1,base_map = al_shape)

crs_all = crs(all_env_sum)

prep_convert_df =  function(dat,crs_use = crs_all){
  dat%>%filter(reps==1)%>%mutate(Presence = factor(Presence))%>%
     st_as_sf(coords = c('lon','lat'),crs = crs_all) 
 }

#  ----- the base map -----------------------------------------------------
occ_data_plot_noThin = tibble(species_season = species_season_no_thin,
                                sub_dat = sub_occ_noThin)
occ_data_plot_thin = tibble(species_season = species_season_no_thin,
                              sub_dat = sub_occ_thin)

# -----  the dist maps ----------------------------------------------------
occ_data_plot_thin <- occ_data_plot_thin%>%
  mutate(prep_dat = purrr::map(.x = sub_dat,.f = prep_convert_df))
occ_data_plot_thin = occ_data_plot_thin%>%
  mutate(bkgnd_plts = purrr::map(.x = prep_dat,
                                 .f = function(x = .x) basemap_al+
                                   geom_sf(data = x, aes(color=Presence),alpha=0.7)+facet_wrap(~Model)))
occ_data_plot_noThin = occ_data_plot_noThin%>%
  mutate(prep_dat =  purrr::map(.x = sub_dat,prep_convert_df),
         bkgnd_plts =  purrr::map(.x = prep_dat,
                                  .f = function(x = .x) basemap_al+
                                    geom_sf(data = x,aes(color = Presence),alpha=0.7)+facet_wrap(~Model)))

occ_data_plot_thin <- occ_data_plot_thin%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)
occ_data_plot_noThin <- occ_data_plot_noThin%>%
  separate(col = species_season,into = c('species','season'),sep = '_',remove = FALSE)

# ----------------- write processed background data -----------------------
save(list = c('occ_data_plot_noThin','occ_data_plot_thin'),
     file = 'analysis/out_tables/processed_all_background_data.RData')