#### From Lagrangian to Eulerian perspectives
#### Script: Simulate species movements
#### Last modification: 13 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)


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



#### ---------------------------------- ####
#### --- Generate species movements --- ####
#### ------ Free ranging species ------ ####
#### ---------------------------------- ####

## set up
# consider the species move at 30-50 km/h when travelling (so 0.5-0.8 km/min) [starting bout == travelling bout]
# at 10-20 km/h when foraging (so 0.16-0.33 km/min)


## Create all possible starting locations, bearings and steps
possible_starting_positions <- expand.grid(Lon = seq(5, 85, by = 0.1), Lat = seq(5, 85, by = 0.1))
possible_start_bearing <- randomdir(n = 1000, mu = 0, kappa = 90)
possible_start_step <- randomdist(1000, shape = 0.7, rate = 3)

## Generate movements for all individuals
traj_FR <- NULL
for(i in cli::cli_progress_along(1:8000, "Processing")){
  day1 <- simulate_trajectory_FR(initial_position = possible_starting_positions[ sample(nrow(possible_starting_positions), 1), 
                                                                                 c("Lon", "Lat")],
                                resource_layer = rsce_2$rasters,
                                starting_hour = ymd_hms("2022-08-02 06:00:00"),
                                starting_bearing = c(sample(possible_start_bearing, 1), 10), 
                                starting_step = c(sample(possible_start_step, 1), 3),
                                travel_bearing = c(0, 20), 
                                travel_step = c(0.7, 3), 
                                foraging_bearing = c(0, 0.5), 
                                foraging_step = c(0.2, 3),
                                minx = 0, maxx = 90, 
                                miny = 0, maxy = 90,
                                step_duration = 1,
                                activity_threshold = 0.7,
                                max_duration = 720, # 12h * 60 min
                                residence_time = 15,
                                residence_radius = 2)
  day2 <- simulate_trajectory_FR(initial_position = day1[ nrow(day1) , c("Lon", "Lat")],
                                 resource_layer = rsce_2$rasters,
                                 starting_hour = ymd_hms("2022-08-03 06:00:00"),
                                 starting_bearing = c(sample(possible_start_bearing, 1), 10), 
                                 starting_step = c(sample(possible_start_step, 1), 3),
                                 travel_bearing = c(0, 20), 
                                 travel_step = c(0.7, 3), 
                                 foraging_bearing = c(0, 0.5), 
                                 foraging_step = c(0.2, 3),
                                 minx = 0, maxx = 90, 
                                 miny = 0, maxy = 90,
                                 step_duration = 1,
                                 activity_threshold = 0.7,
                                 max_duration = 720, # 12h * 60 min
                                 residence_time = 15,
                                 residence_radius = 2)
  res <- rbind(day1, day2)
  res$Day <- as_date(res$ymd_hms)
  res$ind_id <- i
  
  traj_FR <- rbind(traj_FR, res)
}

save.image("Simulate_species.RData")

