aix_2024_autumn

——————————- Loading libraries and files ——————————-

library(sp)
library(ctmm)

## 
## Attachement du package : 'ctmm'

## L'objet suivant est masqué depuis 'package:sp':
## 
##     plot

library(move)

## Le chargement a nécessité le package : geosphere

## Le chargement a nécessité le package : raster

## 
## Attachement du package : 'raster'

## Les objets suivants sont masqués depuis 'package:ctmm':
## 
##     projection, projection<-

## 
## Attachement du package : 'move'

## L'objet suivant est masqué depuis 'package:ctmm':
## 
##     speed

library(lattice)
library(solartime) # calcul de la position du soleil
library(readxl)
library(adehabitatLT) # pour plot les points et trajets

## Le chargement a nécessité le package : ade4

## Le chargement a nécessité le package : adehabitatMA

## Registered S3 methods overwritten by 'adehabitatMA':
##   method                       from
##   print.SpatialPixelsDataFrame sp  
##   print.SpatialPixels          sp

## 
## Attachement du package : 'adehabitatMA'

## L'objet suivant est masqué depuis 'package:raster':
## 
##     buffer

## 
## Attachement du package : 'adehabitatLT'

## L'objet suivant est masqué depuis 'package:move':
## 
##     burst

library(geosphere)
library(terra)     # pour l'image satellite

## terra 1.9.11

## 
## Attachement du package : 'terra'

## L'objet suivant est masqué depuis 'package:adehabitatMA':
## 
##     buffer

## Les objets suivants sont masqués depuis 'package:ctmm':
## 
##     distance, meta

library(basemaps)  # pour l'image satellite
library(sf)        # pour st_bbox / st_crs

## Linking to GEOS 3.14.1, GDAL 3.12.1, PROJ 9.7.1; sf_use_s2() is TRUE

PATH <- "C:/Users/ppm/Documents/Victor_2026/Data/"
general <- read_excel(paste(PATH,"Summary GPS data.xlsx",sep=""))

## New names:
## • `TH_s` -> `TH_s...23`
## • `TH_s` -> `TH_s...28`

# Tous les fichiers
NAMES0 <- paste("cleaned data/",na.omit(general$File_clean),sep="")
NAMES1 <- paste("filtered_data1/",na.omit(general$File_filter1),sep="")
NAMES3 <- paste("filtered_data3/",na.omit(general$File_filter3),sep="")
NAMES <- c(NAMES0,NAMES1,NAMES3)
NAMES

