Purpose

Provide some examples of how to visualize the results of track annotations from EnvDATA using base R some popular packages. While there are so many packages and options we don’t attempt to cover everything, this will help get you started, especially newer R users!

Also consider exploring your annotated tracks in DynamoVis (https://doi.org/10.13020/D6PH49).

Preamble

Notes about the libraries:

Load libraries

library(ggmap)
library(ggplot2)
library(leaflet)
library(move)
library(rgdal)
library(rgl)
options(max.print=100)
options(width=95)

Prepare the tracking data.

Load the Svalbard geese dataset annotated with NDVI from Movebank study 31888520, “Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014)”.

These tracking data are published as Griffin (2014) https://doi.org/10.5441/001/1.5k6b1364 and described in Költzsch et al. (2014) https://doi.org/10.1111/1365-2656.12281 and Shariatinajafabadi et al. (2014) https://doi.org/10.1371/journal.pone.0108331

geese.move <- move("annotations/Svalbard geese 1k 16d NDVI-3281758398705165226/Svalbard geese 1k 16d NDVI-3281758398705165226.csv")

See what’s there.

geese.move
## class       : MoveStack 
## features    : 24488 
## extent      : -4.029, 22.421, 54.15566, 78.71467  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +ellps=WGS84 +datum=WGS84 +towgs84=0,0,0 
## variables   : 6
## names       :  event.id,           timestamp,  comments, height.raw, MODIS.Land.Vegetation.Indices.1km.16d.Aqua.NDVI, MODIS.Land.Vegetation.Indices.1km.16d.Terra.NDVI 
## min values  : 369331563, 2006-04-05 07:00:00, 170563-07,       -0.1,                                    1.000006e-02,                                     1.000000e-02 
## max values  : 369356755, 2011-06-18 17:00:00,  86828-09,       NULL,                                   -9.997746e-03,                                    -9.998666e-02 
## timestamps  : 2006-04-05 07:00:00 ... 2011-06-18 17:00:00 Time difference of 1900 days  (start ... end, duration) 
## sensors     : gps 
## indiv. data : visible, algorithm.marked.outlier, sensor.type, individual.taxon.canonical.name, tag.local.identifier, individual.local.identifier, study.name 
## min ID Data : true, NA, gps, Branta leucopsis,  33102, X170563, Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## max ID Data : true, NA, gps, Branta leucopsis, 186827,  X86828, Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## individuals : X33102, X33103, X33104, X33145, X33953, X33954, X64685, X64687, X70564, X70565, X70566, X70567, X70568, X70618, X70619 
## study name  : Migration timing in barnacle geese (Svalbard) (data from Kölzsch et al. and Shariatinajafabadi et al. 2014) 
## date created: 2018-09-12 11:39:55
head(geese.move)
##    event.id           timestamp comments height.raw
## 1 369336493 2011-02-24 15:00:00 33102-11          9
## 2 369336494 2011-02-24 16:00:00 33102-11     2D fix
## 3 369336495 2011-02-24 18:00:00 33102-11          9
## 4 369336514 2011-02-25 06:00:00 33102-11         10
## 5 369336512 2011-02-25 12:00:00 33102-11         12
## 6 369336513 2011-02-25 18:00:00 33102-11         12
##   MODIS.Land.Vegetation.Indices.1km.16d.Aqua.NDVI
## 1                                       0.5143402
## 2                                       0.5151412
## 3                                       0.5150480
## 4                                       0.5150984
## 5                                       0.5159023
## 6                                       0.5160872
##   MODIS.Land.Vegetation.Indices.1km.16d.Terra.NDVI
## 1                                        0.4171148
## 2                                        0.4220628
## 3                                        0.4205474
## 4                                        0.4120822
## 5                                        0.4134778
## 6                                        0.4103168
slotNames(geese.move)
##  [1] "trackId"                 "timestamps"              "idData"                 
##  [4] "sensor"                  "data"                    "coords.nrs"             
##  [7] "coords"                  "bbox"                    "proj4string"            
## [10] "trackIdUnUsedRecords"    "timestampsUnUsedRecords" "sensorUnUsedRecords"    
## [13] "dataUnUsedRecords"       "dateCreation"            "study"                  
## [16] "citation"                "license"
plot(geese.move, type = "b")

