Setting Up the Data
01 Setting Up the Data.RmdSetting 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 begins 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.
Formatting the Data
Before moving into the data, we could first format the data so that we can go through the workflow smoothly using EmigD_format. This allows us to have a consistent column names as we move through the workflow. The arguments in EmigD_format consist of the necessary columns needed.
deer <- EmigD_format(deer, id_ = "ID", dt_ = "date", x_ = "utm_x", y = "utm_y",
year_ = "year", month_ = "month", jdate_ = "julian")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(dt_, "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.
mutate(start = if_else(day(start) == 31, start - days(10), start)) %>%
group_by(id_, start) %>%
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(dt_, "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.
mutate(start = if_else(day(start) == 31, start - days(10), start)) %>%
group_by(id_, start) %>%
filter(month_ == "7") %>%
group_split()Elements in the list for one individual in January:
head(jan[1:3])
#> <list_of<
#> tbl_df<
#> id_ : character
#> dt_ : date
#> x_ : double
#> y_ : double
#> year_ : double
#> month_: double
#> jdate_: double
#> start : date
#> >
#> >[3]>
#> [[1]]
#> # A tibble: 17 × 8
#> id_ dt_ x_ y_ year_ month_ jdate_ start
#> <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <date>
#> 1 A 2013-01-01 394353. 4222160. 2013 1 1 2013-01-01
#> 2 A 2013-01-04 425764. 4269724. 2013 1 4 2013-01-01
#> 3 A 2013-01-07 426686. 4244820. 2013 1 7 2013-01-01
#> 4 A 2013-01-10 385763. 4176823. 2013 1 10 2013-01-01
#> 5 A 2013-01-03 444235. 4136620. 2013 1 3 2013-01-01
#> 6 A 2013-01-06 281915. 4194640. 2013 1 6 2013-01-01
#> 7 A 2013-01-09 316328. 4247506. 2013 1 9 2013-01-01
#> 8 A 2013-01-02 444606. 4220499. 2013 1 2 2013-01-01
#> 9 A 2013-01-05 413878. 4105887. 2013 1 5 2013-01-01
#> 10 A 2013-01-08 436514. 4239651. 2013 1 8 2013-01-01
#> 11 A 2013-01-01 447905. 4140465. 2013 1 1 2013-01-01
#> 12 A 2013-01-04 399696. 4271535. 2013 1 4 2013-01-01
#> 13 A 2013-01-07 270352. 4163344. 2013 1 7 2013-01-01
#> 14 A 2013-01-10 403482. 4166581. 2013 1 10 2013-01-01
#> 15 A 2013-01-03 343761. 4121445. 2013 1 3 2013-01-01
#> 16 A 2013-01-06 371862. 4129096. 2013 1 6 2013-01-01
#> 17 A 2013-01-09 239170. 4216780. 2013 1 9 2013-01-01
#>
#> [[2]]
#> # A tibble: 16 × 8
#> id_ dt_ x_ y_ year_ month_ jdate_ start
#> <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <date>
#> 1 A 2013-01-13 243076. 4183446. 2013 1 13 2013-01-11
#> 2 A 2013-01-16 277679. 4168675. 2013 1 16 2013-01-11
#> 3 A 2013-01-19 306373. 4188378. 2013 1 19 2013-01-11
#> 4 A 2013-01-12 442920. 4195152. 2013 1 12 2013-01-11
#> 5 A 2013-01-15 441155. 4203151. 2013 1 15 2013-01-11
#> 6 A 2013-01-18 369790. 4240183. 2013 1 18 2013-01-11
#> 7 A 2013-01-11 303512. 4209028. 2013 1 11 2013-01-11
#> 8 A 2013-01-14 308519. 4221374. 2013 1 14 2013-01-11
#> 9 A 2013-01-17 374148. 4162971. 2013 1 17 2013-01-11
#> 10 A 2013-01-20 397508. 4097868. 2013 1 20 2013-01-11
#> 11 A 2013-01-13 320893. 4225999. 2013 1 13 2013-01-11
#> 12 A 2013-01-16 241850. 4209139. 2013 1 16 2013-01-11
#> 13 A 2013-01-19 262219. 4270200. 2013 1 19 2013-01-11
#> 14 A 2013-01-12 257607. 4269495. 2013 1 12 2013-01-11
#> 15 A 2013-01-15 371683. 4170130. 2013 1 15 2013-01-11
#> 16 A 2013-01-18 259174. 4101256. 2013 1 18 2013-01-11
#>
#> [[3]]
#> # A tibble: 19 × 8
#> id_ dt_ x_ y_ year_ month_ jdate_ start
#> <chr> <date> <dbl> <dbl> <dbl> <dbl> <dbl> <date>
#> 1 A 2013-01-22 257936. 4152519. 2013 1 22 2013-01-21
#> 2 A 2013-01-25 379787. 4205195. 2013 1 25 2013-01-21
#> 3 A 2013-01-28 352528. 4174215. 2013 1 28 2013-01-21
#> 4 A 2013-01-31 426091. 4202516. 2013 1 31 2013-01-21
#> 5 A 2013-01-21 412662. 4219820. 2013 1 21 2013-01-21
#> 6 A 2013-01-24 374377. 4181390. 2013 1 24 2013-01-21
#> 7 A 2013-01-27 383701. 4112081. 2013 1 27 2013-01-21
#> 8 A 2013-01-30 268919. 4211473. 2013 1 30 2013-01-21
#> 9 A 2013-01-23 371091. 4139747. 2013 1 23 2013-01-21
#> 10 A 2013-01-26 290570. 4188197. 2013 1 26 2013-01-21
#> 11 A 2013-01-29 384754. 4269268. 2013 1 29 2013-01-21
#> 12 A 2013-01-22 346652. 4207052. 2013 1 22 2013-01-21
#> 13 A 2013-01-25 257383. 4203716. 2013 1 25 2013-01-21
#> 14 A 2013-01-28 251854. 4162917. 2013 1 28 2013-01-21
#> 15 A 2013-01-31 334382. 4218896. 2013 1 31 2013-01-21
#> 16 A 2013-01-21 375032. 4117669. 2013 1 21 2013-01-21
#> 17 A 2013-01-24 299099. 4242995. 2013 1 24 2013-01-21
#> 18 A 2013-01-27 313090. 4107096. 2013 1 27 2013-01-21
#> 19 A 2013-01-30 302671. 4257404. 2013 1 30 2013-01-21We then attach names for each list element in jan and july. In this case, we include the id_ and the year_ for each element. This step is important and it plays an important role when using functions from the EmigD package.
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 x_ and 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 = x_, .y = y_, .t = dt_,
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 = x_, .y = y_, .t = dt_,
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 = "EMigD"))
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 = "EMigD"),
pattern=c("RAP","*.tif"), all.files= FALSE,
full.names = TRUE)
basename(list_names)
#> [1] "2013_RAP.tif" "2014_RAP.tif"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 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 on 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)
#> $year
#> [1] "2013" "2014"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 year 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")