Classifying Bald Eagle Flight Behavior from High-Frequency GPS Data

Introduction

Bald eagles have large ranges and soaring flight, so they are ideal.
PCA bi plots were checked per bird to avoid over plotting.
Data collected at High-frequency GPS: every 1–11 s.
4 states: Perching, Ascending, Flapping, Gliding.
Flight often flap → gain altitude → glide.

Data

2,093,022 GPS from 57 bald eagles over ~4 years
Tagged 100 eagles (2013–2019) across MO (1), OK (14), IL (33), IA (52)

# A tibble: 6 × 15
  Animal_ID TimeDiff segment_id segment_length LocalTime           Latitude Longitude       X        Y   KPH    Sn   AGL VerticalRate abs_angle absVR
      <int>    <int>      <dbl>          <int> <dttm>                 <dbl>     <dbl>   <dbl>    <dbl> <dbl> <dbl> <dbl>        <dbl>     <dbl> <dbl>
1       105        5          1             10 2019-06-25 16:44:10     41.6     -92.9 509745. 4609596.  0.07  0.52  9.55        -0.99      0.45  0.99
2       105        6          1             10 2019-06-25 16:44:16     41.6     -92.9 509746. 4609597.  0.07  0.17  9.55         0         0.21  0   
3       105        5          1             10 2019-06-25 16:44:21     41.6     -92.9 509747. 4609597.  0.2   0.14  9.55         0         0.91  0   
4       105        6          1             10 2019-06-25 16:44:27     41.6     -92.9 509747. 4609599.  0.16  0.22  9.55         0         3.14  0   
5       105        5          1             10 2019-06-25 16:44:32     41.6     -92.9 509747. 4609598.  0.11  0.08 11.6          0.39      2.69  0.39
6       105        6          1             10 2019-06-25 16:44:38     41.6     -92.9 509746. 4609599.  0.07  0.22 11.6          0         0.53  0

What are we looking for?

Can we use the movement variables KPH, Sn, AGL, |Angle|, Vertical rate, and |VR|, to distinguish in-flight from perching points?
Of the in-flight points:

a. Are there distinct flight behaviors that can be identified?

b. What are the characteristics of these behaviors?

c. What are some visual examples of flight segments that demonstrate the different types of in-flight behaviors?

Distributions of Movement Features in Bald Eagle GPS Data

These histogram show the movement features are highly skewed

Choosing K (WSS)

Choosing K (Silhouette)

Perching vs In-Flight

In-Flight Behaviors (K = 2-7)

We compare (K = 2-7) to see how many distinct movement modes appear in PCA space

Behavioral States (K = 4)

In-Flight Behaviors: Profiles

Visual Example of In-Flight Behaviors

[1] 4625

Appendix

Data

library(tidyverse)
set.seed(42)
eagle_s <- eagle_data %>% slice_sample(n = min(20000, nrow(eagle_data)))
(eagle_s
  %>% as_tibble() 
  ) %>% head

Histogram

par(
  mfrow = c(2, 3),
  mar  = c(3.5, 3.5, 2.2, 1),
  oma  = c(2, 2.5, 0, 0)
)

h <- function(x, title, col, breaks = 60, xlim = NULL) {
  x <- x[is.finite(x)]          
  if (is.null(xlim)) {
    hist(x, freq = FALSE, breaks = breaks, col = col, border = "black",
         main = title, xlab = "", ylab = "", cex.main = 1)
  } else {
    hist(x, freq = FALSE, breaks = breaks, xlim = xlim, col = col, border = "black",
         main = title, xlab = "", ylab = "", cex.main = 1)
  }
}

h(eagle_data$KPH,        "KPH",            "darkred")
h(eagle_data$Sn,         "Sn",             "pink")
h(eagle_data$AGL,        "AGL",            "darkgreen")
h(eagle_data$abs_angle,  "|angle|",        "lightgreen")
h(eagle_data$VerticalRate, "Vertical Rate","gold")
h(eagle_data$absVR,      "|Vertical Rate|","darkorange")

mtext("Density", side = 2, outer = TRUE, line = 1)

WSS