ggplot(traj_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_wrap("Day")




#### -------------------------------------- ####
#### ----- Generate species movements ----- ####
#### --- Central place foraging species --- ####
#### -------------------------------------- ####

## set up
# consider the species move at 50 km/h when departing colony (so 0.83 km/min)
# at 30-50 km/h when travelling (so 0.5-0.8 km/min)
# at 10-20 km/h when foraging (so 0.16-0.33 km/min)


### 1 - Colony 1 ---------------------------------------------------------------

Col1_tracks <-  NULL
for(i in cli::cli_progress_along(1:8000, "Processing Colony 1")){
  # day 1
  daylength1 <- insol::daylength(long = 25, lat = 25, jd = insol::JD(as.POSIXct("2022-08-02", tz = "UTC")), tmz = 0)
  sunrise1 <- format(as.POSIXct(daylength1[,1]*3600, origin = as.POSIXct("2022-08-02", tz = "UTC"), "%H:%M", tz = "UTC"))
  day1 <- suppressMessages(
    simulate_trajectory_CPF(
     initial_position = coord_col[1,1:2], 
     resource_layer = rsce_1$rasters, 
     starting_hour = sample(seq( ymd_hms(sunrise1),  # depart within 2 hours after sunrise
                                 ymd_hms(sunrise1) + hours(2), 
                                 by = "15 mins"), 1, replace = T),
     starting_bearing = c(270,10), 
     starting_step = c(0.9, 3),
     travel_bearing = c(0, 20), 
     travel_step = c(0.7, 3), 
     foraging_bearing = c(0, 0.5), 
     foraging_step = c(0.2, 3),
     activity_threshold = 0.6,
     step_duration = 1,
     max_dist = 40, # homing bout at max 40 km from col
     max_duration = 720, # homing bout after 12h
     residence_time = 15,
     residence_radius = 1) 
  )
  day1$Day <- as_date(day1$ymd_hms)
  day1$ind_id <- i
  day1$trip_id <- paste(day1$ind_id, "1", sep = "_")
  
  # day 2
  daylength2 <- insol::daylength(long = 25, lat = 25, jd = insol::JD(as.POSIXct("2022-08-03", tz = "UTC")), tmz = 0)
  sunrise2 <- format(as.POSIXct(daylength2[,1]*3600, origin = as.POSIXct("2022-08-03", tz = "UTC"), "%H:%M", tz = "UTC"))
  day2 <- suppressMessages(
    simulate_trajectory_CPF(
      initial_position = coord_col[1,1:2], 
      resource_layer = rsce_1$rasters, 
      starting_hour = sample(seq( ymd_hms(sunrise2),  # depart within 2 hours after sunrise
                                  ymd_hms(sunrise2) + hours(2), 
                                  by = "15 mins"), 1, replace = T),
      starting_bearing = c(270,10), 
      starting_step = c(0.9, 3),
      travel_bearing = c(0, 20), 
      travel_step = c(0.7, 3), 
      foraging_bearing = c(0, 0.5), 
      foraging_step = c(0.2, 3),
      activity_threshold = 0.6,
      step_duration = 1,
      max_dist = 40, 
      max_duration = 720, 
      residence_time = 15,
      residence_radius = 1) 
  )
  day2$Day <- as_date(day2$ymd_hms)
  day2$ind_id <- i
  day2$trip_id <- paste(day2$ind_id, "2", sep = "_")
  
  # merge
  Col1_tracks <- rbind(Col1_tracks, day1, day2)
} 

save.image("Simulate_species.RData")

ggplot(Col1_tracks) +
  geom_spatraster(data = rsce_1$rasters, alpha = 0.5) + 
  geom_point(aes(x = Lon, y = Lat, color = activity)) +
  scale_fill_viridis(option = "H", limits = c(0.6,1)) +
  facet_wrap("Day")




### 2 - Colony 2 ---------------------------------------------------------------

Col2_tracks <-  NULL
for(i in cli::cli_progress_along(1:10000, "Processing Colony 2")){
  # day 1
  daylength1 <- insol::daylength(long = 25, lat = 25, jd = insol::JD(as.POSIXct("2022-08-02", tz = "UTC")), tmz = 0)
  sunrise1 <- format(as.POSIXct(daylength1[,1]*3600, origin = as.POSIXct("2022-08-02", tz = "UTC"), "%H:%M", tz = "UTC"))
  day1 <- suppressMessages(
    simulate_trajectory_CPF(
      initial_position = coord_col[2,1:2], 
      resource_layer = rsce_1$rasters, 
      starting_hour = sample(seq( ymd_hms(sunrise1),  # depart within 2 hours after sunrise
                                  ymd_hms(sunrise1) + hours(2), 
                                  by = "15 mins"), 1, replace = T),
      starting_bearing = c(90,10), 
      starting_step = c(0.9, 3),
      travel_bearing = c(0, 20), 
      travel_step = c(0.7, 3), 
      foraging_bearing = c(0, 0.5), 
      foraging_step = c(0.2, 3),
      activity_threshold = 0.6,
      step_duration = 1,
      max_dist = 40, # homing bout at max 40 km from col
      max_duration = 720, # homing bout after 12h
      residence_time = 15,
      residence_radius = 1) 
  )
  day1$Day <- as_date(day1$ymd_hms)
  day1$ind_id <- i
  day1$trip_id <- paste(day1$ind_id, "1", sep = "_")
  
  # day 2
  daylength2 <- insol::daylength(long = 25, lat = 25, jd = insol::JD(as.POSIXct("2022-08-03", tz = "UTC")), tmz = 0)
  sunrise2 <- format(as.POSIXct(daylength2[,1]*3600, origin = as.POSIXct("2022-08-03", tz = "UTC"), "%H:%M", tz = "UTC"))
  day2 <- suppressMessages(
    simulate_trajectory_CPF(
      initial_position = coord_col[2,1:2], 
      resource_layer = rsce_1$rasters, 
      starting_hour = sample(seq( ymd_hms(sunrise2),  # depart within 2 hours after sunrise
                                  ymd_hms(sunrise2) + hours(2), 
                                  by = "15 mins"), 1, replace = T),
      starting_bearing = c(90,10), 
      starting_step = c(0.9, 3),
      travel_bearing = c(0, 20), 
      travel_step = c(0.7, 3), 
      foraging_bearing = c(0, 0.5), 
      foraging_step = c(0.2, 3),
      activity_threshold = 0.6,
      step_duration = 1,
      max_dist = 40, 
      max_duration = 720, 
      residence_time = 15,
      residence_radius = 1) 
  )
  day2$Day <- as_date(day2$ymd_hms)
  day2$ind_id <- i
  day2$trip_id <- paste(day2$ind_id, "2", sep = "_")
  
  # merge
  Col2_tracks <- rbind(Col2_tracks, day1, day2)
} 
save.image("Simulate_species.RData")

ggplot(Col2_tracks) +
  geom_spatraster(data = rsce_1$rasters, alpha = 0.5) + 
  geom_point(aes(x = Lon, y = Lat, color = activity)) +
  scale_fill_viridis(option = "H", limits = c(0.6,1)) +
  facet_wrap("Day")





### 4 - Merge the tracks -------------------------------------------------------
hist(Col1_tracks$step)
hist(Col2_tracks$step)

traj_CPF <- rbind( Col1_tracks %>% mutate(Colony = 1),
                   Col2_tracks %>% mutate(Colony = 2)) %>%
  mutate(BirdId = paste0("Col", Colony, "_", trip_id),
         Colony2 = case_when(Colony == 1 ~ "Colony 1",
                             Colony == 2 ~ "Colony 2")
  ) 

# check the numbers
traj_CPF %>% group_by(Colony, Day) %>% summarize(n = length(unique(trip_id)))



#### -------------------------------------- ####
#### ----- Save true and sampled traj ----- ####
#### -------------------------------------- ####

### 1 - Save the true movement files -------------------------------------------
save(traj_CPF, file = "Simulated_tracks_CPF.rda")


### 2 - Sample the free ranger -------------------------------------------------
traj_FR <- traj_FR |> mutate(trip_id = paste0(ind_id, "_", Day))
save(traj_FR, file = "Simulated_tracks_FR.rda")


### 3 - Sample the central place forager ---------------------------------------
## sample 1% of individuals per colonies
spl_col1 <- sample(x = unique(traj_CPF$ind_id[which(traj_CPF$Colony == 1)]), 
                   size = 80, replace = F)
spl_col2 <- sample(x = unique(traj_CPF$ind_id[which(traj_CPF$Colony == 2)]), 
                   size = 120, replace = F)
sample_traj_CPF_12 <- traj_CPF %>% subset(ind_id %in% c(spl_col1, spl_col2))


## downsample temporal resolution : one fix per 15 minutes
sample_traj_CPF_12 <- sample_traj_CPF_12 |> mutate(mins = lubridate::minute(ymd_hms)) |> filter(mins %in% c(00, 15, 30, 45)) |> select(-mins)
sample_traj_CPF_12 |> group_by(Colony, Day, BirdId)|> tally() |> ggplot() + geom_boxplot(aes(y = n, x = Colony, group = Day)) # 49 reloc par ind

write.table(sample_traj_CPF_12, "Simulated_tracks_CPF_sampled_12.csv", row.names = F, dec = ".", sep = ";")

ggplot(sample_traj_CPF_12) +
  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)


