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

knitr::opts_chunk$set(echo = TRUE)

WindR

windR is a R package to connect animal tracking data with wind data (or sea current data) and allows to visualize animal movements within the flow of the medium that they moved in.

The package was written to compile functions used to analyse pectoral sandpiper movements in the wind. A full animation of these tracks can be found on youtube: Moving though the Arctic: Pectoral sandpipers in the wind

What can you do with windR?

windR uses wind data from ERA-Interim (a global atmospheric reanalysis model) described in detail in Dee et al. 2011 and connects them wind tracking data (your own). It can be used to create particle flow animations of wind or sea current data. To connect them with tracking data, it is necessary to calculate the bearing (ground direction), ground speed, wind support and cross winds from animal tracks using an equal area map projection (e.g.Lambert azimuthal equal-area projection). The wind support represents the length of the wind vector in the direction of the birds’ flight and the cross wind represents the length of the wind vector perpendicular to ground vector (see Safi et al. 2013 for schematic representation). For a detailed description see the workflow example described below.

# Packages
sapply(c('magrittr', 'data.table', 'windR', 'raster', 'ggplot2', 'stringr', 'grid', '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 particles and calculate the wind speed
particle_data = system.file('Map', 'Wind_particles.RDS', package = 'windR')
ob = readRDS(particle_data)
setkey(ob, datetime_)
ob[, w_speed := sqrt(u^2 + v^2)]

# Load tracks with wind support
d = PESA_example_ws
d[, datetime_ := as.POSIXct(datetime_)]

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

# 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()

### 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))


# 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')

# 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)

# 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_)

# Example plot
i = 45

# add track for particles
tail     = 20           # lenght of the running tail
tmp_date = ts[i]$date   # current date

# subset particles
tmp_date_sub = seq(tmp_date  - 1800 * tail, tmp_date, by = '30 mins')
pi = ob[J(tmp_date_sub), nomatch=0L]

# add wind particles to base map
if (nrow(pi) > 0) pi[, a :=  alphaAlong(datetime_, head = 70, skew = -2), by = particle_id] # alpha

bm_w = bm +
  geom_path(data = pi, aes(x = x, y = y, group = particle_id, color = w_speed), alpha = pi$a) +
  scale_colour_gradientn(colours = col_Ws, values = col_WS_v,
                         limits = c(min(d$Ws, na.rm = T), max(c(d$Ws, ob$w_speed), na.rm = T)), name = 'Wind support (m/s)') +
  annotate('text', x = 0, y = -1740000,
           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'))

# add 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.2, 3)) , 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)


suppressWarnings(print(gi)); grid.draw(legend)

The figure shows an example snapshot of the particle flow animation including two male Pectoral Sandpipers (thick comets in light green) that left the area around Barrow (northern tip of Alaska) coloured with the wind support (m/s) and wind particles coloured with the wind speed (m/s; same scale as wind support) flying around based on the wind data at this time. Map projection: polar Lambert azimuthal equal-area with longitude origin 156.65° W (Barrow); Map data from Natural Earth

Workflow example

The vignettes give a small step by step example of what one has to do to reach the final result (a particle flow animation with animal tracks). The single vignettes are built up on each other, but each can be run independently (since the output data of each step can be loaded from the package data). The vignettes as html can be found at OSF and viewed in a browser after downloading.

The first vignette A_ERA_Interim_data_download describes how to download the ERA-Interim data using a python script. Note that single month can also be downloaded via the website directly.

The second vignette B_Wind_data_manipulation describes how to open the wind data, interpolate them to a higher resolution and transform them in a data.table including the date and u- & v-wind component.

The third vignette C_Wind_particle_flow describes on how to use wind data to calculate a particle flow (how to create particles) and how to create an animation with them.

The forth vignette D_Wind_support_and_track_animation describes how to connect animal tracks (using a subset of data from Kempenaers & Valcu 2017) with wind data and how to calculate the bearing, ground speed, wind support and cross winds from the tracks. Afterwards, it gives an example of how to plot the tracking data in a simple ggplot and how to do a comet plot animation using tracking data.

The fifth vignette F_Wind_animation_with_tracks brings everything together. It combines the particle flow animation of the wind data with the tracking data.

Tips to get started with your own data

• Follow the vignettes and run the examples.
• Start small with your own data! Take a small subset of your data to go thru the whole process (of whatever you want to do). Starting with a big data set makes everything slow and will make you soon frustrated. If everything is as you which with your small subset – go big and create particles for the whole period of your interest.
• Find the optimal spatial and temporal scale fitting to your tracking data and question.
• Have fun!

Installation

install.packages('devtools')
devtools::install_github('mpio-be/windR')

# install with vignettes
devtools::install_github("mpio-be/windR", build_vignettes = TRUE, force = TRUE)
vignette(package = "windR")

Limitations

All analyses are constrained by the spatio-temporal resolution of the used wind and tracking data.

Inspiration & Background

This project was inspired by the awesome earth project from Cameron Beccario, which was itself inspired by the wind map project from HINT.FM.

We saw these particle flow maps of the wind and wanted to see our bird tracks within such a visualization. Both earth and wind map use fixed times of global wind data sets (one wind layer). Particles are randomly thrown in the map and move based on the wind speed and direction. To plot our bird tracks within the wind, we needed to find a way to continuously change the wind data with the time that the birds moved. We did so by always using the closest wind data in time (continuously changing the wind layers), resulting in a dynamic flow of the wind particles changing over time. We did our best to make this workflow fast in R, but know that using other programming languages (i.e. JavaScript) could improve the speed of these analyses. We are happy if somebody wants to improve (speed-up) this script or translates parts (esp. the particle creation) into another programming language.

mpio-be/windR documentation built on Feb. 2, 2020, 10:04 a.m.