# 20,000 rows from the whole data set
set.seed(42)
eagle_s <- eagle_data %>% slice_sample(n = min(20000, nrow(eagle_data)))

# building the feature matrix
df1 <- eagle_s %>%
  transmute(
    KPH = sqrt(KPH),
    Sn  = sqrt(Sn),
    AGL = sqrt(AGL),
    abs_angle = sqrt(abs_angle),
    VerticalRate = VerticalRate,
    absVR = sqrt(absVR)
  ) %>%
  filter(if_all(everything(), is.finite)) %>%
  scale()

# wss plot
library(factoextra)

fviz_nbclust(df1, FUNcluster = kmeans, method = "wss", k.max = 7, nstart = 20)

Silhouette

# silhouette plot
library(factoextra)

fviz_nbclust(df1, FUNcluster = kmeans, method = "silhouette", k.max = 7, nstart = 20)

Perching vs In-Flight

library(dplyr)

# 20,000 rows from the whole data set 
set.seed(42)
eagle_s <- eagle_data %>% slice_sample(n = min(20000, nrow(eagle_data)))

# building the feature matrix
X <- eagle_s %>%
  mutate(
    KPH       = sqrt(KPH),
    Sn        = sqrt(Sn),
    AGL       = sqrt(AGL),
    abs_angle = sqrt(abs_angle),
    absVR     = sqrt(absVR)
  ) %>%
  select(KPH, Sn, AGL, abs_angle, VerticalRate, absVR) %>%
  filter(if_all(everything(), is.finite))

X_scaled <- scale(X)

Seperating into 2 main groups

library(factoextra)
library(ggplot2)

# spliting into 2 main groups

pca <- prcomp(X_scaled, center = FALSE, scale. = FALSE)

set.seed(42)
km <- kmeans(X_scaled, centers = 2, nstart = 20, iter.max = 100, algorithm = "Lloyd") 
cluster <- factor(km$cluster)

perch_id <- which.min(tapply(X$KPH + X$Sn, km$cluster, median))

state <- ifelse(km$cluster == perch_id, "Perching", "In-flight")
state <- factor(state, levels = c("Perching", "In-flight"))

Biplot

# creating plot

fviz_pca_biplot( 
  pca,
  axes = c(1,2),
  habillage = state,
  geom.ind = "point",
  alpha.ind = 0.55,
  pointsize = 0.6,
  label = "var",
  repel = TRUE,
  axes.linetype = 0
) +
  scale_color_manual(values = c("Perching" = "#009E73",
                                "In-flight" = "darkorange"),
                     name = "State") +
  guides(shape = "none") +    
  labs(x = "PC1", y = "PC2", title = "PCA Biplot") +
  theme_classic() + labs(title = "PCA of Perching vs In-Flight Behaviour")

In-Flight Behaviors (K = 2-7)

#creating X_scaled

library(dplyr)

#20000 rows from the whole data set
set.seed(42)
eagle_s <- eagle_data %>% slice_sample(n = min(20000, nrow(eagle_data))) 

# building the feature matrix
X <- eagle_s %>%
  mutate(
    KPH = sqrt(KPH),
    Sn  = sqrt(Sn),
    AGL = sqrt(AGL),
    abs_angle = sqrt(abs_angle),
    absVR = sqrt(absVR)
  ) %>%
  select(KPH, Sn, AGL, abs_angle, VerticalRate, absVR) %>%
  filter(if_all(everything(), is.finite))

X_scaled <- scale(X)

Seprating perching vs In-flight

#splitting perching vs In-flight
set.seed(42)

km2 <- kmeans(X_scaled, centers = 2, nstart = 20)

perch_id <- which.min(tapply(X$KPH + X$Sn, km2$cluster, median))

state <- factor(ifelse(km2$cluster == perch_id, "Perching", "In-flight"))

X_flight <- X_scaled[state == "In-flight", , drop = FALSE]

Biplot

library(ggplot2)
library(purrr)
library(dplyr)
library(patchwork)

stopifnot(exists("X_flight"))

pca_f <- prcomp(X_flight, center = FALSE, scale. = FALSE)
scores <- as.data.frame(pca_f$x[, 1:2])
names(scores) <- c("PC1", "PC2")

