#### From Lagrangian to Eulerian perspectives
#### Script: Setting and launching surveys
#### Last modification: 15 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)

workdir <- "your working directory"
  
#### ---------------------------------- ####
#### --- 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 survey designs ---- ####
#### --------------------------------- ####

### 1 - Plane survey design ----------------------------------------------------
# generate the layout
plane_survey <- generate_survey_plan( bbx_xmin = 0, bbx_xmax = 90, bbx_ymin = 0, bbx_ymax = 90,
                                      start_x = 5, end_x = 85, start_y = 5, end_y = 85,
                                      space_out_factor = 4, segmentize = TRUE, seg_length = 1,
                                      buffer = TRUE, buffer_width = 1)

ggplot() +
  geom_spatraster(data = rsce_2$rasters) +
  geom_sf(data = plane_survey$buffered_segments)

# assign the flight plan
plane_survey_ready <- assign_flight_plan(sf_segments = plane_survey$buffered_segments,
                                         flight_id = c(1:80, 160:81, 161:240, 320:241, 321:400, 480:401, 481:560, 640:561, 641:720, 800:721, 801:880, 960:881, 961:1040,
                                                       1120:1041, 1121:1200, 1280:1201, 1281:1360, 1440:1361, 1441:1520, 1600:1521, 1601:1680),
                                         flight_day = "2022-08-02",
                                         survey_start_hour = "06:00:00",
                                         flight_speed = 160,
                                         intertransect_gap_duration = 22,
                                         col_trans_id = "transect")


ggplot(plane_survey_ready) +
  geom_sf(aes(fill = start_time),
          color = "black", size = 0.2) +
  scale_fill_datetime(low = "#8073ac", high = "#fdb863") +
  theme_bw() + theme(panel.grid = element_blank()) +
  labs(fill = "Hour", limits = c("06:00", "17:00"), x = "x", y = "y")+
  coord_sf(xlim = c(0,90), ylim = c(0,90))

ggplot(plane_survey_ready) +
  geom_sf(aes(fill = flight_id),
          color = "black", size = 0.2) +
  theme_bw() + theme(panel.grid = element_blank()) +
  coord_sf(xlim = c(0,90), ylim = c(0,90))






### 2 - Boat survey design ----------------------------------------------------
# generate the layout
boat_survey <- generate_survey_plan( bbx_xmin = 15, bbx_xmax = 75, bbx_ymin = 15, bbx_ymax = 75,
                                      start_x = 20, end_x = 70, start_y = 20, end_y = 70,
                                      space_out_factor = 7, segmentize = TRUE, seg_length = 1,
                                      buffer = TRUE, buffer_width = 1) 

ggplot() +
  geom_spatraster(data = rsce_2$rasters) +
  geom_sf(data = boat_survey$buffered_segments)

# assign the obs plan (two days)
boat_survey_ready_Day1 <- assign_flight_plan(sf_segments = boat_survey$buffered_segments |> subset(transect %in% c(1:4)),
                                         flight_id = c(1:50, 100:51, 101:150, 200:151),
                                         flight_day = "2022-08-02",
                                         survey_start_hour = "06:00:00",
                                         flight_speed = 18.5,
                                         intertransect_gap_duration = 68,
                                         col_trans_id = "transect")
boat_survey_ready_Day2 <- assign_flight_plan(sf_segments = boat_survey$buffered_segments |> subset(transect %in% c(5:8)),
                                             flight_id = c(1:50, 100:51, 101:150, 200:151),
                                             flight_day = "2022-08-03",
                                             survey_start_hour = "06:00:00",
                                             flight_speed = 18.5,
                                             intertransect_gap_duration = 68,
                                             col_trans_id = "transect")

boat_survey_ready <- rbind(boat_survey_ready_Day1, 
                           boat_survey_ready_Day2) |> 
  mutate(Day = as_date(start_time))