##  [1] "cleaned data/34620_Hope_Lotekv3"                
##  [2] "cleaned data/34621_Alvin_Lotekv3"               
##  [3] "cleaned data/ID002_Victory_Xeriusv3"            
##  [4] "cleaned data/ID003_Valentine_Xeriusv3"          
##  [5] "cleaned data/34619_Feeby_Lotekv3"               
##  [6] "cleaned data/ID001_Feeby_Xeriusv3"              
##  [7] "cleaned data/ID045_Rocky_Xeriusv3"              
##  [8] "cleaned data/ID042_Rainy_Xeriusv3"              
##  [9] "cleaned data/ID044_Crunchy_Xeriusv3"            
## [10] "cleaned data/ID043_Cloudy_Xeriusv3"             
## [11] "cleaned data/ID044_CrunchyAurore_Xeriusv3"      
## [12] "cleaned data/ID004_Wisk_Xeriusv3"               
## [13] "cleaned data/ID009_Wisk_Xeriusv3"               
## [14] "cleaned data/ID012_Wisk_Xeriusv3"               
## [15] "cleaned data/ID005_Fingers_Xeriusv3"            
## [16] "cleaned data/ID007_Sleepy_Xeriusv3"             
## [17] "cleaned data/ID021_Sleepy_Xeriusv3"             
## [18] "cleaned data/ID006_Corty_Xeriusv3"              
## [19] "cleaned data/ID008_Tosca_Xeriusv3"              
## [20] "cleaned data/ID010_Rox_Xeriusv3"                
## [21] "cleaned data/ID016_Crispy_Xeriusv3"             
## [22] "cleaned data/ID011_Rouky_Xeriusv3"              
## [23] "cleaned data/ID017_Roger_Xeriusv3"              
## [24] "cleaned data/ID014_Roger_Xeriusv3"              
## [25] "cleaned data/ID013_Scar_Xeriusv3"               
## [26] "cleaned data/ID022_Igor_Xeriusv3"               
## [27] "cleaned data/57406_Igor_Lotekv3"                
## [28] "cleaned data/ID015_Ciri_Xeriusv3"               
## [29] "cleaned data/57405_Ciri_Lotekv3"                
## [30] "cleaned data/ID025_Ciri_Xeriusv3"               
## [31] "cleaned data/ID018_Sandra_Xeriusv3"             
## [32] "cleaned data/ID015_Papillon_Xeriusv3"           
## [33] "cleaned data/ID028_Papillon_Xeriusv3"           
## [34] "cleaned data/ID022_Squeeky_Xeriusv3"            
## [35] "cleaned data/ID023_Squeeky_Xeriusv3"            
## [36] "cleaned data/ID029_Fly_Xeriusv3"                
## [37] "cleaned data/ID024_Crepe_Xeriusv3"              
## [38] "cleaned data/ID030_Maria_Xeriusv3"              
## [39] "cleaned data/ID031_Cally_Xeriusv3"              
## [40] "cleaned data/ID027_Peanut_Xeriusv3"             
## [41] "cleaned data/ID026_Adele_Xeriusv3"              
## [42] "cleaned data/ID037_Calico_Xeriusv3"             
## [43] "cleaned data/ID032_Noisette_Xeriusv3"           
## [44] "cleaned data/ID019_Pomme_Xeriusv3"              
## [45] "cleaned data/67377_Speedy_Lotekv3"              
## [46] "cleaned data/67378_Lucky_Lotekv3"               
## [47] "cleaned data/ID019_PommeLucky_Xeriusv3"         
## [48] "filtered_data1/-"                               
## [49] "filtered_data1/-"                               
## [50] "filtered_data1/ID002_Victory_Xeriusv3_Alt"      
## [51] "filtered_data1/ID003_Valentine_Xeriusv3_Alt"    
## [52] "filtered_data1/-"                               
## [53] "filtered_data1/-"                               
## [54] "filtered_data1/ID045_Rocky_Xeriusv3_Alt"        
## [55] "filtered_data1/ID042_Rainy_Xeriusv3_Alt"        
## [56] "filtered_data1/ID044_Crunchy_Xeriusv3_Alt"      
## [57] "filtered_data1/ID043_Cloudy_Xeriusv3_Alt"       
## [58] "filtered_data1/ID044_CrunchyAurore_Xeriusv3_Alt"
## [59] "filtered_data1/ID004_Wisk_Xeriusv3_Alt"         
## [60] "filtered_data1/ID009_Wisk_Xeriusv3_Alt"         
## [61] "filtered_data1/ID012_Wisk_Xeriusv3_Alt"         
## [62] "filtered_data1/ID005_Fingers_Xeriusv3_Alt"      
## [63] "filtered_data1/ID007_Sleepy_Xeriusv3_Alt"       
## [64] "filtered_data1/-"                               
## [65] "filtered_data1/ID006_Corty_Xeriusv3_Alt"        
## [66] "filtered_data1/ID008_Tosca_Xeriusv3_Alt"        
## [67] "filtered_data1/ID010_Rox_Xeriusv3_Alt"          
## [68] "filtered_data1/ID016_Crispy_Xeriusv3_Alt"       
## [69] "filtered_data1/ID011_Rouky_Xeriusv3_Alt"        
## [70] "filtered_data1/ID017_Roger_Xeriusv3_Alt"        
## [71] "filtered_data1/ID014_Roger_Xeriusv3_Alt"        
## [72] "filtered_data1/ID013_Scar_Xeriusv3_Alt"         
## [73] "filtered_data1/-"                               
## [74] "filtered_data1/57406_Igor_Lotekv3_Alt"          
## [75] "filtered_data1/ID015_Ciri_Xeriusv3_Alt"         
## [76] "filtered_data1/57405_Ciri_Lotekv3_Alt"          
## [77] "filtered_data1/-"                               
## [78] "filtered_data1/ID018_Sandra_Xeriusv3_Alt"       
## [79] "filtered_data1/ID015_Papillon_Xeriusv3_Alt"     
## [80] "filtered_data1/-"                               
## [81] "filtered_data1/-"                               
## [82] "filtered_data1/-"                               
## [83] "filtered_data1/-"                               
## [84] "filtered_data1/-"                               
## [85] "filtered_data3/57406_Igor_Lotekv3_Alt_SpAn"     
## [86] "filtered_data3/57405_Ciri_Lotekv3_Alt_SpAn"

# Uniquement les fichiers de Aix - 2024 - Autumn
NAMES_sub <- NAMES[c(59,60,62,63,17,65,66,67,68,69,70,72)]
NAMES_sub

