#### From Lagrangian to Eulerian perspectives
#### Script: Functions to format and analyse tracking data 
#### Last modification: 04 04 2024
#### R version 4.2.3 (2023-03-15 ucrt) -- "Shortstop Beagle"



### 1 - Format the trajectories for RSF analysis ------------------------------
format_track_data <- function(traj, coord_col = NULL, col_sizes = c(8000, 12000), type = "CPF"){
  if(type == "CPF"){
    ### 1 - Define available locations  
    ## A - define buffer around colonies based on observed max distance
    smry_dist <- traj |> 
      group_by(BirdId) |> 
      dplyr::summarize(dist_col = max(dist_col, na.rm = T),
                       Colony = first(Colony2)) 
    ggplot(smry_dist) + geom_boxplot(aes(y = dist_col, fill = Colony))
    
    col_pts <- coord_col[c(1,2),] |>
      mutate(size = col_sizes, Colony = c("Colony 1", "Colony 2"),
             x = Lon, y = Lat) |>
      st_as_sf(coords = c("x","y"), crs = virtualecologist::custom_proj())
    
    col_buf_1 <- col_pts[1,] |> 
      st_buffer(dist = max(smry_dist[which(smry_dist$Colony == "Colony 1"), "dist_col"])+5) |> # create buffer from the max dist around colony + 5 additional km
      st_crop(., y = st_bbox(dist_col)) # crop to the study area
    col_buf_2 <- col_pts[2,] |>
      st_buffer(dist = max(smry_dist[which(smry_dist$Colony == "Colony 2"), "dist_col"])+5) |>
      st_crop(., y = st_bbox(dist_col))
    
    
    ## B - sample points within buffers 
    # sample 10 avail points for each positions
    # -> sample nbr of positions * 10, for each ind
    n_loc_ind <- traj |> 
      group_by(BirdId) |> 
      dplyr::summarize(n = n(), Colony = first(Colony2)) |> ungroup()
    bird.id <- unique(traj$BirdId)
    
    avail <- pbapply::pblapply(1:length(bird.id), function(i) {
      sub <- n_loc_ind |> subset(BirdId == bird.id[i])
      
      if(sub$Colony == "Colony 1") { buf <- col_buf_1 }
      if(sub$Colony == "Colony 2") { buf <- col_buf_2 }
      
      avail <- st_sample(x = buf, size = sub$n*10) |> 
        st_as_sf() %>%
        mutate(Lon = st_coordinates(.)[,1],
               Lat = st_coordinates(.)[,2],
               Colony = if(sub$Colony == "Colony 1") {1} else {2},
               Colony2 = sub$Colony,
               BirdId = bird.id[i])
      return(avail)
    }) |> 
      dplyr::bind_rows() |> 
      mutate(Used = FALSE) 
    
    
    ### 2 - Extract env values
    traj[, c("Resource1", "Resource2")] <- terra::extract(x = env$rasters, y = traj[, c("Lon", "Lat")])[, 2:3]
    avail[, c("Resource1", "Resource2")] <- terra::extract(x = env$rasters, y = avail[, c("Lon", "Lat")])[, 2:3]
    
    
    ### 3 - Extract dist col 
    coord_col_sf <- st_as_sf(x = coord_col, 
                             coords = c("Lon", "Lat"),
                             crs = virtualecologist::custom_proj())
    
    traj <- traj |> 
      mutate(x = Lon, y = Lat) |> 
      st_as_sf(coords = c("x", "y"), crs = virtualecologist::custom_proj()) %>%
      mutate(dist_col = case_when(
        Colony == 1 ~ st_distance(., coord_col_sf[1, 1:2])[,1],
        Colony == 2 ~ st_distance(., coord_col_sf[2, 1:2])[,1]
      ))
    avail <- avail %>%
      mutate(dist_col = case_when(
        Colony == 1 ~ st_distance(., coord_col_sf[1, 1:2])[,1],
        Colony == 2 ~ st_distance(., coord_col_sf[2, 1:2])[,1]
      ))
    
    
    ### 4 - Bind data
    traj$weight <- 1
    avail$weight <- 1e3
    final_df <- plyr::rbind.fill(st_drop_geometry(traj),
                                 st_drop_geometry(avail)) |> units::drop_units()
    
    return(list(buffers = rbind(col_buf_1, col_buf_2), 
                final_table = final_df))
  }
  if(type == "FR"){
    ### 1 - Define available locations  
    n_loc_ind <- traj |> 
      group_by(ID) |> 
      dplyr::summarize(n = n()) |> ungroup()
    bird.id <- unique(traj$ID)
    
    avail <- pbapply::pblapply(1:length(bird.id), function(i) {
      sub <- n_loc_ind |> subset(ID == bird.id[i])
      
      avail <- st_sample(x = env$rasters$sim1 |> ext() |> as.polygons() |> as_sf(), 
                         size = sub$n*10) |> 
        st_as_sf() %>%
        mutate(Lon = st_coordinates(.)[,1],
               Lat = st_coordinates(.)[,2],
               ID = bird.id[i])
      return(avail)
    }) |> 
      dplyr::bind_rows() |> 
      mutate(Used = FALSE) 
    
    
    ### 2 - Extract env values
    traj[, c("Resource1", "Resource2")] <- terra::extract(x = env$rasters, y = traj[, c("Lon", "Lat")])[, 2:3]
    avail[, c("Resource1", "Resource2")] <- terra::extract(x = env$rasters, y = avail[, c("Lon", "Lat")])[, 2:3]
    
    
    ### 3 - Bind data
    traj$weight <- 1
    avail$weight <- 1e3
    final_df <- plyr::rbind.fill(st_drop_geometry(traj),
                                 st_drop_geometry(avail)) |> units::drop_units()
    
    return(final_df)
  }
}