### 4 - Sample the free ranger -------------------------------------------------
## sample 1% of individuals per colonies
spl_fr <- sample(x = unique(traj_FR$ind_id), 
                   size = 80, replace = F)
sample_traj_FR <- traj_FR %>% subset(ind_id %in% spl_fr)


## downsample temporal resolution : one fix per 15 minutes
sample_traj_FR <- sample_traj_FR |> mutate(mins = lubridate::minute(ymd_hms)) |> filter(mins %in% c(00, 15, 30, 45)) |> select(-mins)
sample_traj_FR |> group_by(Day, ind_id)|> tally() |> ggplot() + geom_boxplot(aes(y = n, group = Day)) # 49 reloc par ind

write.table(sample_traj_FR, "Simulated_tracks_FR_sampled.csv", row.names = F, dec = ".", sep = ";")


ggplot(sample_traj_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)

save.image("Simulate_species.RData")



#### ---------------------------- ####
#### ---- Visualise the data ---- ####
#### ---------------------------- ####
source("functions_track_summary.R")
library(ggpubr)
wd <- "your working directory where to save the plots"


#### 1 - Central place forager -------------------------------------------------
# check the typology of movements
cpf_sum <- tracksummary_func(traj_CPF,  n_reloc_per_hour = 60)
cpf_sum$plot_distrib ; ggsave("Simulated_tracks_CPF_distributions.png", height = 8, width = 5, dpi = 350)
cpf_sum$plot_nfix
cpf_sum$plot_proptrack ; ggsave("Simulated_tracks_CPF_activities.png", height = 3, width = 5, dpi = 350)  # changement dans les proportions d'activité
cpf_sum$plot_tripduration # mais pas dans les autres parmaètres de trajets

