#### From Lagrangian to Eulerian perspectives
#### Script: Analysing telemetry data
#### Last modification: 16 02 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"
setwd(workdir)

#### ---------------------------------- ####
#### --- Create virtual environment --- ####
#### ---------------------------------- ####

## 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"

## create prediction table
env_stack <- raster::stack(env$rasters)
env_stack <- raster::addLayer(env_stack, dist_col)
env_stack <- raster::addLayer(env_stack, raster::stack(rsce_1$rasters), raster::stack(rsce_2$rasters))

env_df <- as.data.frame(env_stack, xy = T) |> 
  rename(Resource1 = sim1, Resource2 = sim2,
         Suitability_CPF = suitability.1, Suitability_FR = suitability.2) |> 
  mutate(area = 0.5*0.5) 
 
 
ggplot(env_df)+ geom_tile(aes(x = x, y = y, fill = Resource1)) + scale_fill_viridis(option = "H")
ggplot(env_df)+ geom_tile(aes(x = x, y = y, fill = Resource2)) + scale_fill_viridis(option = "H")
ggplot(env_df)+ geom_tile(aes(x = x, y = y, fill = Suitability_CPF)) + scale_fill_viridis(option = "H")
ggplot(env_df)+ geom_tile(aes(x = x, y = y, fill = Suitability_FR)) + scale_fill_viridis(option = "H")
ggplot(env_df)+ geom_tile(aes(x = x, y = y, fill = dist_col)) + scale_fill_viridis(option = "H")


## True maps -------------------------------------------------------------------
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"))



#### ------------------------------- ####
#### ---- Open and prepare data ---- ####
#### ------------------------------- ####

### 1 - Open data --------------------------------------------------------------
## CPF
cpf <- read.table(paste0(workdir, "/Simulated_tracks_CPF_sampled_12.csv"), h = T, dec = ".", sep = ";") |> 
  mutate(Day = as_date(Day),
         ymd_hms = ymd_hms(ymd_hms),
         Used = TRUE,
         ID = BirdId)

table(cpf$Colony)
cpf |> group_by(Colony2) |> summarize(n = length(unique(ID)))
unique_col <- c("Colony 1", "Colony 2")

## FR
fr <- read.table(paste0(workdir, "/Simulated_tracks_FR_sampled.csv"), h = T, dec = ".", sep = ";") |> 
  mutate(Day = as_date(Day),
         ymd_hms = ymd_hms(ymd_hms),
         Used = TRUE,
         ID = ind_id)

fr |> summarize(n = length(unique(ID)))


### 2 - Visualise the data -----------------------------------------------------
source("functions_track_summary.R")

## A - CPF
ggplot(cpf) +
  geom_spatraster(data = rsce_1$rasters, alpha = 0.5) + 
  geom_point(aes(x = Lon, y = Lat, color = activity)) +
  scale_fill_viridis(option = "H") +
  facet_grid(Colony~Day)

cpf_sum <- tracksummary_func(cpf, n_reloc_per_hour = 4)
cpf_sum$plot_distrib ; ggsave("Simulated_tracks_CPF_sampled_distribution.png", height = 8, width = 5, dpi = 350)
cpf_sum$plot_proptrack ; ggsave("Simulated_tracks_CPF_sampled_activities.png", height = 3, width = 5, dpi = 350)



## B - FR
ggplot(fr) +
  geom_spatraster(data = rsce_2$rasters, alpha = 0.5) + 
  geom_point(aes(x = Lon, y = Lat, color = activity)) +
  scale_fill_viridis(option = "H") +
  facet_grid(.~Day)

fr_sum <- tracksummary_func(trajectories = fr, n_reloc_per_hour = 4, type = "FR")
fr_sum$plot_distrib ; ggsave("Simulated_tracks_FR_sampled_distribution.png", height = 5, width = 5, dpi = 350)
fr_sum$plot_proptrack ; ggsave("Simulated_tracks_FR_sampled_activities.png", height = 3, width = 4, dpi = 350)




### 3 - Format the data --------------------------------------------------------
source("functions_rsf_analysis.R")
cpf_data <- format_track_data(traj = cpf, coord_col = coord_col, 
                             col_sizes = c(8000, 12000), type = "CPF")
fr_data <- format_track_data(traj = fr, coord_col = coord_col, 
                             col_sizes = NULL, type = "FR")



### 4 - Run the RSF ------------------------------------------------------------
source("contrib_function_202303.R")

cpf_rsf <- rsf_fitting_gam(df = cpf_data$final_table, buffers = cpf_data$buffers,
                       Species = "CPF", col_sizes = c(8000, 10000))