## 2 - Compute the RSF version ---------------------------------------------
pred.map.gam <- function(model, data, buffer, type = "CPF") {
  data <- data |> st_drop_geometry() %>% units::drop_units()
  pred.df <- predict(object = model, newdata = data, type = "response", se.fit = T)
  
  # rasterize the prediction (exponentiated to get the prob of use in natural scale)
  rast <- raster::rasterFromXYZ(cbind(data[, c("x", "y")], pred.df$fit ))
  rast.CI.low <- raster::rasterFromXYZ(cbind(data[, c("x", "y")], pred.df$fit - pred.df$se.fit) )
  rast.CI.high <- raster::rasterFromXYZ(cbind(data[, c("x", "y")], pred.df$fit + pred.df$se.fit) )
  
  # crop prediction raster to accessible envelop, and 
  if(type == "CPF"){
    rast <- raster::mask(rast, buffer, updatevalue = 0)
    rast.CI.low <- raster::mask(rast.CI.low, buffer, updatevalue = 0)
    rast.CI.high <- raster::mask(rast.CI.high, buffer, updatevalue = 0)
  }
  
  # normalise (so that the map sums to 1)
  rast <- rast / raster::cellStats(rast, "sum")
  rast.CI.low <- rast.CI.low / raster::cellStats(rast.CI.low, "sum")
  rast.CI.high <- rast.CI.high / raster::cellStats(rast.CI.high, "sum")
  
  # stack all 
  stack <- raster::stack(rast, rast.CI.low, rast.CI.high)
  names(stack) <- c("Mean", "CI low", "CI high")
  return(stack)
}

