Add files via upload

Author: Observarun

01_step-length_turning-angle_ANIMOVER1.R

@@ -0,0 +1,38 @@
+# Code for preprocessing ANIMOV1 simulation data.
+
+
+df_DARs_synth <- read.csv(file=".../2KernelMovement.csv", skip=6)  # first few rows have explanatory text, hence skipped.
+df_DARs_synth <- na.omit(df_DARs_synth)
+
+
+# To bring time in order
+
+df_DARs_synth <- transform(df_DARs_synth, Delta=Delta+99*(Day-1))  # dataframe for synthetic (ANIMOVER1) data; time parameter in ANIMOV1 resets at the turn of each day.
+
+
+# To calculate turning angle given angle of heading.
+
+colnames(df_DARs_synth)[6] <- "Heading"  # Changing name of column 6 to reflect angle of heading.
+df_DARs_synth$Heading <- with(df_DARs_synth, pmin(df_DARs_synth$Heading, 360-df_DARs_synth$Heading))
+
+Turning_angle_deg <- with(df_DARs_synth, diff(df_DARs_synth$Heading))  # Turning angle
+Turning_angle_deg <- abs(Turning_angle_deg)  # To ensure turning angle is in (0,pi]
+df_DARs_synth <- df_DARs_synth[-1, ]  # time parameters in the ANIMOVER\_1 RAMP to produce $99$ points # for each of our "nominal days" other than for the first, which has $100$ points.
+df_DARs_synth$Heading <- Turning_angle_deg
+rm(Turning_angle_deg)
+colnames(df_DARs_synth)[6] <- "Turning_angle_deg"
+df_DARs_synth <- transform(df_DARs_synth, Turning_angle_deg=Turning_angle_deg*pi/180)
+colnames(df_DARs_synth)[6] <- "Turning_angle_rad"
+
+
+# Renaming variables
+
+df_DARs_synth <- df_DARs_synth[, c('X', 'Y', 'Delta', 'Distance', 'Turning_angle_rad')]
+colnames(df_DARs_synth)[1] <- "x"
+colnames(df_DARs_synth)[2] <- "y"
+colnames(df_DARs_synth)[3] <- "delta"
+colnames(df_DARs_synth)[4] <- "speed"
+colnames(df_DARs_synth)[5] <- "turning_angle"
+
+
+saveRDS(df_DARs_synth, file="merged-DARs_ANIMOV.Rds")  # dataframe with step length and turning angle variables

01_step-length_turning-angle_barn-owl.R

@@ -0,0 +1,39 @@
+# Code for preprocessing a single barn owl individual data.
+
+
+list_DARs_indiv <- readRDS(".../mid_to_mid1.Rds")  # see https://github.com/LudovicaLV/DAR_project/blob/main/DAR1_measures.R for the code to extract DARs from multi-DAR time-series data of an individual.
+
+df_DARs_indiv <- do.call(rbind, list_DARs_indiv)  # working w/ all DARs of an individual
+
+
+# Calculate the variables for speed and turning angle, and append to the data frame.
+
+df_DARs_indiv <- df_DARs_indiv[, c('X','Y','dateTime')]
+colnames(df_DARs_indiv)[1] <- "x"
+colnames(df_DARs_indiv)[2] <- "y"
+colnames(df_DARs_indiv)[3] <- "t"
+
+df_DARs_indiv <- transform(df_DARs_indiv, T=as.numeric(T))
+speed <- vector("double", nrow(df_DARs_indiv))
+for (i in 2:nrow(df_DARs_indiv))
+{
+  speed[i] <- sqrt(sum((df_DARs_indiv[i,] - df_DARs_indiv[i-1,])^2)
+                - (df_DARs_indiv[i,3] - df_DARs_indiv[i-1,3])^2)/(df_DARs_indiv[i,3] - df_DARs_indiv[i-1,3])
+}
+df_DARs_indiv$speed <- speed
+df_DARs_indiv <- transform(df_DARs_indiv, T=as.POSIXct(T, origin="1970-01-01", tz="GMT"))
+
+turning_angle <- vector("double", nrow(df_DARs_indiv))
+for (i in 3:nrow(df_DARs_indiv))
+{
+  turning_angle[i] <- abs(atan((df_DARs_indiv[i,2]-df_DARs_indiv[i-1,2])/(df_DARs_indiv[i,1]-df_DARs_indiv[i-1,1]))
+               - atan((df_DARs_indiv[i-1,2]-df_DARs_indiv[i-2,2])/(df_DARs_indiv[i-1,1]-df_DARs_indiv[i-2,1])))
+}
+df_DARs_indiv$turning_angle <- turning_angle
+
+
+df_DARs_indiv <- df_DARs_indiv[-1, ]  # 1st point discarded for lack of turning angle value
+
+
+saveRDS(df_DARs_indiv, file="merged-DARs_barn-owl.Rds")
+