Change the long variable name for NDVI provided by EnvDATA. I requested the two 1-km, 16-d MODIS Land NDVI estimates, one from the Aqua satellite and one from the Terra satellite. We’ll use Aqua (see https://doi.org/10.5067/MODIS/MYD13A2.006).

names(geese.move)[names(geese.move)=="MODIS.Land.Vegetation.Indices.1km.16d.Aqua.NDVI"] <- "NDVI"

Check the distribution of NDVI values. Remember there is no vegetation over the ocean. We can use this info to define break points for coloring by value ranges in some of the examples below.

summary(geese.move$NDVI)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##  -0.153   0.067   0.412   0.327   0.529   0.853    4607
hist(geese.move$NDVI)

Create a version of the moveStack as a data frame.

geese <- as.data.frame(geese.move)

Shorten the variable name for animal identifiers.

names(geese)[names(geese)=="individual.local.identifier"] <- "ID"

Map the annotated tracks using base R.

These are scatterplots with no projections or background but a quick way to review results. First, define a color palette. See a nice summary of color palettes at https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/colorPaletteCheatsheet.pdf

cols = colorRampPalette(c("wheat", "forestgreen"))

Rank records by NDVI.

geese$order <- findInterval(geese$NDVI, sort(geese$NDVI))

Plot.

plot(geese$location.long, geese$location.lat, pch = 20, 
     col = cols(nrow(geese))[geese$order],
     xlab = "Longitude (degrees)", ylab = "Latitude (degrees)")

This time let’s instead bin the NDVI values and color by bin, and also add a legend. To do this, define the break points, labels, and then assign the NDVI values to bins.

breaks <- c(-0.2, 0, 0.2, 0.4, 0.6, 0.8, 1.0)
labels <- c("<0", "0-0.2", "0.2-0.4", "0.4-0.6", "0.6-0.8", "0.8-1.0")
bins <- cut(geese$NDVI, breaks, include.lowest = TRUE, labels = labels)

Next add a column assigning a hexadecimal color per bin using our chosen color palette.

geese$color <- cols(7)[as.numeric(bins)]

Plot and add a legend.

plot(geese$location.long, geese$location.lat, col = geese$color,
     xlab = "Longitude (degrees)", ylab = "Latitude (degrees)")
legend("topleft", title = "NDVI", legend = c(labels), col = cols(7), pch = 1, cex = 0.9, 
       box.lty = 0)

Map the annotated tracks using ggplot2 and ggmap.

These packages allow you to work with projections, rasters and web-based background maps, and provide lots of options for customizing graphics. Let’s use the bounding box from our gridded requests to define the map area. Define a lower zoom value (i.e. more zoomed out) to require fewer map tiles to download. Use a black and white background so we can see our NDVI values more clearly.

basemap <- get_stamenmap(bbox = c(left = -6, bottom = 50, right = 24, top = 80), 
               zoom = 4, maptype = "terrain", color = c("bw"))

You could also use the bounding box from your moveStack to define the map extent like this:
basemap <- get_stamenmap(bbox(geese), zoom = 4, maptype = “terrain”, color = c(“bw”))
Plot the locations colored by NDVI.

ggmap(basemap) + geom_point(data = geese, 
                            aes(x = location.long, y = location.lat, color = NDVI))

Customize the coloring and labels.

ggmap(basemap) + 
  geom_point(data = geese, aes(x = location.long, y = location.lat, 
                               color = NDVI), size = 1, alpha = 0.3) +
  scale_color_gradient(low = "wheat", high = "forestgreen") + 
  labs(title = "Svalbard barnacle geese", x = "Longitude", y = "Latitude",
       caption = "Tracking data from Griffin (2014).")

Plot each individual separately.

ggmap(basemap) + 
  geom_point(data = geese, aes(x = location.long, y = location.lat, 
                                  color = NDVI), size = 1, alpha = 0.3) +
  scale_color_gradient(low = "wheat", high = "forestgreen") + 
  labs(x = "Longitude", y = "Latitude") +
  facet_wrap(~ID, nrow = 3)

Map the annotated tracks using leaflet.

Leaflet is a JavaScript library for interactive maps and this package lets you create them. Here are our locations.

leaflet(geese) %>% addTiles() %>%
  addCircles(~location.long, ~location.lat)

As you can see, this package lets you zoom in and out. Helpful! Let’s color by NDVI. Leaflet has pretty specific requirements for color palettes. See https://rdrr.io/cran/leaflet/man/colorNumeric.html and https://rstudio.github.io/leaflet/colors.html.
First define a leaflet color palette for a continuous variable.

cols.leaflet <- colorNumeric(palette = c("wheat", "forestgreen"), domain = geese$NDVI, 
                             na.color = "#808080")

Plot.

map <- leaflet(geese) %>% addTiles() %>%
  addCircles(~location.long, ~location.lat, color = ~cols.leaflet(NDVI))
map

Add a legend.

map %>%
  addLegend("bottomright", pal = cols.leaflet, values = ~NDVI, na.label = "n/a",
            title = "NDVI")

Looks good, but let’s order the legend values from high to low. First we need to make a color palette reversed from the one above.

cols.leaflet.rev <- colorNumeric(palette = c("forestgreen", "wheat"), domain = geese$NDVI, 
                             na.color = "#808080")

Make a function to sort the NDVI values in decreasing order.

values.rev <- function(NDVI) (sort(NDVI, decreasing = TRUE))

Plot.

map %>%
  addLegend("bottomright", pal = cols.leaflet.rev, values = ~NDVI, na.label = "n/a",
            labFormat = labelFormat(transform = values.rev), title = "NDVI")

Now let’s color using the binned NDVI values. We need to define another color palette for a continuous variable with binned values.

cols.leaflet.bin <- colorBin(c("wheat", "forestgreen"), ~NDVI, bins = breaks, 
                             na.color = "#808080", pretty = FALSE)

Plot.

map <- leaflet(geese) %>% addTiles() %>%
  addCircles(~location.long, ~location.lat, color = ~cols.leaflet.bin(NDVI))
map

Add a legend.

map %>%
  addLegend("bottomright", pal = cols.leaflet.bin, values = ~NDVI, na.label = "n/a",
            title = "NDVI")

To order the legend values from high to low, we need to make yet another reversed color palette (geez).

cols.leaflet.bin.rev <- colorBin(c("forestgreen", "wheat"), ~NDVI, bins = breaks, 
                             na.color = "#808080", pretty = FALSE)

Plot.

map %>%
  addLegend("bottomright", pal = cols.leaflet.bin.rev, values = ~NDVI, na.label = "n/a",
            labFormat = labelFormat(transform = values.rev))

3D space-time cube

The rgl package allows you to make interactive 3D graphics. We can use this to explore our results by location and time. First, create a timestamp field that R understands.

geese$time <-as.POSIXct(geese$timestamp, format="%Y-%m-%d %H:%M:%OS", tz="UTC")

Order the data by timestamp.

geese <- geese[order(geese$time),]

Make a space-time cube, colored using our palette from the first example above.

with(geese, plot3d(location.long, location.lat, time, type = "p", 
                   col = cols(nrow(geese))[geese$order]))
rglwidget()

You can click and scroll within the cube to rotate, shrink and enlarge it. Because our color assignments have focused on continuous variables, and because this plot is less interesting because there is no NDVI over the ocean, let’s do one last example and color by individual. Since we have 22 geese, I’ll pull 22 colors from R’s “rainbow” palette.

cols.id <- rainbow(22)

Make a space-time cube colored by animal ID.

with(geese, plot3d(location.long, location.lat, time, xlab = "Longitude", ylab = "Latitude",
                   zlab = "Time", type = "p", col = cols.id[c(ID)]))
rglwidget()

In the next lesson (EnvDATAvizRaster.R) we’ll look at how to work with our rasters and deal with projections.