#### From Lagrangian to Eulerian perspectives
#### Script: Summarising results
#### Last modification: 04 04 2024
#### R version 4.2.3 (2023-03-15 ucrt) -- "Shortstop Beagle"

library(tidyverse)
library(sf)
library(terra)
library(virtualecologist)
library(ggspatial)
library(tidyterra)
library(ggthemes)
library(viridis)
library(mgcv)
library(Distance)
library(ggpubr)

workdir <- "your working directory"
layout <- theme_bw() + 
  theme(plot.title.position = "panel",
        strip.text = element_text(hjust = 0, face = "italic", size = 12),
        strip.background = element_blank(),
        panel.grid = element_blank(),
        axis.title = element_blank(),
        axis.text = element_blank())

### 0 - Set up resource layers -------------------------------------------------
## Create virtual space
grid <- create_grid()

## Generate env conditions
env <- generate_env_layer(norm = TRUE, return_rasters = TRUE, grid = grid,
                          seed = 2)
plot(env$rasters)

## Generate resource layer for central place forager
rsce_1 <- generate_resource_layer(env_layers = env$rasters,
                                  beta = c(2, -1.5))
plot(rsce_1$rasters)

## Generate resource layer for free ranger
rsce_2 <- generate_resource_layer(env_layers = env$rasters,
                                  beta = c(-2, 0.5))
plot(rsce_2$rasters)

## Compute distance from colonies 
coord_col <- data.frame(Lon = c(30,60),
                        Lat = c(60,30),
                        Col = c("1", "2"))

dist_col <- raster::distanceFromPoints(raster::raster(rsce_1$rasters), xy = coord_col[,1:2])
names(dist_col) <- "dist_col"

## Plot
env_df <- list(env$dataframe |> rename(Resource1 = sim1, Resource2 = sim2), 
               rsce_1$dataframe |> rename(Suitability_CPF = suitability),
               rsce_2$dataframe |> rename(Suitability_FR = suitability))
env_long <- env_df %>% reduce(full_join, by =  c("x", "y")) |> 
  pivot_longer(cols = all_of(c("Resource1", "Resource2", "Suitability_CPF", "Suitability_FR"))) |> 
  mutate(name = case_when(name == "Resource1" ~ "Resource 1",
                          name == "Resource2" ~ "Resource 2",
                          name == "Suitability_CPF" ~ "Suitability (CPF)",
                          name == "Suitability_FR" ~ "Suitability (FR)"))



### 1 - Open true distribution -------------------------------------------------
# true distributions
CPF_true <- terra::rast(paste0(workdir, "/Simulated_tracks_CPF_KDE_totpop.grd"))
CPF_true_contour <- st_read(paste0(workdir, "/Simulated_tracks_CPF_KDE_totpop_poly.shp"))

FR_true <- terra::rast(paste0(workdir, "/Simulated_tracks_FR_KDE_totpop.grd"))
FR_true_contour <- st_read(paste0(workdir, "/Simulated_tracks_FR_KDE_totpop_poly.shp"))

st_crs(CPF_true_contour) <- virtualecologist::custom_proj()
st_crs(FR_true_contour) <- virtualecologist::custom_proj()


### 2 - Open prediction maps ---------------------------------------------------
cpf_rsf <- rast(paste0(workdir, "/CPF_RSF_abundance_map.grd"))
fr_rsf <- rast(paste0(workdir, "/FR_RSF_abundance_map.grd"))

cpf_dsm_boat <- rast(paste0(workdir, "/CPF_boat_abundance_map.grd"))
cpf_dsm_plane <- rast(paste0(workdir, "/CPF_plane_abundance_map.grd"))
fr_dsm_boat <- rast(paste0(workdir, "/FR_boat_abundance_map.grd"))
fr_dsm_plane <- rast(paste0(workdir, "/FR_plane_abundance_map.grd"))


### 3 - Compute 80% contour of prediction maps ---------------------------------

cpf_rsf_contour <- raster::rasterToContour(raster::cut(raster::raster(cpf_rsf), 
                               breaks = c(0, 
                                          quantile(raster::raster(cpf_rsf), probs = 0.5),
                                          quantile(raster::raster(cpf_rsf), probs = 0.8), 
                                          quantile(raster::raster(cpf_rsf), probs = 0.9), 
                                          quantile(raster::raster(cpf_rsf), probs = 1))), levels = c(1,2,3,4)) %>%
  st_as_sf() |> mutate(Class = case_when(level == "2" ~ "[50,80%)", level == "3" ~ "[80,90%)", level == "4" ~ "[90,100%)"))