ggplot(boat_survey_ready) +
  geom_sf(aes(fill = start_time),
          color = "black", size = 0.2) +
  scale_fill_datetime(low = "#8073ac", high = "#fdb863") +
  theme_bw() + theme(panel.grid = element_blank()) +
  labs(fill = "Hour", limits = c("06:00", "17:00"), x = "x", y = "y")+
  coord_sf(xlim = c(0,90), ylim = c(0,90)) + 
  facet_wrap("Day")

ggplot(boat_survey_ready) +
  geom_sf(aes(fill = flight_id),
          color = "black", size = 0.2) +
  theme_bw() + theme(panel.grid = element_blank()) +
  coord_sf(xlim = c(0,90), ylim = c(0,90)) + 
  facet_wrap("Day")


save.image("Surveys.RData")



#### ------------------------- ####
#### ---- Run the surveys ---- ####
#### ------------------------- ####

### 1 - import trajectories ----------------------------------------------------
load(file = "Simulated_tracks_CPF.rda") 
load(file = "Simulated_tracks_FR.rda") 



### 2 - Plane survey with CPF --------------------------------------------------
plane_cpf <- launch_survey_on_movement(
  survey_data_buffered = plane_survey_ready,
  survey_data_linear = plane_survey$segments,
  traj_data = traj_CPF |> subset(Day == "2022-08-02"),
  track_id_col = "BirdId",
  seg_id_col = "seg_id",
  line_transect = TRUE, 
  detection_function = "hn",
  sigma = 0.2)




### 3 - Boat survey with CPF ---------------------------------------------------
boat_cpf <- launch_survey_on_movement(
  survey_data_buffered = boat_survey_ready,
  survey_data_linear = boat_survey$segments,
  traj_data = traj_CPF,
  track_id_col = "BirdId",
  seg_id_col = "seg_id",
  line_transect = TRUE, 
  detection_function = "hn",
  sigma = 0.5)

boat_cpf$obs_table$Day <- as_date(boat_cpf$obs_table$start_time_seg)


save.image("Surveys.RData")





### 4 - Plane survey with FR ---------------------------------------------------
plane_FR <- launch_survey_on_movement(
  survey_data_buffered = plane_survey_ready,
  survey_data_linear = plane_survey$segments,
  traj_data = traj_FR |> subset(Day == "2022-08-02"),
  track_id_col = "trip_id",
  seg_id_col = "seg_id",
  line_transect = TRUE, 
  detection_function = "hn",
  sigma = 0.2)


### 5 - Boat survey with FR ----------------------------------------------------
boat_FR <- launch_survey_on_movement(
  survey_data_buffered = boat_survey_ready,
  survey_data_linear = boat_survey$segments,
  traj_data = traj_FR,
  track_id_col = "trip_id",
  seg_id_col = "seg_id",
  line_transect = TRUE, 
  detection_function = "hn",
  sigma = 0.5)

boat_FR$obs_table$Day <- as_date(boat_FR$obs_table$start_time_seg)


save.image("Surveys.RData")




#### 6 - Plot ------------------------------------------------------------------
g1 <- ggplot() + geom_sf(data = plane_survey$transects) +
  geom_point(data = plane_cpf$obs_table |> 
               subset(detected==1) ,
             aes(x = Lon, y = Lat), color = "black", size = 0.7) +
  theme_bw() + theme(plot.title.position = "plot") +
  scale_color_colorblind(name = "Day") +
  geom_point(data = coord_col, aes(x = Lon, y = Lat), col = "red") +
  geom_text(data = coord_col, aes(x = Lon-2.6, y = Lat-2.6, label = Col, fontface = 2), col = "red") +
  labs(title = "Central-place forager - plane survey",
       subtitle = paste0("N = ", sum(plane_cpf$effort_table$N_ind_detected)),
       x = "x", y = "y") + 
  coord_sf(xlim = c(0,90), ylim = c(0,90))