rsf_fitting_gam <- function(df, Species = "CPF", buffers = NULL, col_sizes = c(8000, 12000)){
  if(Species == "CPF"){
    ## 1 - Fit the GAM 
    rsf.mod <- gam(Used ~ s(Resource1 , k = 4) + s(Resource2, k = 4) + s(dist_col, k = 4), weights = weight, 
                   family = "binomial", data = df)
    
    gam_results <- contrib_func(model = rsf.mod, data = df) |> 
      as.data.frame() |> 
      mutate(Deviance = summary(rsf.mod)$dev.expl)
    
    gam_curves <- predict(rsf.mod, newdata = df, type = "terms")
    gam_curves <- cbind(df, gam_curves)
    
    ## 2 - Predict 
    coord_col_sf <- st_as_sf(x = coord_col, 
                             coords = c("Lon", "Lat"),
                             crs = virtualecologist::custom_proj())
    
    env_df_Col1 <- env_df %>%
      mutate(Lon = x, Lat = y) %>%
      st_as_sf(.,coords = c("Lon","Lat"), crs = virtualecologist::custom_proj())%>% 
      mutate(dist_col = st_distance(., coord_col_sf[1, 1:2])[,1],
             weight = 1)
    env_df_Col2 <- env_df %>%
      mutate(Lon = x, Lat = y) %>%
      st_as_sf(.,coords = c("Lon","Lat"), crs = virtualecologist::custom_proj())%>% 
      mutate(dist_col = st_distance(., coord_col_sf[2, 1:2])[,1],
             weight = 1)
    
    RSF.map1 <- pred.map.gam(model = rsf.mod, 
                         data = env_df_Col1, 
                         buffer = buffers |> subset(Colony == "Colony 1"),
                         type = "CPF")
    RSF.map2 <- pred.map.gam(model = rsf.mod, 
                         data = env_df_Col2, 
                         buffer = buffers |> subset(Colony == "Colony 2"),
                         type = "CPF")
    
    Ab.map <- (RSF.map1*col_sizes[1]) + (RSF.map2*col_sizes[2])
    
    writeRaster(Ab.map[[1]], paste0(workdir, "/", Species, "_RSF_abundance_map.grd"), overwrite = T)
    writeRaster(Ab.map[[2]], paste0(workdir, "/", Species, "_RSF_abundance_map_CIlow.grd"), overwrite = T)
    writeRaster(Ab.map[[3]], paste0(workdir, "/", Species, "_RSF_abundance_map_CIhigh.grd"), overwrite = T)
    
    ## 3 - Pearson's coef
    df_list <- list(as.data.frame(Ab.map, xy = T), 
                    as.data.frame(rsce_1$rasters, xy = T), 
                    as.data.frame(CPF_true, xy = T))
    binded_tables <- df_list %>% reduce(full_join, by =  c("x", "y"))
    
    corr_kde <- Hmisc::rcorr(as.matrix(binded_tables[, c("Mean", "mean")]), type = "pearson")
    corr_suit <- Hmisc::rcorr(as.matrix(binded_tables[, c("Mean", "suitability")]), type = "pearson")
    
    gam_results$Pearson_kde <- as.vector(corr_kde$r)[2]
    gam_results$Pearson_suitability <- as.vector(corr_suit$r)[2]
    
    
    ## 4 - Export
    return(list(RSF_results = gam_results,
                Model = rsf.mod,
                Curves = gam_curves,
                Abundance_map_df = as.data.frame(Ab.map, xy = T),
                Abundance_map = Ab.map,
                RSF_Col1_df = as.data.frame(RSF.map1, xy = T),
                RSF_Col2_df = as.data.frame(RSF.map2, xy = T),
                RSF_Col1 = RSF.map1,
                RSF_Col2 = RSF.map2))
  }
  if(Species == "FR"){
    ## 1 - Fit the GAM 
    rsf.mod <- gam(Used ~ s(Resource1 , k = 4) + s(Resource2, k = 4), weights = weight, 
                   family = "binomial", data = df)
    
    gam_results <- contrib_func(model = rsf.mod, data = df) |> 
      as.data.frame() |> 
      mutate(Deviance = summary(rsf.mod)$dev.expl)
    
    gam_curves <- predict(rsf.mod, newdata = df, type = "terms")
    gam_curves <- cbind(df, gam_curves)
    
    
    ## 2 - Predict 
    env_df <- env_df |> mutate(weight = 1)
    
    RSF.map <- pred.map.gam(model = rsf.mod, 
                        data = env_df, 
                        buffer = NULL,
                        type = "FR")
    
    Ab.map <- RSF.map*col_sizes
    
    writeRaster(Ab.map[[1]], paste0(workdir, "/", Species, "_RSF_abundance_map.grd"), overwrite = T)
    writeRaster(Ab.map[[2]], paste0(workdir, "/", Species, "_RSF_abundance_map_CIlow.grd"), overwrite = T)
    writeRaster(Ab.map[[3]], paste0(workdir, "/", Species, "_RSF_abundance_map_CIhigh.grd"), overwrite = T)
    
    ## 3 - Pearson's coef
    df_list <- list(as.data.frame(Ab.map, xy = T), 
                    as.data.frame(rsce_2$rasters, xy = T), 
                    as.data.frame(FR_true, xy = T))
    binded_tables <- df_list %>% reduce(full_join, by =  c("x", "y"))
    
    corr_kde <- Hmisc::rcorr(as.matrix(binded_tables[, c("Mean", "mean")]), type = "pearson")
    corr_suit <- Hmisc::rcorr(as.matrix(binded_tables[, c("Mean", "suitability")]), type = "pearson")
    
    gam_results$Pearson_kde <- as.vector(corr_kde$r)[2]
    gam_results$Pearson_suitability <- as.vector(corr_suit$r)[2]
    
    ## 4 - Export
    return(list(RSF_results = gam_results,
                Model = rsf.mod,
                Curves = gam_curves,
                Abundance_map_df = as.data.frame(Ab.map, xy = T),
                Abundance_map = Ab.map,
                RSF_df = as.data.frame(RSF.map, xy = T),
                RSF = RSF.map))
  }
}