##  [1] "filtered_data1/ID004_Wisk_Xeriusv3_Alt"   
##  [2] "filtered_data1/ID009_Wisk_Xeriusv3_Alt"   
##  [3] "filtered_data1/ID005_Fingers_Xeriusv3_Alt"
##  [4] "filtered_data1/ID007_Sleepy_Xeriusv3_Alt" 
##  [5] "cleaned data/ID021_Sleepy_Xeriusv3"       
##  [6] "filtered_data1/ID006_Corty_Xeriusv3_Alt"  
##  [7] "filtered_data1/ID008_Tosca_Xeriusv3_Alt"  
##  [8] "filtered_data1/ID010_Rox_Xeriusv3_Alt"    
##  [9] "filtered_data1/ID016_Crispy_Xeriusv3_Alt" 
## [10] "filtered_data1/ID011_Rouky_Xeriusv3_Alt"  
## [11] "filtered_data1/ID017_Roger_Xeriusv3_Alt"  
## [12] "filtered_data1/ID013_Scar_Xeriusv3_Alt"

——————————- HELPER FUNCTION ——————————-

Find the actual file path from the filename

find_file_path <- function(short_name, name_list, PATH) {
  # Searches which element in name_list contains the string "short_name"
  # Example: short_name="Fingers", name_list="filtered_data1/ID005_Fingers_Xeriusv3_Alt" -> Found
  
  if (!grepl("/$", PATH)) {
    PATH <- paste0(PATH, "/")
  }
  
  matches <- grep(short_name, name_list, value = TRUE)
  
  if (length(matches) == 0) {
    stop(paste("Error: No file found containing the name '", short_name, "' in the provided list.", sep=""))
  }
  
  if (length(matches) > 1) {
    warning(paste("Warning: Multiple files found for '", short_name, "'. The first one will be used: ", matches[1], sep=""))
  }
  
  # Add the .R extension if not already present in the short name provided
  if (!grepl("\\.R$", matches[1])) {
    path <- paste0(matches[1], ".R")
  } else {
    path <- matches[1]
  }
  
  # 4. Concaténation finale : PATH + Sous-chemin/Fichier
  full_path <- paste0(PATH, path)
  
  return(full_path)
}

# Exemple d'utilisation
find_file_path("Crispy", NAMES_sub, PATH)

## [1] "C:/Users/ppm/Documents/Victor_2026/Data/filtered_data1/ID016_Crispy_Xeriusv3_Alt.R"

find_file_path("Rainy", NAMES, PATH)

## Warning in find_file_path("Rainy", NAMES, PATH): Warning: Multiple files found
## for 'Rainy'. The first one will be used: cleaned data/ID042_Rainy_Xeriusv3

## [1] "C:/Users/ppm/Documents/Victor_2026/Data/cleaned data/ID042_Rainy_Xeriusv3.R"

——————————- SIMPLE FUNCTIONS ——————————-

# Plot all GPS points and trajectories
plot_all_points <- function(animal_name, NAMES_sub, PATH) {
  # animal_name: The short name (e.g., "ID005_Fingers" or just "Fingers")
  
  # 1. Find the full file path
  full_path <- find_file_path(animal_name, NAMES_sub, PATH)

  # 2. Load the data and extract sub-objects
  # Verify the file exists physically on the disk before loading
  if (!file.exists(full_path)) {
    stop(paste("Error: The file '", full_path, "' does not exist on the disk.", sep=""))
  }
  
  loaded_objects <- load(full_path)
  data <- obj[[1]]       # Raw dataframe
  data_sq <- obj[[2]]    # ltraj object (trajectory)
  
  # 3. Plotting
  # x11()
  
  # Dynamic title
  main_title <- paste("Trajectory:", animal_name, "\n(N =", nrow(data_sq[[1]]), "points)")

  plot(data_sq, main = main_title, xlab = "Longitude",  ylab = "Latitude") # Début = trianble bleu, Fin = carré rouge
}

