In jhnhng/trekR: Migration Classifier

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

library(trekR)
library(lubridate)
library(dplyr)
library(tibble)

Setting Up the Data

The example shown below is one of many ways the data could be processed. The important components required is a list object containing home range kernels for each interval that is created for at least one individual.

For example, for a single individual if I wanted to break a one month period into 10-day intervals I would get a total of three 10-day intervals. These intervals should be in a list object, and each element within the list should represent a single interval.

The selected periods to be divided into intervals should be chosen so that they fall between migratory trips and don't overlap them. For example if individuals typically begin their migratory trip in spring and fall, the selected periods should not overlap spring and fall, and should be during the winter and summer seasons.

Create Intervals

We created 10-day intervals for deer in the January and July. We first created a start column using the floor_Date() function to obtain the 10-day extents within the months. We included the 31st day (if a month had 31 days) in the third 10-day interval to prevent the creation of a 4th 10-day interval with only a single day in it.

Intervals for January

jan <- deer %>%
  # Creates a start column assigning the first day in the 10-day interval in which
  # the date falls under (e.g., 01-03-2021 would be in the first 10-day interval
  # so the `floor_date` assigned to it would be 01-01-2021)
  mutate(start = floor_date(date, "10 days")) %>%
  # For any months that has 31 days, the 31st day would normally be assigned its 
  # own interval. The code below takes the 31st day and joins it with the 
  # previous interval. 
  group_by(ID) %>% 
  mutate(start = if_else(day(start) == 31, start - days(10), start)) %>% 
  group_by(start, .add = TRUE) %>%
  filter(month == "1") %>% 
  group_split()

Intervals for July

july <- deer %>%
  # Creates a start column assigning the first day in the 10-day interval in which
  # the date falls under (e.g., 01-03-2021 would be in the first 10-day interval
  # so the `floor_date` assigned to it would be 01-01-2021)
  mutate(start = floor_date(date, "10 days")) %>%
  # For any months that has 31 days, the 31st day would normally be assigned its 
  # own interval. The code below takes the 31st day and joins it with the 
  # previous interval. 
  group_by(ID) %>% 
  mutate(start = if_else(day(start) == 31, start - days(10), start)) %>% 
  group_by(start, .add = TRUE) %>%
  filter(month == "7") %>% 
  group_split()

Elements in the list for one individual in January:

head(jan[1:3])

Create Tracks for Each Interval

Before we estimate LoCoH home range for each of these 10-day intervals, we first create objects of class track_xyt using make_track() from the amt package. This function creates a track using the utm_x and utm_y columns from jan and july. lapply() allows us to apply the specified function across a list object (e.g., jan and july)

track_list_sum <- lapply(july, function(x) {
  amt::make_track(tbl = x, .x = utm_x, .y = utm_y, .t = date,
             uid = ID,
             # lat/long: 4326 (lat/long, WGS84 datum).
             # utm: crs = sp::CRS("+init=epsg:32612"))
             crs = 32612)
})
track_list_win <- lapply(jan, function(x) {
  amt::make_track(tbl = x, .x = utm_x, .y = utm_y, .t = date,
             uid = ID,
             crs = 32612)
})

Creates LoCoHs Using the Tracks

Using the track created for the jan and july, we estimate type "a" LoCoH home range using hr_locoh() from the amt package. If you would like to learn more about the arguments of this function please refer to ?amt::hr_locoh for more details.

sum_locoh <- lapply(track_list_sum, function(tracks){
    # Calculate LoCoH a*
    dmat <- dist(tracks[, c("x_", "y_")])
    a <- max(dmat)

    # Fit LoCoH
    locoh <- try(amt::hr_locoh(
      x = tracks,
      levels = seq(0.1, 1, by = 0.1),
      keep.data = TRUE,
      n = a,
      type = "a",
      rand_buffer = 1e-05
    ))


})
win_locoh <- lapply(track_list_win, function(tracks){
    # Calculate LoCoH a*
    dmat <- dist(tracks[, c("x_", "y_")])
    a <- max(dmat)

    # Fit LoCoH
    locoh <- try(amt::hr_locoh(
      x = tracks,
      levels = seq(0.1, 1, by = 0.1),
      keep.data = TRUE,
      n = a,
      type = "a",
      rand_buffer = 1e-05
    ))


})

Creates Isopleths from the LoCoHs

We then created isopleths from the LoCoH home range estimates using hr_isopleths()from the amt package.

sum_iso <- lapply(sum_locoh, function(x){
  amt::hr_isopleths(x)
})
win_iso <- lapply(win_locoh, function(x){
  amt::hr_isopleths(x)
})

Rasterizes the Isopleths

After creating isopleths for each of the intervals, we have to rasterize them. Similar to the the previous functions, we use lapply() to apply our function across our list object. We first select the resolution for that we want for our rasters (res). We then create a template raster with a resolution of res (keep in mind that the speed the isopleths are rasterize() will be dependent on the specified resolution, with time of processing increasing with smaller resolutions). The projection of the template raster will match that of the specified isopleth.

After creating rasters for each interval, we then subtract all the isopleth values from 1, then normalize them.

res <- 250