scores$PC2 <- -scores$PC2

#zoom and avoid out liers

xlim <- quantile(scores$PC1, c(.01, .99), na.rm = TRUE)
ylim <- quantile(scores$PC2, c(.01, .99), na.rm = TRUE)

set.seed(42)
plots <- map(2:7, function(k) {
  km <- kmeans(X_flight, centers = k, nstart = 20)
  df <- transform(scores, cluster = factor(km$cluster))

  ggplot(df, aes(PC1, PC2, color = cluster)) +
    geom_point(size = 0.3, alpha = 0.35) +
    scale_color_manual(values = c(
  "1" = "#009E73",
  "2" = "gold",
  "3" = "dodgerblue",
  "4" = "darkorange2",
  "5" = "grey3",
  "6" = "orchid3",
  "7" = "tomato"
)) +
    coord_cartesian(xlim = xlim, ylim = ylim) +
    labs(title = paste0("K=", k), x = "PC1", y = "PC2") +
    theme_bw() +
    theme(
      panel.grid = element_blank(),
      legend.position = "none",
      plot.title = element_text(hjust = 0.5, face = "bold"),
      panel.border = element_rect(color = "black", fill = NA, linewidth = 0.8),
      axis.text = element_text(color = "black"),
      axis.ticks = element_line(color = "black")
    )
})

wrap_plots(plots, ncol = 3)

Behavioral States (K = 4)

library(dplyr)
library(ggplot2)

#20000 rows from the whole data set
set.seed(42)
eagle_s <- eagle_data %>% slice_sample(n = min(20000, nrow(eagle_data)))
 
# building the feature matrix
X_trans <- eagle_s %>%
  transmute(
    KPH          = sqrt(KPH),
    Sn           = sqrt(Sn),
    AGL          = sqrt(AGL),
    abs_angle    = sqrt(abs_angle),
    VerticalRate = VerticalRate,
    absVR        = sqrt(absVR)
  ) %>%
  filter(if_all(everything(), is.finite))

X_scaled <- scale(X_trans)

Runs k-means and adds the cluster label

set.seed(42)
km4 <- kmeans(X_scaled, centers = 4, nstart = 50)

df <- as.data.frame(X_trans) %>% mutate(cluster = factor(km4$cluster))

Summarizing each cluster and relabeling into behaviors

#summarizing each cluster and relabeling into behaviors
summ <- df %>%
  group_by(cluster) %>%
  summarise(
    medKPH   = median(KPH),
    medSn    = median(Sn),
    medVR    = median(VerticalRate),
    .groups  = "drop"
  )

perch_id <- summ$cluster[which.min(summ$medKPH + summ$medSn)]
asc_id   <- summ %>% filter(cluster != perch_id) %>% slice_max(medVR, n = 1, with_ties = FALSE) %>% pull(cluster)
rem      <- setdiff(summ$cluster, c(perch_id, asc_id))
glide_id <- summ %>% filter(cluster %in% rem) %>% slice_max(medKPH + medSn, n = 1, with_ties = FALSE) %>% pull(cluster)
flap_id  <- setdiff(rem, glide_id)

label_map <- setNames(
  c("Perching","Ascending","Flapping","Gliding"),
  c(as.character(perch_id), as.character(asc_id), as.character(flap_id), as.character(glide_id))
)

# creating behavior factor
df$behavior <- factor(label_map[as.character(df$cluster)],
                      levels = c("Perching","Ascending","Flapping","Gliding"))

Extracting and Flipping

pca <- prcomp(X_trans, center = TRUE, scale. = TRUE)

#extracting PC1/PC2 scores
scores <- as.data.frame(pca$x[, 1:2])
names(scores) <- c("PC1","PC2")
scores$behavior <- df$behavior

#extracting variable scores
load <- as.data.frame(pca$rotation[, 1:2])
names(load) <- c("PC1","PC2")
load$var <- rownames(load)

#flipping PC axes
if (load$PC1[load$var == "KPH"] < 0) { scores$PC1 <- -scores$PC1; load$PC1 <- -load$PC1 }
if (load$PC2[load$var == "VerticalRate"] < 0) { scores$PC2 <- -scores$PC2; load$PC2 <- -load$PC2 }