# Plot all GPS points colored by the time of the day (day/night/crep)
plot_day_night_points <- function(animal_name, NAMES_sub, PATH){
  # animal_name: The short name (e.g., "ID005_Fingers" or just "Fingers")
  
  # 1. Find the full file path
  full_path <- find_file_path(animal_name, NAMES_sub, PATH)
  
  # 2. Load the data and extract sub-objects
  # Verify the file exists physically on the disk before loading
  if (!file.exists(full_path)) {
    stop(paste("Error: The file '", full_path, "' does not exist on the disk.", sep=""))
  }
  
  loaded_objects <- load(full_path)
  data <- obj[[1]]       # Raw dataframe
  data_sq <- obj[[2]]    # ltraj object (trajectory)

  # 3. Calculate Day/Night Periods
  
  # Calculate the sun's altitude in the sky
  # Here, the elevation is calculated only from the first position only. It doesn’t matter, as the critical angle is small (a few kilometres makes no difference in relation to the sun)
  elev <- computeSunPosition(data_sq[[1]]$date,data_sq[[1]]$y,data_sq[[1]]$x)[,"elevation"]
  
  # Retrieve the animal's raw file name from the 'general' summary dataframe
  filename   <- basename(full_path)
  parts      <- strsplit(filename, "_")[[1]]
  name_short <- paste(parts[1], parts[2], sep = "_")
  
  # Retrieve the sun elevation thresholds (day_start, night_start) from 'general' 
  day_night_param <- as.numeric(general[which(general$File_raw==name_short),c("elev_day", "elev_night")])
  
  # Assign "day", "night", or "crep" (crepuscule) based on sun elevation
  # Logic: 
  # - If elevation > day_threshold -> "day"
  # - If elevation < night_threshold -> "night"
  # - Otherwise -> "crep"
  data_sq[[1]]$day_night <- c("crep","day")[1+(elev>day_night_param[1])]
  data_sq[[1]]$day_night[elev<( day_night_param[2])] <- "night"  
  
  # 4. Plotting
  # Plot the animal's positions colored by day/night status  
  # plot(data_sq[[1]]$x,data_sq[[1]]$y,col=as.numeric(as.factor(data_sq[[1]]$day_night)),pch=20,cex=2,main=name_short)
  
  plot(data_sq[[1]]$x, data_sq[[1]]$y, 
       col = as.numeric(as.factor(data_sq[[1]]$day_night)), 
       pch = 20,      # Solid circles
       cex = 2,       # Point size
       main = name_short, xlab = "Longitude", ylab = "Latitude")
  legend("bottomright", legend = levels(as.factor(data_sq[[1]]$day_night)),
         col = 1:3, 
         pch = 20, 
         title = "Period")
}

——————————- MAIN FUNCTION ——————————-

Identifies potential nesting sites by clustering GPS points recorded during “night” periods