fr_rsf_contour <- raster::rasterToContour(raster::cut(raster::raster(fr_rsf), 
                                             breaks = c(0, 
                                                        quantile(raster::raster(fr_rsf), probs = 0.5),
                                                        quantile(raster::raster(fr_rsf), probs = 0.8), 
                                                        quantile(raster::raster(fr_rsf), probs = 0.9), 
                                                        quantile(raster::raster(fr_rsf), probs = 1))), levels = c(1,2,3,4)) %>%
  st_as_sf() |> mutate(Class = case_when(level == "2" ~ "[50,80%)", level == "3" ~ "[80,90%)", level == "4" ~ "[90,100%)"))


cpf_dsm_boat_contour <- raster::rasterToContour(raster::cut(raster::raster(cpf_dsm_boat), 
                                              breaks = c(0, 
                                                         quantile(raster::raster(cpf_dsm_boat), probs = 0.5),
                                                         quantile(raster::raster(cpf_dsm_boat), probs = 0.8), 
                                                         quantile(raster::raster(cpf_dsm_boat), probs = 0.9), 
                                                         quantile(raster::raster(cpf_dsm_boat), probs = 1))), levels = c(1,2,3,4)) %>%
  st_as_sf() |> mutate(Class = case_when(level == "2" ~ "[50,80%)", level == "3" ~ "[80,90%)", level == "4" ~ "[90,100%)"))
cpf_dsm_plane_contour <- raster::rasterToContour(raster::cut(raster::raster(cpf_dsm_plane), 
                                                            breaks = c(0, 
                                                                       quantile(raster::raster(cpf_dsm_plane), probs = 0.5),
                                                                       quantile(raster::raster(cpf_dsm_plane), probs = 0.8), 
                                                                       quantile(raster::raster(cpf_dsm_plane), probs = 0.9), 
                                                                       quantile(raster::raster(cpf_dsm_plane), probs = 1))), levels = c(1,2,3,4)) %>%
  st_as_sf() |> mutate(Class = case_when(level == "2" ~ "[50,80%)", level == "3" ~ "[80,90%)", level == "4" ~ "[90,100%)"))


fr_dsm_boat_contour <- raster::rasterToContour(raster::cut(raster::raster(fr_dsm_boat), 
                                                            breaks = c(0, 
                                                                       quantile(raster::raster(fr_dsm_boat), probs = 0.5),
                                                                       quantile(raster::raster(fr_dsm_boat), probs = 0.8), 
                                                                       quantile(raster::raster(fr_dsm_boat), probs = 0.9), 
                                                                       quantile(raster::raster(fr_dsm_boat), probs = 1))), levels = c(1,2,3,4)) %>%
  st_as_sf() |> mutate(Class = case_when(level == "2" ~ "[50,80%)", level == "3" ~ "[80,90%)", level == "4" ~ "[90,100%)"))
fr_dsm_plane_contour <- raster::rasterToContour(raster::cut(raster::raster(fr_dsm_plane), 
                                                             breaks = c(0, 
                                                                        quantile(raster::raster(fr_dsm_plane), probs = 0.5),
                                                                        quantile(raster::raster(fr_dsm_plane), probs = 0.8), 
                                                                        quantile(raster::raster(fr_dsm_plane), probs = 0.9), 
                                                                        quantile(raster::raster(fr_dsm_plane), probs = 1))), levels = c(1,2,3,4)) %>%
  st_as_sf() |> mutate(Class = case_when(level == "2" ~ "[50,80%)", level == "3" ~ "[80,90%)", level == "4" ~ "[90,100%)"))



contours_cpf <- rbind(cpf_rsf_contour |> mutate(Species = "Central-place forager", name = "RSF"),
                  cpf_dsm_boat_contour |> mutate(Species = "Central-place forager", name = "DSM - boat"),
                  cpf_dsm_plane_contour |> mutate(Species = "Central-place forager", name = "DSM - plane"))
contours_fr <- rbind(fr_rsf_contour |> mutate(Species = "Free ranger", name = "RSF"),
                     fr_dsm_boat_contour |> mutate(Species = "Free ranger", name = "DSM - boat"),
                     fr_dsm_plane_contour |> mutate(Species = "Free ranger", name = "DSM - plane"))

st_crs(contours_cpf) <- virtualecologist::custom_proj()
st_crs(contours_fr) <- virtualecologist::custom_proj()