g2 <- ggplot() + geom_sf(data = boat_survey$transects) +
  geom_point(data = boat_cpf$obs_table |> 
               subset(detected == 1),
             aes(x = Lon, y = Lat, color = as.factor(Day)), size = 0.7) +
  theme_bw() + theme(plot.title.position = "plot") +
  scale_color_colorblind(name = "Day") +
  geom_point(data = coord_col, aes(x = Lon, y = Lat), col = "red") +
  geom_text(data = coord_col, aes(x = Lon-2.6, y = Lat-2.6, label = Col, fontface = 2), col = "red") +
  labs(title = "Central-place forager - boat survey",
       subtitle = paste0("N = ", sum(boat_cpf$effort_table$N_ind_detected)),
       x = "x", y = "y") + 
  coord_sf(xlim = c(0,90), ylim = c(0,90))

g3 <- ggplot() + geom_sf(data = plane_survey$transects) +
  geom_point(data = plane_FR$obs_table |> 
               subset(detected==1) ,
             aes(x = Lon, y = Lat), color = "black", size = 0.7) +
  theme_bw() + theme(plot.title.position = "plot") +
  scale_color_colorblind(name = "Day") +
  labs(title = "Free ranger - plane survey",
       subtitle = paste0("N = ", sum(plane_FR$effort_table$N_ind_detected)),
       x = "x", y = "y") + 
  coord_sf(xlim = c(0,90), ylim = c(0,90))
g4 <- ggplot() + geom_sf(data = boat_survey$transects) +
  geom_point(data = boat_FR$obs_table |> 
               subset(detected == 1),
             aes(x = Lon, y = Lat, color = as.factor(Day)), size = 0.7) +
  theme_bw() + theme(plot.title.position = "plot") +
  scale_color_colorblind(name = "Day") +
  labs(title = "Free ranger - boat survey",
       subtitle = paste0("N = ", sum(boat_FR$effort_table$N_ind_detected)),
       x = "x", y = "y") + 
  coord_sf(xlim = c(0,90), ylim = c(0,90))

library(ggpubr)
ggarrange(g1 + ggtitle("Plane survey") + theme(plot.title.position = "panel"), 
          g2 + ggtitle("Boat survey") + theme(plot.title.position = "panel"), 
          ncol = 2, nrow = 1,
          common.legend = T, legend = "none") + bgcolor("white") + border("white")
ggsave(paste0(workdir, "/Virtual_survey_with_sightings_CPF.png"), height = 11.5, width = 20, units = "cm", dpi = 350)

ggarrange(g3 + ggtitle("Plane survey") + theme(plot.title.position = "panel"), 
          g4 + ggtitle("Boat survey") + theme(plot.title.position = "panel"), 
          ncol = 2, nrow = 1,
          common.legend = T, legend = "none") + bgcolor("white") + border("white")
ggsave(paste0(workdir, "/Virtual_survey_with_sightings_FR.png"), height = 11.5, width = 20, units = "cm", dpi = 350)




## plot detection functions as well
det_func_df <- rbind(
  boat_cpf$obs_table |> st_drop_geometry() |> 
    mutate(Species = "Central-place forager",
           Platform = "Boat"),
  plane_cpf$obs_table |> st_drop_geometry() |> 
    mutate(Day = as_date(start_time_seg),
           Species = "Central-place forager",
           Platform = "Plane"),
  
  boat_FR$obs_table |> st_drop_geometry() |> 
    mutate(Species = "Free ranger",
           Platform = "Boat"),
  plane_FR$obs_table |> st_drop_geometry() |> 
    mutate(Day = as_date(start_time_seg),
           Species = "Free ranger",
           Platform = "Plane")
  )

ggplot(det_func_df) +
  geom_line(aes(x = dist_seg, y = prob_dist)) +
  facet_grid(Platform~Species) +
  labs(x = "Distance to transect line", y = "Detection probability")+
  theme_bw() +
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold.italic", hjust = 0))
ggsave(paste0(workdir, "/Detection_functions.png"), 
       height = 10, width = 12, units = "cm", dpi = 350)