Setting, Scaling and relabling

#setting plot limits
xlim <- quantile(scores$PC1, c(.01, .99))
ylim <- quantile(scores$PC2, c(.01, .99))

#scaling the loading arrows
arrow_scale <- 0.9 * min(diff(xlim), diff(ylim))
load$PC1a <- load$PC1 * arrow_scale
load$PC2a <- load$PC2 * arrow_scale

pretty <- c(abs_angle="|angle|", VerticalRate="VR", absVR="|VR|")
load$lab <- ifelse(load$var %in% names(pretty), pretty[load$var], load$var)

Find where each behavior cluster sits

cent <- scores %>%
  group_by(behavior) %>%
  summarise(cx = median(PC1), .groups = "drop") %>%
  arrange(cx)

cent$cy <- ylim[1] - 0.12 * diff(ylim)

x_targets <- seq(xlim[1] + 0.15*diff(xlim),
                 xlim[2] - 0.15*diff(xlim),
                 length.out = nrow(cent))
cent$cx_lab <- x_targets

Biplot

ggplot(scores, aes(PC1, PC2, color = behavior)) +
  geom_point(size = 0.35, alpha = 0.28) +
  scale_color_manual(values = c(
    "Perching"  = "#009E73",
    "Ascending" = "gold",
    "Flapping"  = "darkorange2",
    "Gliding"   = "dodgerblue"
  )) +
  geom_segment(
    data = load,
    aes(x = 0, y = 0, xend = PC1a, yend = PC2a),
    inherit.aes = FALSE,
    color = "steelblue",
    linewidth = 0.7,
    arrow = arrow(length = grid::unit(0.25, "cm"))
  ) +
  geom_text(
    data = load,
    aes(x = PC1a, y = PC2a, label = lab),
    inherit.aes = FALSE,
    size = 4,
    fontface = "bold",
    hjust = ifelse(load$PC1a >= 0, -0.1, 1.1),
    vjust = ifelse(load$PC2a >= 0, -0.2, 1.2)
  ) +
  geom_text(
  data = cent,
  aes(x = cx_lab, y = cy, label = behavior, color = behavior),
  inherit.aes = FALSE,
  fontface = "bold",
  size = 5
)+
  coord_cartesian(xlim = xlim, ylim = c(cent$cy[1], ylim[2])) +
  theme_classic(base_size = 18) +
  theme(
    legend.position = "none",
    panel.border = element_rect(color="black", fill=NA, linewidth=0.8)
  ) +
  labs(title = "K = 4", x = "PC1", y = "PC2")

In-Flight Behaviors: Profiles

library(dplyr)
library(tidyr)
library(ggplot2)

set.seed(42)

eagle_s <- eagle_data %>% slice_sample(n = min(20000, nrow(eagle_data)))

X <- eagle_s %>%
  mutate(
    KPH       = sqrt(KPH),
    Sn        = sqrt(Sn),
    AGL       = sqrt(AGL),
    abs_angle = sqrt(abs_angle),
    absVR     = sqrt(absVR)
  ) %>%
  select(KPH, Sn, AGL, abs_angle, VerticalRate, absVR) %>%
  filter(if_all(everything(), is.finite))

X_scaled <- scale(X)

Splitting and finding behaviours

km2 <- kmeans(X_scaled, centers = 2, nstart = 20)
perch_id <- which.min(tapply(X$KPH + X$Sn, km2$cluster, median, na.rm = TRUE))
state <- ifelse(km2$cluster == perch_id, "Perching", "In-flight")

X_flight     <- X_scaled[state == "In-flight", , drop = FALSE]  
X_flight_raw <- X[state == "In-flight", , drop = FALSE]

Runs k-means and adds the cluster label

k_beh <- 3
set.seed(42)
km_beh <- kmeans(X_flight, centers = k_beh, nstart = 50)

df_flight <- as.data.frame(X_flight_raw) %>%
  mutate(cluster = factor(km_beh$cluster))

Summarise