02_segmentation.R

@@ -0,0 +1,53 @@
+# Code for segmentation and clustering of either a barn owl individual or
+# ANIMOVER_1 synthetic data
+
+
+library(dplyr)
+library(tseries)
+
+
+df_DARs <- data.frame(readRDS(".../merged-DARs_ANIMOV.Rds"))  # for ANIMOV data; line 9 must be uncommented and line 10 commented for ANIMOV data
+#df_DARs <- data.frame(readRDS(".../merged-DARs_barn-owl.Rds"))  # for barn owl data; line 9 must be commented and line 10 uncommented for barn owl data
+
+
+df_DARs <- df_DARs[df_DARs_indiv$speed < 70,]  # to remove obviously bogus points (unrealistically high barn owl speeds)
+
+df_DARs$speed <- (df_DARs$speed)/max(df_DARs_indiv$speed)
+df_DARs$turning_angle <- (df_DARs$turning_angle)/pi
+# Scaling to bring the variables to range on [0,1].
+
+
+# Method to create segments and construct a representation.
+
+create_chunks <- function(df)
+{
+  df_chunks_DARs <- df |>
+    group_by(seg=cut(T, breaks="60 sec")) |>  # for barn owl; line 25 must be uncommented and line 26 commented for barn owl data 
+    #group_by(seg=cut(T, breaks=seq(0, max(T), by=30))) |>  # for ANIMOV; line 25 must be commented and line 26 uncommented for ANIMOV data 
+    summarise(
+      speed_mean = mean(speed, na.rm = TRUE),
+      speed_std = sd(speed, na.rm = TRUE),
+      turning_angle_mean = mean(turning_angle, na.rm = TRUE),
+      turning_angle_std = sd(turning_angle, na.rm = TRUE),
+      disp_ends = sqrt((x[n()] - x[1])^2 + (y[n()] - y[1])^2),
+      disp_consec = sum(sqrt(diff(x)^2 + diff(y)^2)),
+      n_pts = n(),
+      fraction_WP = sum(grepl("^WP", Kernel))/n(),  # Only for ANIMOV data; comment line 35 and comma at the end of line 34 for barn owl data
+    ) |>
+    filter(n_pts >= 12) |>  # Only for the case of barn owls
+    summarise(
+      seg,
+      speed_mean, speed_std, turning_angle_mean, turning_angle_std,
+      net_displacement = disp_ends/disp_consec,
+      fraction_WP  # Only for ANIMOV data; comment line 42 and comma at the end of line 41 for barn owl data
+    )
+  
+  df_chunks_DARs <- na.omit(df_chunks_DARs)
+  
+  return(df_chunks_DARs)
+}
+
+df_chunks_DARs = create_chunks(df_DARs)
+
+
+saveRDS(df_chunks_DARs, file="segmented-DARs.Rds")

03_clustering_visualisation.R