### 4 - Correlation ------------------------------------------------------------
# Pearson's coefficients - function
correlations <- function(predictors, effort, method, plot = T, ...) {
  corr <- Hmisc::rcorr(as.matrix(effort[!is.na(effort[,predictors]), predictors]), type = "pearson")
  # symnum(corr, abbr.colnames=FALSE)
  if(plot == T){
    corrplot::corrplot(corr$r, type="upper", order="original", tl.col="black", tl.srt=45, method = method, diag = F,
                     addCoefasPercent = F, tl.cex = 0.7, is.corr = F, number.cex = .5,...)}
  if(plot == F){ return(corr$r) }
}

# CPF 
cpf_corr <- cbind(as.data.frame(cpf_dsm_boat), as.data.frame(cpf_dsm_plane), 
                 as.data.frame(cpf_rsf),
                 as.data.frame(rsce_1$rasters), as.data.frame(CPF_true))
names(cpf_corr) <- c("DSM - boat", "DSM - plane", "RSF", "Suitability", "True ab.")

correlations(effort = cpf_corr, predictors = c("DSM - boat", "DSM - plane", "RSF", "Suitability", "True ab."),
             method = "number", plot = T)

# FR
fr_corr <- cbind(as.data.frame(fr_dsm_boat), as.data.frame(fr_dsm_plane), 
                  as.data.frame(fr_rsf),
                  as.data.frame(rsce_2$rasters), as.data.frame(FR_true))
names(fr_corr) <- c("DSM - boat", "DSM - plane", "RSF", "Suitability", "True ab.")

correlations(effort = fr_corr, predictors = c("DSM - boat", "DSM - plane", "RSF", "Suitability", "True ab."),
             method = "number", plot = T) # "pie"

# plot both
png("Correlation_coef.png", height = 7.2, width = 14, res = 350, units = "cm")
par(mfrow = c(1,2))
correlations(effort = cpf_corr, predictors = c("RSF", "DSM - boat", "DSM - plane", "Suitability", "True ab."),
             method = "color", plot = T, cl.pos = "n", addCoef.col = "black", col.lim = c(0,1))
mtext("Central-place forager", side = 3, line = 3.1, adj = 0)
correlations(effort = fr_corr, predictors = c("RSF", "DSM - boat", "DSM - plane", "Suitability", "True ab."),
             method = "color", plot = T, cl.pos = "n", addCoef.col = "black", col.lim = c(0,1))
mtext("Free ranger", side = 3, line = 3.1, adj = 0)
dev.off()



### 5 - Plot predictions maps --------------------------------------------------
cpf_df <- cbind(as.data.frame(CPF_true, xy = T)[,1:2], cpf_corr) |> 
  mutate(Species = "Central-place forager") |> 
  pivot_longer(cols = c("RSF", "DSM - boat", "DSM - plane", "True ab.")) |> 
  mutate(name = case_when(name == "True ab." ~ "True abundance", TRUE ~ name))
cpf_df$name <- factor(cpf_df$name, levels = c("True abundance", "RSF", "DSM - plane", "DSM - boat"))

fr_df <- cbind(as.data.frame(FR_true, xy = T)[,1:2], fr_corr) |> 
  mutate(Species = "Free ranger") |> 
  pivot_longer(cols = c("RSF", "DSM - boat", "DSM - plane", "True ab.")) |> 
  mutate(name = case_when(name == "True ab." ~ "True abundance", TRUE ~ name))
fr_df$name <- factor(fr_df$name, levels = c("True abundance", "RSF", "DSM - plane", "DSM - boat"))


rsce <- ggplot() +  
  geom_sf(data = FR_true_contour)+
  geom_tile(data = env_long, aes(x = x, y = y, fill = value))+
  facet_grid(name~.)+
  scale_fill_viridis(option = "H", name = "Value") +
  layout 

pred_cpf <- ggplot() + 
  geom_tile(data = cpf_df, aes(x = x, y = y, fill = value)) + 
  facet_grid(name~.) +
  scale_fill_viridis(option = "H", name = "N. ind") +
  geom_point(data = coord_col, aes(x = Lon, y = Lat), color = "pink") +
  geom_sf(data = CPF_true_contour, color = "pink", linetype = "dotted") + 
  labs(title = "B - Central-place forager") +
  layout

pred_fr <- ggplot() + 
  geom_tile(data = fr_df, aes(x = x, y = y, fill = value)) + 
  facet_grid(name~.) +
  scale_fill_viridis(option = "H", name = "N. ind") +
  geom_sf(data = FR_true_contour, color = "pink", linetype = "dotted") + 
  labs(title = "C - Free ranger") +
  layout

ggarrange(rsce + labs(title = "A - Resource layers"), 
          pred_cpf, 
          pred_fr, common.legend = F, ncol = 3, align = "hv") +
  bgcolor("white")+border("white")
ggsave(paste0(workdir, "/Final_map.png"), height = 8, width = 10, dpi = 450)