summ <- df_flight %>%
  group_by(cluster) %>%
  summarise(
    med_VR  = median(VerticalRate, na.rm = TRUE),
    med_KPH = median(KPH, na.rm = TRUE),
    .groups = "drop"
  )

asc_id   <- summ$cluster[which.max(summ$med_VR)]
rem      <- setdiff(summ$cluster, asc_id)
glide_id <- summ %>% filter(cluster %in% rem) %>% slice_max(med_KPH, n = 1, with_ties = FALSE) %>% pull(cluster)
flap_id  <- setdiff(rem, glide_id)

map <- setNames(c("Ascending","Flapping","Gliding"),
                c(as.character(asc_id), as.character(flap_id), as.character(glide_id)))

df_flight <- df_flight %>%
  mutate(
    behavior = factor(map[as.character(cluster)],
                      levels = c("Ascending","Flapping","Gliding"))
  )

pretty <- c(abs_angle="|angle|", VerticalRate="VR", absVR="|VR|")
load$lab <- ifelse(load$var %in% names(pretty), pretty[load$var], load$var)

Boxplot

long <- df_flight %>%
  pivot_longer(
    cols = c(abs_angle, absVR, AGL, KPH, Sn, VerticalRate),
    names_to = "variable",
    values_to = "value"
  )

ggplot(long, aes(behavior, value, fill = behavior)) +
  geom_boxplot(width = 0.6, outlier.shape = NA) +
  scale_fill_manual(values = c(
    "Ascending" = "gold",
    "Flapping"  = "darkorange2",
    "Gliding"   = "dodgerblue"
  )) +
  facet_wrap(~variable, scales = "free_y", ncol = 3) +
  theme_classic() +
  theme(legend.position = "none") +
  labs(x = "In-flight behavior", y = "Value")

Visual Example of In-Flight Behaviors

library(dplyr)
library(ggplot2)
library(patchwork)
 
vars <- c("KPH", "Sn", "AGL", "VerticalRate")

df <- eagle_data %>%
drop_na(all_of(vars))
 
X <- scale(df %>% select(all_of(vars)))
 
set.seed(123)
km <- kmeans(X, centers = 4)

df$cluster <- factor(km$cluster)

perching_cluster <- df %>%
group_by(cluster) %>%
summarise(mean_speed = mean(KPH, na.rm = TRUE)) %>%
arrange(mean_speed) %>%
slice(1) %>% pull(cluster)

bird_id <- 107
 
flight_segment <- df %>%
  filter(Animal_ID == bird_id, cluster != perching_cluster) %>%
  group_by(segment_id) %>%
  summarise(kph_range = max(KPH) - min(KPH), .groups="drop") %>%
  arrange(desc(kph_range)) %>%
  slice(1) %>%
  pull(segment_id)
 
flight_segment   
 
bird_df <- df %>%
filter(Animal_ID == bird_id) %>%
filter(segment_id == flight_segment) %>% 
filter(cluster != perching_cluster) %>% 
arrange(LocalTime) %>%
mutate(seconds = cumsum(TimeDiff))
 
remaining_clusters <- sort(unique(bird_df$cluster))
 
custom_colors <- setNames(
c("gold", "darkorange2", "dodgerblue"),
remaining_clusters)

Plot KPH

# Plot KPH
p_kph <- ggplot(bird_df, aes(x = seconds, y = KPH, color = cluster)) +
geom_line(alpha = 0.5) +
geom_point(size = 1.2) +
scale_color_manual(values = custom_colors) +
theme_minimal(base_size = 14) +
labs(title = paste("Speed (KPH)", "Segment", flight_segment),
x = "Seconds",
y = "KPH" )

Plot AGL

# Plot AGL
p_agl <- ggplot(bird_df, aes(x = seconds, y = AGL, color = cluster)) +
geom_line(alpha = 0.5) +
geom_point(size = 1.2) +
scale_color_manual(values = custom_colors) +
theme_minimal(base_size = 14) +
labs(title = paste("Altitude (AGL)", "Segment", flight_segment),
x = "Seconds",
y = "Altitude (m)")
 
final_plot <- p_kph / p_agl
final_plot