In mpio-be/windR: Analysis of tracking data in relation to wind

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

This file describes how to connect tracking data with wind data, calculate the wind support and cross wind, plot these data with a map and as a comet plot animation.

Summary

Connect tracks with wind data
Calculate wind support & cross wind
Plot tracks colored with wind support
Comet plot animation of tracks

Load packages, tracking data and set working directory

# Packages
sapply(c('magrittr', 'data.table', 'foreach', 'windR', 'raster', 'foreach', 'doParallel', 'ggplot2', 'stringr', 'rgeos'),
       function(x) suppressPackageStartupMessages(require(x , character.only = TRUE, quietly = TRUE) ) )

# Projection (polar Lambert azimuthal equal-area with longitude origin 156.65° W)
PROJ   = '+proj=laea +lat_0=90 +lon_0=-156.653428 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 '

# load data
d = PESA_example_tracks
d[, datetime_ := as.POSIXct(datetime_)]

# Specify directory
wd = paste0(getwd(), '/temp')
# wd = tempdir() # temporary in this case (choose yourself)

1. Connect tracks with wind data

# load wind data
wind_data = system.file('ERA_Interrim', 'ERA_Interrim_850mb_4_7_June_2014_10km.RDS', package = 'windR')
w = readRDS(wind_data)

# select layer of closet wind data
w_datetime_ = unique(w$datetime_)
d[, w_date := closestDatetime(datetime_, w_datetime_), by = 1:nrow(d)]
d[, id     := 1:.N]

# should not be over 3 hours difference (wind data are in 6-hour steps)
hist(x  = as.numeric(difftime(d$datetime_,  d$w_date, units = "hours")))

# get wind data for each point
setkey(d, id)
setkey(w, datetime_)

d[, c('u', 'v') := getWind(x, y, w = w[J(w_date), nomatch=0L], PROJ), by = id]

2. Calculate wind support & cross wind

# distance to next point
d[, x2   := data.table::shift(x, type = 'lead'), by = bird_patch_id]
d[, y2   := data.table::shift(y, type = 'lead'), by = bird_patch_id]
d[, dist := sqrt(sum((c(x, y) - c(x2, y2))^2)) , by = 1:nrow(d)]

# time to the next point
d[, datetime_2  := data.table::shift(datetime_, type = 'lead'), by = bird_patch_id]
d[, TbtwPoints  := as.numeric(difftime(datetime_2, datetime_, units = 'sec'))]

# ground speed & direction
d[, g_speed := dist/TbtwPoints]
d[, g_direction := bearing(x , y, x2, y2), by = 1:nrow(d)]

# Calculate wind direction & speed
d[, w_direction := atan2(u, v)]
d[, w_speed     := sqrt(u^2 + v^2)]

# Calculate the wind support
d[, Ws := as.numeric(windSupport(g_direction = g_direction, g_speed = g_speed, w_direction = w_direction, w_speed = w_speed, crosswind = FALSE)),
  by = row.names(d)]

# Calculate the relative wind support
d[, rWs := Ws / w_speed, by = 1:nrow(d)]

# Calculate the cross wind
d[, Wc := as.numeric(windSupport(g_direction = g_direction, g_speed = g_speed, w_direction = w_direction, w_speed = w_speed, crosswind = TRUE)),
  by = row.names(d)]

3. Plot tracks colored with wind support

# Get borders of the wind data
w1 = w[datetime_ == unique(w$datetime_[1])]
wp = SpatialPoints(cbind(w1$x, w1$y), proj4string = CRS(PROJ))
datBdry = gEnvelope(wp)

# load high resolution map of area around Barrow
map_data = system.file('Map', 'BRW_map.RDS', package = 'windR')
BRW = readRDS(map_data)

# Create simple base map
b  = datBdry %>% fortify %>% data.table
bb = gBuffer(datBdry, width = 30000) %>% gEnvelope %>% fortify %>% data.table
wr = BRW %>% fortify %>% data.table

bm =
  ggplot() + coord_equal() +
  labs(x = NULL, y = NULL) +
  geom_polygon(  data = bb, aes(long, lat), fill = NA, colour = 1, size = .1) +
  geom_polygon(  data = b,  aes(long, lat), fill = 'dodgerblue4', colour = 1, size = .1) +
  geom_polygon(  data = wr, aes(long, lat, group = group), fill   = 'steelblue4' , colour = 'white', size = .2) +
  # ggsn::scalebar(data = b,  dist = 100, st.size = 3 ,height = .01) +
  rangeMapper::theme_rangemap()

print(bm)


