Given a flight trajectory [(t, x, y)] and co-located wind vectors [(u, v)], calculate the relative angle between flight and wind along the trajectory

relative_wind.R

library(dplyr)
library(geosphere)
library(glue)
library(ggplot2)
library(mapproj)

# Example data frame
set.seed(1)
tracks <- tibble(
  time = 1:5,
  lon = -122 + cumsum(rnorm(5)),
  lat = 37 + cumsum(rnorm(5)),
  u = cumsum(rnorm(5)),
  v = cumsum(rnorm(5))
)

# Visualize tracks and wind vectors. Wind vectors in blue
(p <- ggplot(tracks) +
  geom_segment(aes(x = lon, y = lat, xend = lead(lon), yend = lead(lat)),
               arrow = arrow(angle = 15, 
                             length = unit(10, "points"), 
                             type = "closed")) +
  geom_segment(aes(x = lon, y = lat,
                   xend = lon + u / 4, yend = lat + v / 4),
               arrow = arrow(angle = 15, 
                             length = unit(10, "points"), 
                             type = "closed"),
               color = "blue") +
  geom_label(aes(x = lon, y = lat, label = time), nudge_x = -0.15) +
  coord_map() +
  theme_classic())

# Calculate wind direction relative to travel
relative_direction <- function(x, y, u, v) {
  # Divide the uv vector by its norm to get the unit vector along uv
  uv <- cbind(u, v)
  norm_uv <- apply(uv, 1, function(row) sqrt(sum(row^2)))
  uv_unit <- uv / cbind(norm_uv, norm_uv)
  # Track unit vectors calculated from sequential bearing
  # Bearing by default is in "degrees east of north", but we need "radians north 
    # of east" to compare against uv
  track_bearing <- (90 - bearing(cbind(x, y))) * pi / 180
  track_unit <- cbind(cos(track_bearing), sin(track_bearing))
  # Angle between track and wind is the difference of the arctangents
  atan2(track_unit[,1], track_unit[,2]) - atan2(uv_unit[,1], uv_unit[,2])
}
annotated_tracks <- tracks %>% 
  mutate(winddir = relative_direction(lon, lat, u, v))

# Show plot with wind angle
p + geom_label(aes(x = lon, y = lat,
                   label = glue("{format(winddir*180/pi, digits=2)}º")),
               annotated_tracks,
               nudge_x = 0.3, nudge_y = 0.1)

# Sanity checks
# A track and wind in the same direction should get 0 rad
relative_direction(x = c(0, 1), y = c(0, 0),
                   u = c(1, 1), v = c(0, 0))
# A track and wind in opposite directions should get pi rad
relative_direction(x = c(0, 1), y = c(0, 0),
                   u = c(-1, -1), v = c(0, 0))
# Wind *from* the left should get -pi/2 rad
relative_direction(x = c(0, 0), y = c(0, 1),
                   u = c(1, 1), v = c(0, 0))
# Wind *from* the right should get pi/2 rad
relative_direction(x = c(0, 0), y = c(0, 1),
                   u = c(-1, -1), v = c(0, 0))