plot_night_clusters <- function(animal_name, NAMES_sub, PATH, distance_threshold, nb_points_to_cluster) {
  
  
  # ---------------------------------------------------------------------------
  # 1. DATA LOADING
  # ---------------------------------------------------------------------------
  
  # 1.1 Find the full file path based on the animal's name
  full_path <- find_file_path(animal_name, NAMES_sub, PATH)
  
  # 1.2 Verify the file exists physically on the disk before attempting to load
  if (!file.exists(full_path)) {
    stop(paste("Error: The file '", full_path, "' does not exist on the disk."))
  }
  
  # 1.3 Load the data and extract sub-objects from the loaded list
  loaded_objects <- load(full_path)
  data <- obj[[1]]       # Raw data
  data_sq <- obj[[2]]    # 'ltraj' object
  
  # ---------------------------------------------------------------------------
  # 2. DAY/NIGHT CLASSIFICATION (Astronomical)
  # ---------------------------------------------------------------------------
  
  # Calculate the sun's elevation (angle above the horizon) for every GPS point.
  # Note: The function 'computeSunPosition' expects a single latitude/longitude (scalar) 
  # and a vector of dates. Although we pass vectors, it only uses the first coordinate.
  # This is still valid because the home range is small (few km), so the 
  # sun's position relative to the animal is effectively constant across the path.
  elev <- computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y, data_sq[[1]]$x)[,"elevation"]
  
  # Retrieve the animal's raw file name from the 'general' summary datafram
  filename   <- basename(full_path)
  parts      <- strsplit(filename, "_")[[1]]
  name_short <- paste(parts[1], parts[2], sep = "_")
  
  # Retrieve the sun elevation thresholds (day_start, night_start) specific to this animal.
  # These values define the angle at which day transitions to crepuscule or night.
  day_night_param <- as.numeric(general[which(general$File_raw == name_short), c("elev_day", "elev_night")])
  
  # Assign "day", "night", or "crep" (crepuscule) based on sun elevation.
  # Logic:
  # - If elevation > day_threshold -> "day"
  # - If elevation < night_threshold -> "night"
  # - Otherwise -> "crep"
  data_sq[[1]]$day_night <- c("crep", "day")[1 + (elev > day_night_param[1])]
  data_sq[[1]]$day_night[elev < day_night_param[2]] <- "night"
  
  # ---------------------------------------------------------------------------
  # 3. DATA CONVERSION (adehabitatLT -> ctmm)
  # ---------------------------------------------------------------------------
  # Convert the trajectory from the classic 'adehabitatLT' format to the 'ctmm' format.
  # This requires an intermediate step using the 'move' package, which is the 
  # modern standard for handling animal movement data.
  
  # Convert to 'move' object (standardized format)
  squirrel_move <- move::move(x = data_sq, sensor = "GPS")
  
  # Explicitly set the Coordinate Reference System (CRS) to WGS84 (Global GPS standard).
  # This is crucial because subsequent distance calculations must be in meters, not degrees.
  crs(squirrel_move) <- CRS("+proj=longlat +ellps=WGS84")
  
  # Convert to 'ctmm' telemetry object (Continuous-Time Movement Model).
  # This format prepares the data for advanced statistical modeling.
  squirrel_move_ctmm <- as.telemetry(squirrel_move, timeformat = "auto", timezone = "UTC")
  
  # ---------------------------------------------------------------------------
  # 4. CLUSTERING NIGHT POINTS
  # ---------------------------------------------------------------------------
  
  # 4.1 Filter only "night" and "crep" points
  # We focus exclusively on periods when the animal is likely resting/nesting.
  night_sub <- data_sq[[1]][data_sq[[1]]$day_night %in% c("night", "crep"), ]
  
  # 4.2 Prepare the data frame for spatial analysis
  names(night_sub)[3] <- "time"             # Rename 'date' to 'time' (some spatial functions expect this name)
  night_sub$id <- row.names(night_sub)      # Add a unique row ID for tracking
  
  # 4.2 (suite) Assign a unique night ID to each group of consecutive points.
  # Logic: compute the time gap (in seconds) between each consecutive point.
  # If the gap exceeds 8 hours (28 800 s), it means a new night has started.
  # 'cumsum' converts the TRUE/FALSE gap-break vector into an incrementing integer ID.
  gap_threshold_s <- 8 * 3600  # 8 hours in seconds
  
  time_gaps <- c(0, diff(as.numeric(night_sub$time)))            # gap before each point (0 for the first)
  night_sub$night_id <- cumsum(time_gaps > gap_threshold_s) + 1  # starts at 1
  
  # Propagate night_id into the SpatialPointsDataFrame (built right after)
  # We store it temporarily; it will be added to xy after step 4.3.
  
  # 4.3 Convert to 'SpatialPointsDataFrame'
  # The 'geosphere' package requires spatial objects to calculate accurate geodesic distances.
  # We convert the simple x/y coordinates into a spatial object with WGS84 projection.
  x <- night_sub$x
  y <- night_sub$y
  
  xy <- SpatialPointsDataFrame(matrix(c(x, y), ncol = 2), data.frame(ID = seq(1:length(x))), proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84"))
  
  # Attach the night_id vector to the spatial object so it is available in step 6
  xy$night_id <- night_sub$night_id
  
  # 4.4 Calculate the geodesic distance matrix
  # 'distm' calculates the distance (in meters) between EVERY pair of points on the sphere.
  # Result: A square matrix where cell [i, j] = distance between point i and point j.
  mdist <- distm(xy)
  
  # 4.5 Perform Hierarchical Clustering
  # 'hclust' groups points based on similarity (distance).
  # Method "complete" (complete linkage) defines the distance between two clusters as the distance between their most distant points.
  # Advantage: This creates compact, spherical clusters. A point is only included if it is within the threshold distance 
  # of *all* points in the cluster, preventing elongated or noisy groups.
  hc <- hclust(as.dist(mdist), method = "complete")
  
  # 4.6 Define clusters based on the distance threshold 'cutree' cuts the hierarchical tree at a specific height (distance).
  # Points connected by a path shorter than 'distance_threshold' are assigned the same ID.
  xy$clust <- cutree(hc, h = distance_threshold)
  
  # ---------------------------------------------------------------------------
  # 5. FILTERING VALID NESTS
  # ---------------------------------------------------------------------------
  
  # 5.1 Count the number of points in each cluster
  # 'summary' counts occurrences of each cluster ID.
  HCS <- summary(as.factor(xy$clust))
  HCS <- data.frame(clust_id = as.numeric(names(HCS)), clust_nb = as.vector(HCS))
  
  # 5.2 Apply the minimum points filter
  # Logic: A valid nest must have enough GPS fixes. A single point might be a GPS error; multiple points suggest the animal stayed there.
  # A cluster (= nest) is considered to exist when there are at least “nb_points_to_cluster” points.
  xy$nest <- xy$clust
  xy$nest[xy$nest %in% HCS$clust_id[HCS$clust_nb < nb_points_to_cluster]] <- NA
  
  # ---------------------------------------------------------------------------
  # 6. VISUALIZATION
  # ---------------------------------------------------------------------------
  
  # --- 6.0 PREPARATION: Basemap & Projection ---
  
  # Extract coordinates from the spatial object
  coords <- coordinates(xy)
  
  # Define bounding box with a small margin (0.0005 degrees ~ 55 meters)
  margin <- 0.0005
  ext_bbox <- st_bbox(
    c(xmin = min(coords[,1]) - margin, xmax = max(coords[,1]) + margin,
      ymin = min(coords[,2]) - margin, ymax = max(coords[,2]) + margin),
    crs = st_crs(4326) # WGS84
  )
  
  # Download satellite basemap from Mapbox
  # Token provided in your script
  MAPBOX_TOKEN <- "pk.eyJ1IjoiY2VjYWxiIiwiYSI6ImNscnVvaG9rMjBqcXcyaXRlOTJheWpzc2cifQ.VplfhX9CpEQC-BdD2u815w"
  
  # 'basemap_terra' fetches the map tiles and creates a SpatRaster object
  base <- basemap_terra(ext_bbox, map_service = "mapbox", map_type = "satellite", map_token = MAPBOX_TOKEN)
  
  # Convert points to 'SpatVector' (terra format) for projection
  xy_vect <- vect(coords, crs = "+proj=longlat +datum=WGS84")
  
  # Project points to match the basemap CRS (usually Web Mercator, units in meters)
  # This ensures distances and shapes are accurate on the flat map display
  xy_vect_proj <- project(xy_vect, base)
  
  # Extract projected coordinates matrix for base R plotting functions (lines, symbols)
  coords_proj <- crds(xy_vect_proj)
  
  # Retrieve the day/night classification vector in the same order as coordinates
  # Needed to assign different shapes (pch) to night vs crepuscule
  day_night_sub <- data_sq[[1]][data_sq[[1]]$day_night %in% c("night","crep"), "day_night"]
  
  # --- DYNAMIC WIDTH CALCULATIONS For the legend position ---
  # Retrieve the current graph boundaries (after projection)
  base_ext <- ext(base)
  xmin_map <- base_ext[1]
  xmax_map <- base_ext[2]
  ymin_map <- base_ext[3]
  ymax_map <- base_ext[4]
  
  # Total width of the map
  map_width_m <- xmax_map - xmin_map
  map_height_m <- ymax_map - ymin_map
  
  # Setting positions as a percentage of the width/height
  # X-axisegend: Starts at 80% of the width (i.e. within the right-hand 20%) 
  legend_start_x <- xmin_map + (map_width_m * 0.80)
  # Y-axis legend: At the top, with a 10% margin from the top
  legend_start_y <- ymax_map - (map_height_m * 0.10)
  
  # Scale: Starts at 10% of the width (i.e. within the left-hand 10%)
  scale_start_x <- xmin_map + (map_width_m * 0.10)
  # At the bottom, with a 5% margin from the bottom
  scale_start_y <- ymin_map + (map_height_m * 0.05)
  
  # --- 6.1 PLOT BASEMAP, TRAJECTORIES & POINTS ---
  
  # Initialize the plot with the satellite image
  terra::plot(base)
  
  # Add the title manually AFTER the raster is drawn, so it is never hidden
  # 'line = 1' places it just above the plot area; 'font = 2' = bold
  mtext(paste("Night clusters –", animal_name),
        side = 3, line = 1, cex = 1.2, font = 2, col = "white")
  
  # Draw one trajectory per night rather than a single polyline over all points.
  # Each night gets a distinct color so consecutive nights are easy to tell apart.
  night_ids    <- sort(unique(xy$night_id))
  night_pal    <- rainbow(length(night_ids), alpha = 0.85)  # one color per night
  
  for (n in seq_along(night_ids)) {
    idx_night <- which(xy$night_id == night_ids[n])
    
    # Need at least 2 points to draw a line segment
    if (length(idx_night) >= 2) {
      lines(coords_proj[idx_night, ], col = night_pal[n], lwd = 1.2)
    }
  }
  
  # Define plotting symbols (pch): 16 = solid circle, 17 = solid triangle
  # Logic: If "night" -> circle (16), else (crep) -> triangle (17)
  pch_vec <- ifelse(day_night_sub == "night", 16, 17)
  
  # Plot the points
  # 'col = "white"' ensures visibility on dark backgrounds
  # 'pch = pch_vec' applies the shape logic defined above
  # 'cex = 0.7' controls point size
  terra::points(xy_vect_proj, col = "white", pch = pch_vec, cex = 0.8)
  
  # --- 6.2 PLOT VALID CLUSTERS ---
  
  # Identify unique IDs of valid nests (excluding NA)
  valid_nests <- unique(xy$nest[!is.na(xy$nest)])
  
  # Generate a color palette for the clusters (rainbow)
  # We will apply transparency later using 'adjustcolor'
  nest_colors <- rainbow(length(valid_nests))
  
  # Check the unit of the projected CRS to handle radius conversion correctly
  # Web Mercator (EPSG:3857) uses meters. If degrees, conversion is needed.
  crs_unit <- crs(base, proj = TRUE)
  
  if (length(valid_nests) > 0) {
    for (i in seq_along(valid_nests)) {
      nest_id <- valid_nests[i]
      
      # Find indices of points belonging to this cluster
      idx <- which(xy$nest == nest_id)
      
      # Extract coordinates (both projected for drawing, and lat/lon for distance calc)
      nest_coords_ll <- coords[idx, , drop = FALSE]      # Lat/Lon
      nest_coords_pr <- coords_proj[idx, , drop = FALSE] # Projected (meters)
      
      # Calculate cluster center (mean of projected coordinates)
      cx <- mean(nest_coords_pr[, 1])
      cy <- mean(nest_coords_pr[, 2])
      
      # Calculate radius: Max geodesic distance from center to any point in the cluster
      # 'distm' returns meters. We take the max to ensure all points are inside.
      center_ll <- cbind(mean(nest_coords_ll[,1]), mean(nest_coords_ll[,2]))
      dists <- distm(center_ll, nest_coords_ll)
      radius_m <- max(dists)
      
      # Convert radius to plot units
      # If CRS is in meters (like Web Mercator), use radius_m directly.
      # If CRS is in degrees, approximate conversion (1 deg ~ 111320 m).
      radius_plot <- if (grepl("units=m", crs_unit)) radius_m else radius_m / 111320
      
      # Draw the circle using base R 'symbols' function
      # 'circles': radius value
      # 'inches = FALSE': interpret radius in user coordinates (map units), not inches
      # 'add = TRUE': draw on top of existing plot
      # 'bg': background fill color (transparent)
      # 'fg': border color (more opaque)
      symbols(cx, cy, 
              circles = radius_plot, 
              inches = FALSE, 
              add = TRUE,
              bg = adjustcolor(nest_colors[i], alpha.f = 0.25), # 25% opaque fill
              fg = adjustcolor(nest_colors[i], alpha.f = 0.90), # 90% opaque border
              lwd = 2)
      
      # Add cluster ID label in the center
      text(cx, cy, labels = paste0("N", nest_id), cex = 0.9, font = 2, col = "white")
    }
  }
  
  # --- 6.3 ADD SCALE BAR (Dynamic Position) ---
  
  # Define a target scale length (~20% of plot width)
  scale_m_raw <- map_width_m * 0.20
  
  # Round to a "nice" number (1, 2, 5, 10, 20, 50, 100, etc.)
  # This makes the scale bar easy to read mentally (e.g., "100 m" instead of "87.3 m")
  magnitude <- 10^floor(log10(scale_m_raw))
  scale_m <- magnitude * round(scale_m_raw / magnitude)
  
  # Dynamic X position (defined above: 10% from the left)
  x0_scale <- scale_start_x
  x1_scale <- x0_scale + scale_m
  
  # Draw the scale bar rectangle
  # Height is dynamically calculated as 0.5% of plot height to remain proportional
  bar_height <- map_height_m * 0.005
  rect(x0_scale, scale_start_y - bar_height/2, 
       x1_scale, scale_start_y + bar_height/2, 
       col = "white", border = "black", lwd = 1.5)
  
  # Add text label above the bar
  # Format: "X m" or "X km" if >= 1000
  scale_label <- if (scale_m >= 1000) paste0(scale_m/1000, " km") else paste0(scale_m, " m")
  text((x0_scale + x1_scale) / 2, scale_start_y + bar_height * 1.5, 
       labels = scale_label, col = "white", cex = 0.8, font = 2)
  
  # --- 6.4 LEGEND (Dynamic Position) ---
  
  # Build legend entries in 4 groups:
  # (1) Point types: night circle + crepuscule triangle
  # (2) One colored line per night
  # (3) One colored circle per valid nest
  
  legend_labels <- c(
    "Night", "Crepuscule",            # group 1 : point types
    paste0("Night ", night_ids),              # group 2 : one entry per night trajectory
    paste0("Nest ",  valid_nests)             # group 3 : one entry per nest (empty if none)
  )
  
  legend_pch <- c(
    16, 17,                                   # group 1 : solid circle, solid triangle
    rep(NA, length(night_ids)),               # group 2 : lines, no point symbol
    rep(1,  length(valid_nests))              # group 3 : open circle for nest border
  )
  
  legend_lty <- c(
    NA, NA,                                   # group 1 : no line
    rep(1, length(night_ids)),                # group 2 : solid line
    rep(NA, length(valid_nests))              # group 3 : no line
  )
  
  legend_lwd <- c(
    NA, NA,                                   # group 1
    rep(1.5, length(night_ids)),              # group 2 : same lwd as trajectories
    rep(2,   length(valid_nests))             # group 3 : same lwd as circle borders
  )
  
  legend_col <- c(
    "black", "black",                         # group 1 : black on white legend bg
    night_pal,                                # group 2 : one color per night (same palette as lines)
    adjustcolor(nest_colors, alpha.f = 0.90)  # group 3 : same colors as nest circles
  )
  
  # Draw the legend
  # 'bg': White background with 85% opacity (makes text readable regardless of map darkness)
  # 'box.col': Dark gray border for the legend box
  # 'text.col': Black text for maximum contrast
  # x = 80% of the width from the left
  # y = 85% of the height from the bottom (i.e. 15% from the top)
  legend(x = legend_start_x, y = legend_start_y,
         legend   = legend_labels,
         pch      = legend_pch, 
         lty      = legend_lty, 
         lwd      = legend_lwd, 
         col      = legend_col,
         pt.cex   = 1.2,      # Slightly larger point symbols in legend
         bg       = adjustcolor("white", alpha.f = 0.95),
         box.col  = "gray40", 
         text.col = "black",
         bty      = "o",      # Box type: ordinary rectangle
         cex      = 0.85,     # Text size
         xjust    = 0,        # Align to the left (the x-value is the left edge of the legend)
         yjust    = 1)        # Top alignment (the y-value corresponds to the top edge of the legend)
}