# check the distribution of locations and activities
cpf_d1_nfix <- nfixpercell_function(traj_CPF |> subset(Day == "2022-08-02"), type = "CPF", outdir = wd, outname = "Simulated_CPF_Day1", base_raster = dist_col)
cpf_d2_nfix <- nfixpercell_function(traj_CPF |> subset(Day == "2022-08-03"), type = "CPF", outdir = wd, outname = "Simulated_CPF_Day2", base_raster = dist_col)

extent_function(traj_CPF |> subset(Day == "2022-08-02"), type = "CPF", outdir = wd, outname = "Simulated_CPF_Day1", base_raster = dist_col)
extent_function(traj_CPF |> subset(Day == "2022-08-03"), type = "CPF", outdir = wd, outname = "Simulated_CPF_Day2", base_raster = dist_col)



#### 2 - Free ranger -----------------------------------------------------------
traj_FR$BirdId <- traj_FR$trip_id

# check the typology of movements
fr_sum <- tracksummary_func(traj_FR,  n_reloc_per_hour = 60, type = "FR")
fr_sum$plot_distrib ; ggsave("Simulated_tracks_FR_distributions.png", height = 5, width = 5, dpi = 350)
fr_sum$plot_nfix # all the same, by construction
fr_sum$plot_proptrack ; ggsave("Simulated_tracks_FR_activities.png", height = 3, width = 4, dpi = 350)  # changement dans les proportions d'activité
fr_sum$plot_tripduration # all the same, by construction

# check the distribution of locations and activities
fr_d1_nfix <- nfixpercell_function(trajectories = traj_FR |> subset(Day == "2022-08-02"), type = "FR", outdir = wd, outname = "Simulated_FR_Day1", base_raster = dist_col)
fr_d2_nfix <- nfixpercell_function(trajectories = traj_FR |> subset(Day == "2022-08-03"), type = "FR", outdir = wd, outname = "Simulated_FR_Day2", base_raster = dist_col)

extent_function(traj_FR |> subset(Day == "2022-08-02"), type = "FR", outdir = wd, outname = "Simulated_FR_Day1", base_raster = dist_col)
extent_function(traj_FR |> subset(Day == "2022-08-03"), type = "FR", outdir = wd, outname = "Simulated_FR_Day2", base_raster = dist_col)

save.image("Simulate_species.RData")