@@ -0,0 +1,100 @@
+library(parallelDist)
+library(fastcluster)
+library(factoextra)
+
+
+df_chunks_DARs <- data.frame(readRDS(".../merged-DARs_ANIMOV.Rds"))
+
+
+# Hierarchical agglomerative clustering w/ Ward's criterion.
+
+d <- parDist(as.matrix(df_chunks_DARs[c("speed_mean","turning_angle_mean","speed_std","turning_angle_std","net_displacement")]),
+             method = "euclidean")
+hierarc_res <- fastcluster::hclust(d, method="ward.D2")
+
+hc_ward_8 <- cutree(hierarc_res, k=8)
+
+centroids <- NULL
+for(k in 1:8){
+  centroids <-
+    rbind(
+      centroids,
+      colMeans(df_chunks_DARs[c("speed_mean","turning_angle_mean","speed_std","turning_angle_std","net_displacement")][hc_ward_8==k, , drop = FALSE])
+    )
+}
+print(noquote(paste0("Centroids of clusters: ",centroids)))
+
+occupancy <-
+  as.numeric(
+    table(hc_ward_8)
+  )
+print(noquote(paste0("Occupancy of clusters: ",occupancy)))
+
+df_chunks_DARs_synth <- cbind(df_chunks_DARs_synth, clustNum=hc_ward_8)
+df_chunks_DARs_synth |>
+  group_by(clustNum) |>
+  summarise(
+    Cluster_avg_fraction_WP_points_in_a_segment=mean(fraction_WP),
+    Cluster_std_fraction_WP_points_in_a_segment=sd(fraction_WP),
+    Cluster_std_net_displacement_of_a_segment=sd(net_displacement)
+    )  # lines 34 to 40 to be commented out for working with barn owl data
+
+
+# Compute principal components.
+
+pca_res <-
+  prcomp(
+    df_chunks_DARs[, c("speed_mean","turning_angle_mean","speed_std","turning_angle_std","net_displacement")],
+    scale=FALSE,
+    center = FALSE
+  )
+
+# Eigenvalues
+get_eigenvalue(pca_res)
+
+fviz_pca_var(pca_res,
+             col.var = "contrib",  # color by contributions to the PC
+             repel = TRUE,  # avoid text overlapping
+             scale = FALSE,
+             center = FALSE
+)
+
+pca_res$x
+
+# eigenvectors
+pca_res$rotation
+
+
+# Plot StaMEs in mean speed - mean turning angle space or PC1-PC2 space w/ different colours representing clusters.
+
+cluster_data_PC <-
+  data.frame(
+    Cluster = as.factor(hc_ward_8),
+    PC1 = pca_res$x[, 1],
+    PC2 = pca_res$x[, 2]
+  )
+var_expl <- round(pca_res$sdev^2 / sum(pca_res$sdev^2) * 100, 2)
+# Create a cluster plot
+ggplot(cluster_data_PC, aes(x = PC1, y = PC2, colour = Cluster)) +
+  geom_point(size = 3) +
+  scale_colour_manual(values = c("green1", "blue4", "lightblue1", "violet", "yellow1", "red1", "orange1", "red4")) +
+  labs(title = "Cluster Plot on Principal Components",
+       x = paste("PC1 (", var_expl[1], "%)", sep = ""),
+       y = paste("PC2 (", var_expl[2], "%)", sep = "")) +
+  theme_minimal()
+
+cluster_data_SpTa <- data.frame(
+  Cluster = as.factor(hc_ward_8),
+  mean_speed = df_chunks_DARs_$speed_mean,
+  mean_turning_angle = df_chunks_DARs$turning_angle_mean
+)
+# Create a cluster plot
+ggplot(
+    cluster_data_SpTa,
+    aes(x = mean_speed, y = mean_turning_angle, colour = Cluster)
+) +
+geom_point(size = 3) +
+  #scale_color_manual(values = "red4")
+  scale_color_manual(values = c("green1", "blue4", "lightblue1", "violet", "yellow1", "red1", "orange1", "red4")) +
+  labs(title = "Cluster Plot on mean speed and mean turning angle") +
+  theme_minimal()

README.md