——————————- Example of use with Crispy ——————————-

plot_all_points("ID016_Crispy", NAMES_sub, PATH)

plot_day_night_points ("ID016_Crispy", NAMES_sub, PATH)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 229.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 229.
## Using the first value.

plot_night_clusters("ID016_Crispy", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 229.
## Using the first value.
## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 229.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 7 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

## Warning: `lift()` was deprecated in purrr 1.0.0.
## ℹ The deprecated feature was likely used in the slippymath package.
##   Please report the issue at
##   <https://www.github.com/milesmcbain/slippymath/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

——————————- All maps for Aix - 2024 - Autumn ——————————-

Certains points méritent d’être examinés dans le détail (par exemple pour Roger)

plot_night_clusters("ID004_Wisk", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 73. Using
## the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 73.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 9 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("ID009_Wisk", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 260.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 260.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 5 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Fingers", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 366.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 366.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 8 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("ID007_Sleepy", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 496.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 496.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 6 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

# plot_night_clusters("ID021_Sleepy", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3) # only 1 point at crep, none during night
plot_night_clusters("Corty", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 202.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 202.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 8 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Tosca", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 222.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 222.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 8 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Rox", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 362.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 362.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 7 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Crispy", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 229.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 229.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 7 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Rouky", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 415.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 415.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 7 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Roger", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 438.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 438.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 7 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...

plot_night_clusters("Scar", NAMES_sub, PATH, distance_threshold = 40, nb_points_to_cluster = 3)

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected latDeg to be a scalar value, but was length 478.
## Using the first value.

## Warning in computeSunPosition(data_sq[[1]]$date, data_sq[[1]]$y,
## data_sq[[1]]$x): Expected longDeg to be a scalar value, but was length 478.
## Using the first value.

## Warning in Move2CSV(object, timeformat = timeformat, timezone = timezone, :
## Move objects in geographic coordinates are automatically projected.

## Minimum sampling interval of 7 minutes in 1

## Warning: Transforming 'ext' to Web Mercator (EPSG: 3857), since 'ext' has a
## different CRS. The CRS of the returned basemap will be Web Mercator, which is
## the default CRS used by the supported tile services.

## Loading basemap 'satellite' from map service 'mapbox'...