fr_rsf <- rsf_fitting_gam(df = fr_data, buffers = NULL,
                       Species = "FR", col_sizes = 8000)


### 5 - Look at results --------------------------------------------------------
cpf_rsf$RSF_results
fr_rsf$RSF_results



### 6 - Plot prediction maps ---------------------------------------------------

# plot
par(mfrow = c(1,2), oma = c(0,0,2,2), mar = c(2.5,2.5,1.5,5))
plot(cpf_rsf$Abundance_map[[1]], col = viridis(20, option = "H"))
lines(CPF_true_contour, col = "pink") ; points(coord_col[, c("Lon", "Lat")], col = "white", pch = 20)
plot(fr_rsf$Abundance_map[[1]], col = viridis(20, option = "H"))
lines(FR_true_contour, col = "pink")




### 7 - Functional responses ---------------------------------------------------
png("CPF_model_curves_RSF.png", height = 8, width = 12, res = 450, units = "cm")
par(mfrow = c(2,3), mar = c(2,2,1.1,0.5), oma = c(0.1,1.5,0.8,0.5), mgp = c(2.5,0.5,0), tck = -0.04, cex = 0.7)
plot(cpf_rsf$Model, shade = T, shade.col = "#147608ff", scale = 0, select = 1, xlim = c(0,1), rug = T)
mtext("log(ind)", side = 2, line = 2, cex = 0.8)
mtext("Resource 1", side = 3, line = 0.3, cex = 0.8, adj = 0)
plot(cpf_rsf$Model, shade = T, shade.col = "#147608ff", scale = 0, select = 2, xlim = c(0,1), rug = T)
mtext("Resource 2", side = 3, line = 0.3, cex = 0.8, adj = 0)
plot.new()

plot.new()
plot.new()
plot.new()
dev.off()

png("FR_model_curves_RSF.png", height = 8, width = 12, res = 450, units = "cm")
par(mfrow = c(2,3), mar = c(2,2,1.1,0.5), oma = c(0.1,1.5,0.8,0.5), mgp = c(2.5,0.5,0), tck = -0.04, cex = 0.7)
plot(fr_rsf$Model, shade = T, shade.col = "#8c3eb2ff", scale = 0, select = 1, xlim = c(0,1), rug = T)
mtext("log(ind)", side = 2, line = 2, cex = 0.8)
mtext("Resource 1", side = 3, line = 0.3, cex = 0.8, adj = 0)
plot(fr_rsf$Model, shade = T, shade.col = "#8c3eb2ff", scale = 0, select = 2, xlim = c(0,1), rug = T)
mtext("Resource 2", side = 3, line = 0.3, cex = 0.8, adj = 0)
plot.new()

plot.new()
plot.new()
plot.new()
dev.off()



### 8 - Plot the data ----------------------------------------------------------
ggplot() + 
  geom_point(data = cpf_data$final_table %>% subset(Used == "FALSE") %>% mutate(Colony = paste("Colony", Colony, sep = " ")),
             aes(x = Lon, y = Lat), color = "grey") +
  geom_point(data = cpf_data$final_table %>% subset(Used == "TRUE") %>% mutate(Colony = paste("Colony", Colony, sep = " ")),
             aes(x = Lon, y = Lat), color = "black") +
  geom_point(data = coord_col %>% mutate(Colony = paste("Colony", Col, sep = " ")), 
             aes(x = Lon, y = Lat), color = "pink", size = 2) +
  facet_wrap("Colony") +
  labs(x = "x", y = "y") +
  theme_bw() + 
  xlim(c(0,90)) + ylim(c(0,90))+
  theme(panel.grid = element_blank(), 
        strip.background = element_blank(),
        strip.text = element_text(face = 2, size = 18, hjust = 0),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 16))
ggsave(paste0(workdir, "/RSF_data_CPF.png"), height = 6.5, width = 12, dpi = 350)


ggplot() + 
  geom_point(data = fr_data %>% subset(Used == "FALSE"),
             aes(x = Lon, y = Lat), color = "grey") +
  geom_point(data = fr_data %>% subset(Used == "TRUE"),
             aes(x = Lon, y = Lat), color = "black") +
  theme_bw() + 
  labs(x = "x", y = "y") +
  xlim(c(0,90)) + ylim(c(0,90))+
  theme(panel.grid = element_blank(), 
        strip.background = element_blank(),
        strip.text = element_text(face = 2, size = 18),
        axis.title = element_text(size = 16),
        axis.text = element_text(size = 14),
        legend.text = element_text(size = 14),
        legend.title = element_text(size = 16))
ggsave(paste0(workdir, "/RSF_data_FR.png"), height = 6, width = 6, dpi = 350)







save.image(paste0(workdir, "/RSF_analysis.RData"))