@@ -0,0 +1,23 @@
+# The Statistical Building Blocks of Animal Movement Simulations
+
+
+This repository shows a procedure to construct Statistical Movement Elements (StaMEs) from animal movement track. Movement data itself is in the form of a relocation time-sequence, given by a walk $`{W} = \{(t;x(t),y(t))|t=0,...,n^{time} \}.`$ The algorithm parses the normalised velocity (equivalent to normalised displacement or step length) and turning-angle time series (derived from the walk) into $z=1,\cdots,n^{\rm seg}$ segments (where $` n^{\rm seg}=\lfloor \frac{n^{\rm time}-1}{\mu} \rfloor `$) of size $\mu$ each. A clustering procedure is then applied on a hand-engineered representation comprising a set ($\cal S_{\mu}$) of the following statistical variables for each segment: mean velocity $V$, mean absolute turning angle $| \Delta \Theta |$, the associated standard deviations ${\rm SD}^{V}$ and ${\rm SD}^{| \Delta \Theta |}$, and finally, a normalized net displacement ($\Delta^{\rho}$). The last one is the Euclidean distance between the end points of each segment divided by a product of the number of points and mean step-length, and is required to pick up any possible circular motion type biases in movement.
+
+Clustering is an unsupervised machine learning procedure which finds an optimum partition of a set of points with two or more variables into subsets (or clusters) of similar points. Here I've used hierarchical clustering approach, which results in a tree structure called dendrogram having a single cluster at the root, with leaf nodes representing data points in the set. Cluster analysis is performed on the segment level data (i.e., with $n_{\rm seg} = \lfloor \frac{t}{\mu} \rfloor +1$ points) using the variables in the se $\cal S_{\mu}$. Specifically, we want to find an optimal partition $\cal P(\cal S_{\mu}) = \{C_1, C_2, ..., C_k\},$ where $\cup_{i=1}^k C_i = \cal S_{\mu}$, and $C_i \cap C_j = \emptyset$ (hard clustering).
+
+Hierarchical agglomerative algorithms implement bottom-up clustering methodology, which starts with each point of $\cal S$ in its own singleton cluster, followed by successive fusions of two clusters at a time, depending on similarity between them, leading to a specified number of clusters (or, alternatively, leading to one cluster followed by cutting the dendrogram at the desired number of clusters). These are deterministic, yet greedy, in the sense that the clusters are merged based entirely on the similarity measure, thereby yielding a local solution. Most of the similarity schemes are specified not in terms of an objective function to be optimized, but procedurally. I use Ward's minimum sum of square scheme, which performs a fusion of clusters while minimizing the intra-cluster variance. The distance metric to quantify dissimilarity is the Euclidean measure.
+
+## Description of scripts
+
+* 01_step-length_turning-angle_barn-owl.R: Computation of step length (or speed) and turning angle time series starting from the multi-DAR relocation time-series data for a barn owl individual obtained using the ATLAS reverse GPS system in north-eastern Israel.
+* 01_step-length_turning-angle_ANIMOVER1.R: Computation of step length and turning angle time series starting from 2-mode simulated multi-DAR relocation data generated from Numerus ANIMOVER_1.
+* 02_segmentation.R: Parsing of multi-DAR speed and turning angle series into segments of $`\mu`$ points each to construct a representation using a set of statistics for each segment for barn owl and ANIMOV data.
+* 03_clustering_visualisation.R: Clustering analysis on the representation previously obtained to perform StaME extraction along with visualisation of results for barn owl and ANIMOV data.
+
+## Results
+
+![ClustersOwl](https://github.com/Observarun/StaME-extraction_CAMs/assets/83636458/eb8c67a1-90ba-4407-8e61-866430789d92)
+*StaMEs obtained by clustering segments from the tracks of two different barn owls.*
+
+![ClustersSim](https://github.com/Observarun/StaME-extraction_CAMs/assets/83636458/dfb5fef1-2af4-401e-b99d-c779e61516e7)
+*Clustering results from 10-point and 30-point segmentation of the simulation data.*