### add tracks coloured with wind support
# Wind support scale
col_Ws = c('firebrick4', 'firebrick3', 'gold', 'springgreen3', 'springgreen4', 'springgreen4')

Ws   = d$Ws %>% na.omit(.)
Ws_n = (-6 - min(Ws)) / (max(Ws) - min(Ws))
Ws_0 = (0  - min(Ws)) / (max(Ws) - min(Ws))
Ws_p = (6  - min(Ws)) / (max(Ws) - min(Ws))
Ws_m = (max(Ws)  - min(Ws)) / (max(Ws) - min(Ws))

col_WS_v  = c(0, Ws_n, Ws_0, Ws_p, Ws_m, 1)

# plots
gi =
  bm +
  geom_path(data = d, aes(x, y, group = bird_patch_id , colour = Ws), lineend = "round", size = 0.7) +
  scale_colour_gradientn(colours = col_Ws, values = col_WS_v, limits=c(min(d$Ws, na.rm = T), max(d$Ws, na.rm = T)), name = 'Wind support (m/s)') +
  theme(legend.direction = 'horizontal', legend.position = c(0.12, 0.13),
        legend.title = element_text(face = 'bold', color = 'white'),
        legend.text = element_text(face = 'bold', color = 'white')) +
  guides(colour = guide_colourbar(title.position = 'top', title.hjust = 0.5, barwidth = 8))

gi


# Function to extract legend (from: https://gist.github.com/crsh/be88be19233f1df4542aca900501f0fb#file-gglegend-r-L7)
# otherwise plots without tracks would have no legend

gglegend <- function(x){
  tmp <- ggplot_gtable(ggplot_build(x))
  leg <- which(sapply(tmp$grobs, function(y) y$name) == "guide-box")
  tmp$grobs[[leg]]
}


legend = gglegend(gi)
# legend position
legend$vp$x <- unit(.13, 'npc')
legend$vp$y <- unit(.08, 'npc')

4. Comet plot animation of tracks

# Set time frame and path
ts = data.table(date = seq('2014-06-05' %>% as.Date %>% as.POSIXct, '2014-06-07' %>% as.Date %>% as.POSIXct, by = '30 mins') )
setkey(ts, date)
tmp_path = wd # change path
ts[, path := paste0(tmp_path, '/', str_pad(1:.N, 4, 'left', pad = '0'), '.png')   ]

# calculate last datetime of each point (to make it disappear at this point)
d[, lastDate := max(datetime_, na.rm = TRUE), by = bird_patch_id]
setkey(d,  datetime_)


# register parallel computing
cl = 20 %>% makePSOCKcluster; registerDoParallel(cl)


# loop that creates pictures for the animation
foreach(i = 1:nrow(ts), .packages = c('scales', 'ggplot2', 'lubridate', 'stringr', 'data.table', 'windR', 'grid') ) %dopar% {

  png(filename = ts[i, path], width = 700, height = 700, units = "px", pointsize = 9, bg = "white")

  tmp_date = ts[i]$date   # current date

  # add time stamp to base map
  bm_w = bm +
    scale_colour_gradientn(colours = col_Ws, values = col_WS_v, limits=c(min(d$Ws, na.rm = T), max(d$Ws, na.rm = T)), name = 'Wind support (m/s)') +
    annotate('text', x = 0, y = -1730000,
             label = paste0('2014 ', format(tmp_date, "%B %d %H:00")),
             color = 'white', size = 8, fontface = 1) + # size = 5 normal resolution, 8 for high
    theme(legend.position = 'none',
          plot.margin = unit(c(-1, -1, -1, -1), 'cm'))

  # subset bird tracks
  ki = d[datetime_ <= tmp_date & tmp_date < lastDate]

  if (nrow(ki) > 0) ki[, a:=  alphaAlong(datetime_, head = 70, skew = -0.1) ,     by = bird_patch_id] # alpha
  if (nrow(ki) > 0) ki[, s:=  sizeAlong( datetime_, head = 70, to = c(0.1,2.5)) , by = bird_patch_id] # size

  gi =
    bm_w +
    geom_path(data = ki, aes(x, y, group = bird_patch_id , colour = Ws), lineend = "round",  alpha = ki$a, size = ki$s)


  if (nrow(ki) == 0) {
    print(bm_w); grid.draw(legend)
  } else {
    print(gi); grid.draw(legend)
  }

  dev.off()

}



stopCluster(cl)
registerDoSEQ()


# merge png into animation using ffmpeg (or with a different programm)

setwd(tmp_path)
system("ffmpeg -framerate 6 -pattern_type glob -i '*.png' -y -c:v libx264 -profile:v high -crf 1 -pix_fmt yuv420p Bird_tracks_animation.mov")