sum_raster <- lapply(sum_iso, function(iso){
# Creates template raster
    template_raster <- raster::raster(iso,
      resolution = res, vals = 0,
      crs = sp::CRS(sf::st_crs(iso)[[2]])
    )
    # Rasterizes the isopleth
    raster_locoh <- raster::rasterize(iso,
      template_raster,
      field = "level",
      fun = "first"
    )

    # Now you want to subtract all the isopleth values from 1
    raster::values(raster_locoh) <- 1 - raster::values(raster_locoh)

    # Now normalize
    x <- raster_locoh / raster::cellStats(raster_locoh, sum)
})

win_raster <- lapply(win_iso, function(iso){
# Creates template raster
    template_raster <- raster::raster(iso,   
      resolution = res, vals = 0,
      crs = sp::CRS(sf::st_crs(iso)[[2]])
    )
    # Rasterizes the isopleth
    raster_locoh <- raster::rasterize(iso,
      template_raster,
      field = "level",
      fun = "first"
    )

    # Now you want to subtract all the isopleth values from 1
    raster::values(raster_locoh) <- 1 - raster::values(raster_locoh)

    # Now normalize
    x <- raster_locoh / raster::cellStats(raster_locoh, sum)
})

names(win_raster) <- sapply(jan, function(x) paste(x$ID[1],
                                            x$year[1], sep = '_'))

names(sum_raster) <- sapply(july, function(x) paste(x$ID[1],
                                              x$year[1], sep = '_'))

Attach attributes to each list of rasters

We attach attributes to our list of rasters to allow for more flexibility in working with our lists of rasters.

# Add an attribute to hold the attributes of each list element
attributes(win_raster) <- data.frame(id = sapply(jan, function(x) paste(x$ID[1])),
                                      interval_start_date = sapply(jan, function(x) paste(x$start[1])),
                                      year = paste0(sapply(jan, function(x) paste(x$year[1])))
)

# Check the attributes
attributes(win_raster)

# Add an attribute "tab" to hold the attributes of each list element
attributes(sum_raster) <- data.frame(id = sapply(july, function(x) paste(x$ID[1])),
                                      interval_start_date = sapply(july, function(x) paste(x$start[1])),
                                      year = paste0(sapply(july, function(x) paste(x$year[1])))
)

# Check the attributes
attributes(sum_raster)

Setting up input for Calculating EMDs in Environmental Space

The data set up above have specifically focused on setting up the data for calculating EMDs in geographical space. The code below demonstrates how you can set up the data to calculate the EMDs in environmental space.

Creating 'RasterStacks' to calculate EMD and 'EM-Speeds in Environmental Space

Environmental Covariates without Temporal Variation

We first import the environmental covariate file (e.g., DEM layer), then we re-project the extent from latitude and longitude to UTMs to match those of the home range rasters.

  #Load raster
  UtahDEM <- raster::raster(system.file("extdata", "Utah_DEM.tif", package = "trekR"))
  out_proj <- "+proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"
  UtahDEM@extent <- raster::extent(raster::projectExtent(UtahDEM, out_proj))

The env_raster_stack() function takes the environmental covariate raster and resamples it to the home range rasters. It then takes the resampled rasters and the home range rasters and constructs a object with class RasterStack. The steps following the creation of these RasterStacks is the steps under "Earth Mover's Distance", where the RasterLayers objects (i.e., sum_raster and win_raster) are replaced with the RasterStack objects (i.e., sum_rst and win_rst).

sum_rst <- env_raster_stack(UtahDEM, sum_raster)
win_rst <- env_raster_stack(UtahDEM, win_raster)

Environmental Covariates with Temporal Variation

For environmental covariates that vary across time (e.g., temperature, snow depth, foraging resources, precipitation, etc.), we first assign the file names of which the covariates are housed to list_names. In our example, we use Rangeland Analysis Platform (RAP) data in 2013 and 2014 that can be an indicator of forage resources on the landscape.

list_names <- list.files(path=system.file("extdata", package = "trekR"),
                 pattern=c("RAP","*.tif"), all.files= FALSE,
                 full.names = TRUE) 

basename(list_names)

We then read in those files into R as RasterLayer objects. For each of those list elements in env_lists we assign the path names for each of those files. Using the assigned names, we created attributes(metadata) for our list with a year column. The year column may change depending on the temporal scale of your environmental covariate. For example, if we had environmental covariates that changed daily, we can first aggregate the daily covariates based the intervals. Instead of using year, we would then use the start date of the intervals.

# Rasterize the environmental covariates
env_list <- lapply(list_names, raster::raster)

names(env_list) <- sapply(list_names, function(x) paste(x))

# Uses the name of the environemntal covariate file to assign attributes 
attributes(env_list) <- data.frame(year = trimws(basename(names(env_list)), 
                                                 whitespace = "_.*"))

attributes(env_list)

We use env_raster_stack with the addition of the attribute argument to stack environmental covariate layers with our home range rasters. The input for the attribute argument should be a attribute component that is the same in both of the list objects (i.e., env_list and sum_raster/win_raster) with the same unique values. For example, both env_list and sum_raster/win_raster have a attribute component labelled year and both have the same unique values (e.g., 2013 and 2014).

NOTE: Currently this feature assumes that only environmental layers relevant to sum_raster or win_raster are being used and are in the same folder(i.e., excluding any that are not present in sum_raster and/or win_raster.)

sum_rst_temporal <- env_raster_stack(env_list, sum_raster, attribute = "year")
win_rst_temporal <- env_raster_stack(env_list, win_raster, attribute = "year")

jhnhng/trekR documentation built on April 25, 2022, 12:43 p.m.