Animating snow leopard telemetry data using the moveVis package to motivate transboundary conservation.


Introduction

Oblivious to the often invisible administrative boundaries delineated by humans, animals move freely through their environments, but such anthropogenic constructs can seriously confine conservation efforts. In recent years, though, there has been a new push for conservation efforts across political boundaries, with the goal of developing initiatives in a contiguous fashion, insuring that protective measures are thorough and robust: transboundary conservation. These efforts have been focused in particular on wide-ranging carnivores, known to inhabit remote areas around the world and decidedly unconstraied from any ecologically-arbitrary borders we might impose.

Explicitly observing an individual crossing a border between conuntries might be nearly impossible, but with the increasing colleciton of animal movement data, it can now be (almost) a reality. My colleague, Ali Nawaz, sent me a dataset he has from a single gps-collared snow leopard, traversing the mountains of the Hindukuush range in northern Pakistan, just west of the town of Chitral.

Below is the process for turning this data set into a highly-informative animation, clearly showing the seasonal elevational migration of the species as well as its exceptional aptititude for long-distance movements, crossing not only borders, but 5,500-meter mountain peaks.


Setup

# Tidyverse
library(lubridate)
library(dplyr)
library(stringr)

# Tidy plotting
library(ggplot2)
library(ggnewscale)

# Animal movement
library(move)
library(moveVis)

# Spatial data
library(sp)
library(raster)
library(leaflet)

# Resolving namespace conflict
select =  dplyr::select


Explore & format data

# Load all GPS data
gps <- read.csv("data/gps/ali_sl_gps_fix.csv", row.names = NULL) %>%
  mutate(timestamps = ymd_hms(paste(Date, Time))) %>% # Convert date to POSIXct
  select(X = Longitude,
         Y = Latitude, 
         timestamps) %>%
  mutate(track_id = "A") %>%
  select(X, Y, timestamps, track_id) %>%
  filter(X < 100) # Remove bad fixes

# Plot gps
ggplot(data = gps, aes(x=X,y=Y)) +
  geom_point(size=1.75, alpha=0.3, pch = 20) +
  coord_fixed() +
  theme_bw()


Assign CRS

# Assign CRS to convert
gps_spdf <- SpatialPointsDataFrame(
  coords = gps[,c("X", "Y")], 
  data = gps,
  proj4string = CRS("+init=epsg:4326"))

# Plot spatial points
plot(gps_spdf)

# Check to make sure it makes sense
leaflet(gps_spdf) %>% addTiles() %>%
  addMarkers()


Import GIS data

# International borders
borders <- readOGR("data/gis/international_border/WorldAdministrativeDivisions.shp") %>%
  crop(., 4*extent(spTransform(gps_spdf, crs(.)))) %>%
  spTransform(., crs(gps_spdf))

# DEM
dem_raw <- raster("data/gis/dem/Elevation_UTM.tif") %>%
  crop(., 5*extent(spTransform(gps_spdf, crs(.)))) # crop to gps pts

# Reproject dem and crop
dem <- projectRaster(dem_raw, crs = crs(gps_spdf)) %>%
  crop(., 4*extent(gps_spdf))


Improve DEM visuals

# Shading
slope <- terrain(dem, opt="slope", unit='radians')
aspect <- terrain(dem, opt="aspect", unit='radians')
hillshade <- hillShade(slope, aspect)
plot(hillshade, col = gray.colors(20))

# Terrain colors
sl.terrain = colorRampPalette(c("#73866f", "#73866f", "#73866f", "#c0b195", "#c0b195", "#c0b195", "white", "white"))

# Test plot
plot(dem, col = sl.terrain(1000), alpha=1.5, zlim = c(1200, 5400))
plot(borders, lwd=2, add=T)
plot(hillshade, col = gray.colors(1000), alpha = 0.2, add=T, legend = F)
plot(slope, col = gray.colors(1000), alpha = 0.1, add=T, legend = F)
plot(spTransform(gps_spdf, crs(dem)), add=T, col = "blue", pch = 1, cex=0.5)


Prepare data for moveVis

# Convert from SPDF to DF
gps4move <- gps_spdf %>%
  as.data.frame()

# Convert to move data
gps_move <- df2move(
  gps4move, x = "X.1", y = "Y.1", time = "timestamps", track_id = "track_id",
  proj = "+init=epsg:4326")

# Align move_data to a uniform time scale
sl_move_alinged <- align_move(
  gps_move, 
  res = 480, # should be 240
  digit = 0, 
  unit = "mins") # should be mins

# Start time for raster basemap
starttime = min(gps4move$timestamps)


Design the animation

# Make the map frames
frames.sp <- frames_spatial(
  m = sl_move_alinged, 
  r_list = dem, r_times = starttime, r_type = "gradient", fade_raster = F,
  path_colours = c("red"), path_alpha = 1, path_legend = FALSE, alpha = 1) %>%
  add_gg(frames = ., gg = expr(scale_fill_gradientn(colours = sl.terrain(1000)))) %>%
  add_gg(frames = ., gg = expr(labs(fill="Elevation (m)"))) %>%
  add_gg(frames = ., gg = expr(new_scale("fill"))) %>%
  add_timestamps(type = "label") %>% 
  add_progress() %>%
  add_gg(frames = ., gg = expr(geom_raster(data = hills.df, aes(lon, lat, fill = hills, group=1), alpha = .2, interpolate=TRUE)), data = hills.df) %>%
  add_gg(frames = ., gg = expr(geom_raster(data = slope.df, aes(lon, lat, fill = slope, group=2), alpha = .1, interpolate=TRUE)), data = slope.df) %>%
  add_gg(frames = ., gg = expr(scale_fill_gradientn(colours = grey.colors(100), guide='none'))) %>%
  add_gg(frames = ., gg = expr(geom_polygon(data = borders.fort, aes(x = long, y = lat, group = group), 
                                              colour = "#262322", fill = NA, size = 1)), data = borders.fort) %>%
  #add_gg(frames = ., gg = expr(geom_text(aes(label = "Pakistan"), x = 71.6, y = 35.95, color = "white", fontface = "bold"))) %>%
  #add_gg(frames = ., gg = expr(geom_text(aes(label = "Afghanistan"), x = 71.25, y = 35.7, color = "white", fontface = "bold"))) %>%
  add_labels(x = "Longitude", y = "Latitude") %>% 
  add_text("Pakistan", x = 71.6, y = 35.95, colour = "#262322", size = 6) %>%
  add_text("Afghanistan", x = 71.25, y = 35.775, colour = "#262322", size = 6) %>%
  add_scalebar(height = 0.01)


Render animation

# Make final animation
animate_frames(frames.sp, out_file = "output/sl_moveVis.mov")