#### ------------------------------------- ####
#### ---- Export surveys for analysis ---- ####
#### ------------------------------------- ####

## export survey as a table
write.table(x = plane_cpf$effort_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_CPF_plane.csv"), dec = ".", sep = ";", row.names = F)
write.table(x = plane_cpf$obs_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_CPF_plane_obstable.csv"), dec = ".", sep = ";", row.names = F)

write.table(x = boat_cpf$effort_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_CPF_boat.csv"), dec = ".", sep = ";", row.names = F)
write.table(x = boat_cpf$obs_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_CPF_boat_obstable.csv"), dec = ".", sep = ";", row.names = F)

write.table(x = plane_FR$effort_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_FR_plane.csv"), dec = ".", sep = ";", row.names = F)
write.table(x = plane_FR$obs_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_FR_plane_obstable.csv"), dec = ".", sep = ";", row.names = F)

write.table(x = boat_FR$effort_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_FR_boat.csv"), dec = ".", sep = ";", row.names = F)
write.table(x = boat_FR$obs_table |> st_drop_geometry(), file = paste0(workdir, "/Virtual_survey_FR_boat_obstable.csv"), dec = ".", sep = ";", row.names = F)


save.image("Surveys.RData")



#### look at the number of time ind are sighted by a single survey
hist_sightings <- function(obs_table, Species, Survey){
  # ind available to detection
  t <- table(obs_table |> pull("track_id"))
  psce_freq <- as.data.frame(t)
  psce_freq$Survey <- Survey
  psce_freq$Species <- Species

  # ind detected
  t <- table(obs_table |> subset(detected == 1) |> pull("track_id"))
  sight_freq <- as.data.frame(t)
  sight_freq$Survey <- Survey
  sight_freq$Species <- Species
  
  # return
  return(list(Presence = psce_freq, Detected = sight_freq))
}

cpf_plane_hist <- hist_sightings(obs_table = plane_cpf$obs_table, Species = "Central-place forager", Survey = "Plane")
cpf_boat_hist <- hist_sightings(obs_table = boat_cpf$obs_table, Species = "Central-place forager", Survey = "Boat")
fr_plane_hist <- hist_sightings(obs_table = plane_FR$obs_table, Species = "Free ranger", Survey = "Plane")
fr_boat_hist <- hist_sightings(obs_table = boat_FR$obs_table, Species = "Free ranger", Survey = "Boat")

hist_freq <- rbind(cpf_plane_hist$Detected, cpf_boat_hist$Detected, fr_plane_hist$Detected, fr_boat_hist$Detected)
hist_freq2 <- rbind(cpf_plane_hist$Presence, cpf_boat_hist$Presence, fr_plane_hist$Presence, fr_boat_hist$Presence)


ggplot(as.data.frame(hist_freq)) + 
  geom_histogram(aes(x = Freq)) +
  facet_grid(Survey~Species)+
  labs(title = "Frequency of the number of time\na unique individual has been sighted",
       y = "Number of individuals", x = "Frequency") +
  theme_bw() + 
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold.italic", hjust = 0),
        plot.title.position = "plot",
        panel.grid = element_blank())
ggsave("Histogram_detected.png", height = 4.5, width = 4.5)


ggplot(as.data.frame(hist_freq2)) + 
  geom_histogram(aes(x = Freq)) +
  facet_grid(Survey~Species)+
  labs(title = "Frequency of the number of time\na unique individual was available to detection",
       y = "Number of individuals", x = "Frequency") +
  theme_bw() + 
  theme(strip.background = element_blank(),
        strip.text = element_text(face = "bold.italic", hjust = 0),
        plot.title.position = "plot",
        panel.grid = element_blank())
ggsave("Histogram_available_to_detection.png", height = 4.5, width = 4.5)
