03a_MoveTraits_database_Tucker.Rmd

---
title: "MoveTraits Database"
author: "Anne G. Hertel"
date: "7/4/2024"
output: 
  html_document:
    theme: cosmo
    toc: yes
    toc_float:
      collapsed: true
    code_folding: show
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r packages, warning=FALSE,message=FALSE}
library(lubridate);library(metafor);library(tidyverse);library(amt);library(adehabitatLT);library(adehabitatHR)
```

# Aim

We here provide code for a first version of the MoveTraits database. MoveTraits uses animal movement data collected from GPS sensors to summarize a suite of movement metrics on the individual level. 

The workflow is as follow:
1. resampling of the raw GPS relocation data to regular time intervals, here: 1hour, 24 hour
2. using resampled data to quantify movement metrics 
  a) using these resampled data to build regular movement trajectories to quantify step lengths of successive locations (at the hourly or 24 hourly rate)
  b) using resampled data to quantify maximum displacement within a set time interval from all pairwise distance comparisons
  c) 
3. summarize each metric to obtain one value per individual (mean, median, cv, 5 & 95 %ile)
4. create a database with metrics summarized at the indivdual level AND provide the raw metrics without spatial information

We here compile a first version from open access data obtained from Tucker et al. 2023 "Behavioral responses of terrestrial mammals to COVID-19 lockdowns" (https://zenodo.org/records/7704108, file ""). 

We only used data from 2019 (i.e., not from 2020 during COVID lockdowns).

## Load a prepare raw spatial data

```{r load data}
movedata <- readRDS("./DATA/Tucker/Tucker_Road_Spatial.rds")
```

```{r}
colnames(movedata)
```

```{r select columns}
movedata <- movedata[,c("Species","ID","TimestampUTC","Latitude", "Longitude","BodyMass_kg","ContactPerson")]

colnames(movedata) <- c("species","individual_id","TimestampUTC","Latitude","Longitude", "animal_mass","contact_person_name")

head(movedata)
```

Add time information to aggregate later; set timezone to UTC

```{r wrangle data}
movedata <-
  movedata %>% 
  mutate(
    OrdinalDay = lubridate::yday(TimestampUTC),
    year = lubridate::year(TimestampUTC),
    month = lubridate::month(TimestampUTC),
    hour = lubridate::hour(TimestampUTC))

```

Keep only data from 2019, i.e. not during Covid lockdowns
```{r select time}
movedata <- movedata %>% 
  filter(year < 2020) 
```

# Resample data to 1h, 24h time scales

We will aggregate movement traits at different temporal time scales. For this we resample the dataset to hourly position and to one position per day. 

Here is the function to resample the data, based on the amt() package.
```{r function.resample}
make_and_resample_track <- function(x, sampling.interval = 1, tolerance = 15){
  require(amt)
  make_track(x, 
             Longitude, 
             Latitude, 
             TimestampUTC, 
             individual_id = individual_id, 
             species = species, 
             contact_person_name = contact_person_name, 
             animal_mass = animal_mass,
             OrdinalDay = OrdinalDay, 
             year = year, 
             month = month, 
             hour = hour, 
             crs = 4326) %>%
    track_resample(., rate = hours(sampling.interval),
                   tolerance = minutes(tolerance)) %>%
    return()
}
```

Resample to 1hr and to 24 hrs. We start with unstratified

```{r resample locs, warning=FALSE,message=FALSE}
animlocs.1hourly <- movedata %>% 
  group_by(individual_id) %>% 
  group_split() %>% 
  map(~ make_and_resample_track(., sampling.interval = 1, 
                                 tolerance = 15)) %>% 
  setNames(unique(movedata$individual_id)) 

animlocs.daily <- movedata %>%
  group_by(individual_id) %>%
  group_split() %>%
  map(~ make_and_resample_track(., sampling.interval = 24,
                                tolerance = 15)) %>%
  setNames(unique(animlocs.1hourly$individual_id)) 

```

# Calculate movement metrics

## Displacement distances
The distance moved over a given time period, based on a specific GPS relocation rate.

### 1hr displacement distance
As input to calculate hourly displacement we use the resampled hourly dataset. We only include individuals with at least 167 hourly relocations, which is equivalent to 1 week of data.

```{r prepare 1h locs}
animlocs.1hourly_sl <- data.frame()

for (i in 1:length(animlocs.1hourly)) {
  temp_df <- animlocs.1hourly[[i]] %>% 
    filter(n() > 167) %>% 
    mutate(d1h = step_lengths(.)*100000) |> 
    mutate(time_diff = as.numeric(difftime(lead(t_),t_, units = "mins"))) |> 
    filter(time_diff >= 45 & time_diff <= 75) |> 
    dplyr::select(individual_id,t_, d1h, x_, y_) |> 
    mutate(ymd = as.character(format(as.Date(t_), "%Y-%m-%d"))) |> 
    filter(!is.na(d1h))

animlocs.1hourly_sl <- rbind(animlocs.1hourly_sl,as.data.frame(temp_df))
}
```

### Maximum 24hr displacement
Group all GPS relocations within a 24hr interval. Include only days with at least 12 relocations. Calculate all pairwise distances and select the maximum distance observed. Keep only individuals with at least 7 days of data

```{r max 24hr displacement, class.source = 'fold-hide'}
locs1h <- animlocs.1hourly %>% 
  do.call("rbind", .) |>  
  mutate(ymd = as.character(format(as.Date(t_), "%Y-%m-%d"))) |> 
  mutate(id.ymd = paste(individual_id,ymd,sep="_")) |> 
  group_by(id.ymd) %>% 
  filter(n() >= 12) %>%
  ungroup() 

mean.coord <- locs1h |> group_by(id.ymd) |> 
  mutate(mean.x = mean(x_, na.rm=T),
         mean.y = mean(y_, na.rm=T)) |> 
  dplyr::select(id.ymd, mean.x, mean.y) |> distinct()

locs1h.sf <- sf::st_as_sf(locs1h,
                          coords = c("x_", "y_"),
                          crs = 4326)

moveObjSplitTime <- split(locs1h.sf, locs1h$id.ymd)
maxNetDispL <- lapply(moveObjSplitTime, function(x){max(sf::st_distance(x))})
maxNetDisp <- do.call("rbind",maxNetDispL)

rm(moveObjSplitTime);rm(maxNetDispL);rm(locs1h);rm(locs1h.sf)

dmax24 <- 
    data.frame(keyName=row.names(maxNetDisp), dmax24h=maxNetDisp[,1], row.names=NULL) |> 
      mutate(ymd = str_sub(keyName, -10),
             individual_id = str_extract(keyName, "^[^_]*_[^_]*")) |> 
  filter(!is.na(dmax24h)) |>  
  group_by(individual_id) |>  
  filter(n() >= 7) %>%
  ungroup() |> 
  mutate(id.ymd = paste(individual_id,ymd,sep="_")) |> 
  left_join(mean.coord, by = "id.ymd") |> 
  dplyr::select(individual_id,ymd,dmax24h, mean.x, mean.y) 

rm(mean.coord);rm(locs1h)
```

### 24hr displacement distance
Straight line displacement distance between 24hr relocations. Include only individual which have relocation data on at least 30 days (1 month of data)

```{r 24hr displacement, class.source = 'fold-hide'}
animlocs.daily_sl <- data.frame()

for (i in 1:length(animlocs.daily)) {
  #i = 1
  temp_df <- animlocs.daily[[i]] %>% 
    filter(n() > 30) %>% 
    mutate(d24h = step_lengths(.)*100000) |> 
    mutate(time_diff = as.numeric(difftime(lead(t_),t_, units = "hours"))) |> 
    filter(time_diff >= 20 & time_diff <= 28) |> 
    dplyr::select(individual_id,t_, d24h, x_, y_) |> 
    mutate(ymd = as.character(format(as.Date(t_), "%Y-%m-%d")),
         month = lubridate::month(t_),
         year = lubridate::year(t_)) |> 
    filter(!is.na(d24h))

animlocs.daily_sl <- rbind(animlocs.daily_sl,as.data.frame(temp_df))
}
```

### Maximum 7day displacement distance
Based on daily (24h) relocations we calculate the maximum weekly displacemnt from all pairwise comparisons including only weeks with at least 5 relocation and individuals with at least 10 weeks of data.

```{r class.source = 'fold-hide'}
locs24h <- animlocs.daily %>% 
  do.call("rbind", .) |>  
  mutate(week = as.numeric(strftime(t_,format="%W")), year = as.numeric(strftime(t_,format="%Y")),
         year_week = paste(year,week, sep="_")) |> 
  mutate(id.year_week = paste(individual_id,year_week,sep="_")) |>  
  group_by(individual_id,year_week) |> filter(n() >= 5) |> ungroup()

mean.coord <- locs24h |> group_by(id.year_week) |> 
  mutate(mean.x = mean(x_, na.rm=T),
         mean.y = mean(y_, na.rm=T)) |> 
  dplyr::select(id.year_week, mean.x, mean.y) |> distinct()

locs24h.sf <- sf::st_as_sf(locs24h,
                          coords = c("x_", "y_"),
                          crs = 4326)

moveObjSplitTime <- split(locs24h.sf, locs24h.sf$id.year_week)
maxNetDispL <- lapply(moveObjSplitTime, function(x){max(sf::st_distance(x))})
maxNetDisp <- do.call("rbind",maxNetDispL)

rm(moveObjSplitTime);rm(maxNetDispL);rm(locs24h);rm(locs24h.sf)

dmax7d <- 
  data.frame(keyName=row.names(maxNetDisp), dmax7d=maxNetDisp[,1], row.names=NULL) |> 
  mutate(week = str_sub(keyName, -2),
         year_week = str_sub(keyName, -7),
         individual_id = str_extract(keyName, "^[^_]*_[^_]*")) |> 
  filter(!is.na(dmax7d))|> 
  group_by(individual_id) |>  filter(n() >= 10) |> ungroup() |>  
  mutate(id.year_week = paste(individual_id,year_week,sep="_")) |> 
  left_join(mean.coord, by = "id.year_week") |> 
  dplyr::select(individual_id,week, year_week,dmax7d,mean.x, mean.y) 

rm(mean.coord)
```

## Range size - MCPs
The are used within a given time period. We opted to calculate area use as simple 95% Minimum Convex Polygons.

### Daily range size (MCPs)
Daily range used based on hourly relocations including only individuals with at least 12 locations on a given day.

```{r,message=FALSE,warning=FALSE}
dat.mcp.daily <- animlocs.1hourly %>% 
  do.call("rbind", .) |>  
  tibble() %>% mutate(ymd = as.character(format(as.Date(t_), "%Y-%m-%d")))  %>%
  filter(!is.na(x_)) %>% filter(!is.na(y_)) %>% 
  mutate(id.day = paste(individual_id,ymd,sep=".")) %>% group_by(id.day) %>% filter(n() > 12) %>% ungroup() %>% 
  dplyr::select(x_,y_,id.day) 

      mean.coord <- dat.mcp.daily |> group_by(id.day) |> 
        mutate(mean.x = mean(x_, na.rm=T),
               mean.y = mean(y_, na.rm=T)) |> 
        dplyr::select(id.day, mean.x, mean.y) |> distinct()

coordinates(dat.mcp.daily) <- c("x_","y_")
proj4string(dat.mcp.daily) <- CRS("+init=epsg:4326")

# Bonne equal area projection - https://spatialreference.org/ref/esri/54024/
# "+proj=bonne +lat_1=60 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"
# Mollweide projection
dat.mcp.daily <- spTransform(dat.mcp.daily,sp::CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"))

mcp.daily <- mcp(dat.mcp.daily, percent = 95, unout = c( "m2")) %>% data.frame() %>% 
    left_join(mean.coord, by = c("id" = "id.day")) |> 
  mutate(ymd = str_split(id, '[.]', simplify = TRUE)[,2],
         individual_id = str_split(id, '[.]', simplify = TRUE)[,1]) %>% 
  dplyr::select(individual_id,ymd,area, mean.x, mean.y)
rownames(mcp.daily) <- NULL  

rm(mean.coord);rm(dat.mcp.daily)
```

### Weekly MCPs
Weekly range size based on 1h data with at least 84 hourly relocations (3.5 days)

```{r class.source = 'fold-hide',message=FALSE,warning=FALSE}
dat.mcp.weekly <- animlocs.1hourly %>% 
  do.call("rbind", .) |>  
 mutate(week = as.numeric(strftime(t_,format="%W")), 
        year = as.numeric(strftime(t_,format="%Y")),
         year_week = paste(year,week, sep="_")) %>% 
  mutate(id.week = paste(individual_id,year_week,sep=".")) %>% 
  filter(!is.na(x_)) %>% filter(!is.na(y_)) %>%
  group_by(id.week) %>% filter(n() > 84) %>% ungroup() %>% 
  dplyr::select(x_,y_,id.week) 

      mean.coord <- dat.mcp.weekly |> group_by(id.week) |> 
        mutate(mean.x = mean(x_, na.rm=T),
               mean.y = mean(y_, na.rm=T)) |> 
        dplyr::select(id.week, mean.x, mean.y) |> distinct()

coordinates(dat.mcp.weekly) <- c("x_","y_")
proj4string(dat.mcp.weekly) <- CRS("+init=epsg:4326")

dat.mcp.weekly <- spTransform(dat.mcp.weekly,
                              sp::CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"))

mcp.weekly <- mcp(dat.mcp.weekly, percent = 95, unout = c("m2")) %>% data.frame() %>%
  left_join(mean.coord, by = c("id" = "id.week")) |> 
  mutate(week = as.numeric(stringr::str_extract(id, "(\\d+$)")),
         year_week = stringr::str_extract(id, "[^.]*$"),
         individual_id = str_extract(id, "[^.]+")) |> 
  dplyr::select(individual_id, week, year_week, area, mean.x, mean.y)

rm(dat.mcp.weekly);rm(mean.coord)
```

### Monthly MCPs
Monthly MCPs with based on daily relocations, include only individulas with at least 14 days of data

```{r class.source = 'fold-hide',message=FALSE,warning=FALSE}
dat.mcp.monthly <- animlocs.daily %>% 
  do.call("rbind", .) |>   
  mutate(year_month = paste(year,month,sep="_")) |> 
  mutate(id.month = paste(individual_id,year_month,sep=".")) %>% 
  filter(!is.na(x_)) %>% filter(!is.na(y_)) %>% group_by(id.month) %>% 
  filter(n() > 14) %>% ungroup() %>% dplyr::select(x_,y_,id.month) 

      mean.coord <- dat.mcp.monthly |> group_by(id.month) |> 
        mutate(mean.x = mean(x_, na.rm=T),
               mean.y = mean(y_, na.rm=T)) |> 
        dplyr::select(id.month, mean.x, mean.y) |> distinct()

coordinates(dat.mcp.monthly) <- c("x_","y_")
proj4string(dat.mcp.monthly) <- CRS("+init=epsg:4326")

# Bonne equal area projection - https://spatialreference.org/ref/esri/54024/
dat.mcp.monthly <- spTransform(dat.mcp.monthly,sp::CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs"))

mcp.monthly <- mcp(dat.mcp.monthly, percent = 95, unout = c("m2")) %>% data.frame() %>%   left_join(mean.coord, by = c("id" = "id.month")) |> 
  mutate(month = as.numeric(stringr::str_extract(id, "(\\d+$)")),
         year_month = stringr::str_extract(id, "[^.]*$"),
         individual_id = str_extract(id, "[^.]+")) %>% 
  dplyr::select(individual_id,month,year_month,area, mean.x, mean.y)

rm(dat.mcp.monthly);rm(mean.coord)
```

## Intensity of use 

Intensity is defined as the cummulative movement distance during a set time period, divided by the square root of the range size during that same time period. Smaller values 

### Daily Intensity of Use (IoU)

```{r}
mcp.daily$id.day<-paste(mcp.daily$individual_id,mcp.daily$ymd,sep=".")
  
df.IoU24h <- animlocs.1hourly_sl %>% 
  mutate(id.day = paste(individual_id,ymd,sep=".")) %>% 
  group_by(id.day) %>% 
  mutate(cumsum.d1h = sum(d1h,na.rm=T),
         mean.x = mean(x_),
         mean.y = mean(y_)) %>% 
  dplyr::select(individual_id,ymd,cumsum.d1h,mean.x,mean.y) %>% distinct() %>%
  left_join(mcp.daily[,c("id.day","area")],by = "id.day") %>% 
  mutate(iou24h = cumsum.d1h/sqrt(area)) %>% 
  filter(!is.na(iou24h))
```

### Monthly Intensity of Use

```{r class.source = 'fold-hide',message=FALSE,warning=FALSE}
mcp.monthly$id_month = paste(mcp.monthly$individual_id,mcp.monthly$year_month,sep=".")

df.IoU1m <-   animlocs.daily_sl  %>%  
  mutate(year_month = paste(year,month,sep="_"),
         id_month = paste(individual_id,year_month,sep="."))  |>  group_by(id_month)  |>  
  mutate(cumsum.d24h = sum(d24h,na.rm=T),
         mean.x = mean(x_),
         mean.y = mean(y_))  |>  
  dplyr::select(individual_id,id_month,month,year,year_month,cumsum.d24h,mean.x,mean.y)  |>
  distinct() |> 
  left_join(mcp.monthly[,c("id_month","area")],by = "id_month")  |>  
  mutate(iou1m = cumsum.d24h/sqrt(area))  |>  
  filter(!is.na(iou1m)) |> ungroup() |> 
  dplyr::select(individual_id,month,year,year_month,mean.x,mean.y,cumsum.d24h,area,iou1m)  
```

## Diurnality Index

Estimating proportiona daily diurnality, corrected for daylight changes,
using movement data. The diurnality index (from here on diurnality)
is based on Hoogenboom, Daan, Dallinga, and Schoenmakers (1984):

[(AD/DD) - (AN/DN)] / [(AD/DD) + (AN/DN)]

where AD and AN are the sums of the movement during the
day and night, respectively, and DD and DN are the durations of the
day and night, respectively. The diurnality index varies between -1
(night active) and 1 (day active).

```{r}

library(suncalc)

diurn.ind <- animlocs.1hourly_sl |> 
  mutate(lat=y_,lon=x_,date=as.Date(t_)) |> 
  dplyr::select(lat,lon,date,individual_id,t_,ymd,d1h) 

sun_all <- getSunlightTimes(data = diurn.ind, 
                            tz = "UTC",
                            keep = c("sunrise", "sunset"))

diurn.ind$sunrise <- sun_all$sunrise
diurn.ind$sunset<- sun_all$sunset

diurn.ind$daytime <- 
  ifelse(diurn.ind$t_ < diurn.ind$sunrise |
         diurn.ind$t_ > diurn.ind$sunset, "night","day")

      mean.coord <- diurn.ind |> 
        mutate(id.day = paste(individual_id,ymd,sep=".")) |> group_by(id.day) |> 
        mutate(mean.x = mean(lon, na.rm=T),
               mean.y = mean(lat, na.rm=T)) |> 
        dplyr::select(id.day, mean.x, mean.y) |> distinct()

DI <- diurn.ind %>%
  data.frame %>% 
      group_by(individual_id, ymd)%>%
      summarise(dist.sum.day = sum(d1h[daytime=="day"]),
                dist.sum.night = sum(d1h[daytime=="night"]),
                dist.length.day = sum(daytime=="day"),
                dist.length.night = sum(daytime=="night"),
                total.daylength = dist.length.day + dist.length.night)

DI <- data.frame(DI)
```

```{r}
Diurnality <- function(a, b, c, d) {((a / b) - (c / d)) / ((a / b) + (c / d))}
```

We compute Diurnality as an Index according to the above explained 
formula using the "Diurnality" function

a = movement activity during the day   = dist.sum.day
b = length of the day                  = dist.length.day
c = movement activity during the night = dist.length.night
d = length of the night                = dist.length.night

```{r}
DI$diurnality <- Diurnality(DI$dist.sum.day, DI$dist.length.day, DI$dist.sum.night, DI$dist.length.night)

DI <- DI |> 
  mutate(id.day = paste(individual_id,ymd,sep=".")) |> 
  left_join(mean.coord,by = "id.day") |> 
  filter(total.daylength > 19) 
```

# Summarize movement metrics per individual

Use a custom function over each movement metric dataframe to extract individual means/medians/cv/05percentile/95percentile of each metric.

```{r function to summarize 1h Displacements}
FDispl1h<-function(x){
  
  individual_id<-with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per individual
  n1h<-as.numeric(with(x, tapply(x$t_,individual_id, length)))
  
   # 1hr Displacement
  d1h.mean <- as.numeric(with(x, tapply(x$d1h+0.001,individual_id, mean, na.rm=T)))
  d1h.median <- as.numeric(with(x, tapply(x$d1h+0.001,individual_id, median, na.rm=T)))
  d1h.cv <- as.numeric(with(x, tapply(x$d1h+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  d1h.95 <- as.numeric(with(x, tapply(x$d1h+0.001,individual_id, quantile,.95, na.rm=T)))
  d1h.05 <- as.numeric(with(x, tapply(x$d1h+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n1h,
                   d1h.mean,d1h.median,d1h.cv,d1h.95,d1h.05)
  
  return(dats)
}
```

```{r summarize 1h displ}
animlocs.1hourly_sl <- data.frame(animlocs.1hourly_sl)
animlocs.1hourly_sl$individual_id <- factor(animlocs.1hourly_sl$individual_id)

Displ1h <- FDispl1h(animlocs.1hourly_sl[,c("individual_id","t_","d1h")])
rownames(Displ1h) <- NULL
head(Displ1h)
```

Function to summarize Max24h Displacements
```{r function to summarize Max24h Displacements, class.source = 'fold-hide'}
FMaxDispl24h<-function(x){
 individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n.dmax24h.days<-as.numeric(with(x, tapply(x$ymd,individual_id, length)))
  
  # 24hr Displacement
  dmax24h.mean<-as.numeric(with(x, tapply(x$dmax24h+0.001,individual_id, mean, na.rm=T)))
  dmax24h.median<-as.numeric(with(x, tapply(x$dmax24h+0.001,individual_id, median, na.rm=T)))
  dmax24h.cv<-as.numeric(with(x, tapply(x$dmax24h+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  dmax24h.95<-as.numeric(with(x, tapply(x$dmax24h+0.001,individual_id, quantile,.95, na.rm=T)))
  dmax24h.05<-as.numeric(with(x, tapply(x$dmax24h+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.dmax24h.days,
                   dmax24h.mean,dmax24h.median,dmax24h.cv,dmax24h.95,dmax24h.05)
  return(dats)
}

MaxDispl24h <- FMaxDispl24h(dmax24)
rownames(MaxDispl24h) <- NULL
```

Function to summarize 24h Displacements
```{r function to summarize 24h Displacements, class.source = 'fold-hide'}
FDispl24h<-function(x){
  
  individual_id<-with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n24h.days<-as.numeric(with(x, tapply(as.character(x$t_),individual_id, length)))
  
  # 24hr Displacement
  d24h.mean <- as.numeric(with(x, tapply(x$d24h+0.001,individual_id, mean, na.rm=T)))
  d24h.median <- as.numeric(with(x, tapply(x$d24h+0.001,individual_id, median, na.rm=T)))
  d24h.cv <- as.numeric(with(x, tapply(x$d24h+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  d24h.95 <- as.numeric(with(x, tapply(x$d24h+0.001,individual_id, quantile,.95, na.rm=T)))
  d24h.05 <- as.numeric(with(x, tapply(x$d24h+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n24h.days,
                   d24h.mean,d24h.median,d24h.cv,d24h.95,d24h.05)

  return(dats)
}

animlocs.daily_sl <- data.frame(animlocs.daily_sl)
animlocs.daily_sl$individual_id <- factor(animlocs.daily_sl$individual_id)

Displ24h <- FDispl24h(animlocs.daily_sl[,c("individual_id","t_","d24h")])
rownames(Displ24h) <- NULL
```

Function to summarize Max7d Displacements
```{r function to summarize Max7d Displacements, class.source = 'fold-hide'}
FMaxDispl7d<-function(x){
  individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n.dmax7d.weeks<-as.numeric(with(x, tapply(x$year_week,individual_id, length)))
  
  # 24hr Displacement
  dmax7d.mean<-as.numeric(with(x, tapply(x$dmax7d+0.001,individual_id, mean, na.rm=T)))
  dmax7d.median<-as.numeric(with(x, tapply(x$dmax7d+0.001,individual_id, median, na.rm=T)))
  dmax7d.cv<-as.numeric(with(x, tapply(x$dmax7d+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  dmax7d.95<-as.numeric(with(x, tapply(x$dmax7d+0.001,individual_id, quantile,.95, na.rm=T)))
  dmax7d.05<-as.numeric(with(x, tapply(x$dmax7d+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.dmax7d.weeks,
                   dmax7d.mean,dmax7d.median,dmax7d.cv,dmax7d.95,dmax7d.05)

  return(dats)
}

MaxDispl7d <- FMaxDispl7d(dmax7d)
rownames(MaxDispl7d) <- NULL
```

Function to summarize 1d MCP
```{r function to summarize 1d MCP, class.source = 'fold-hide'}
FSumMCP1d<-function(x){
   individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n.mcp24h.days<-as.numeric(with(x, tapply(x$individual_id,individual_id, length)))
  
  # 24MCP Displacement
  mcp24h.mean<-as.numeric(with(x, tapply(x$area+0.001,individual_id, mean, na.rm=T)))
  mcp24h.median<-as.numeric(with(x, tapply(x$area+0.001,individual_id, median, na.rm=T)))
  mcp24h.cv<-as.numeric(with(x, tapply(x$area+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  mcp24h.95<-as.numeric(with(x, tapply(x$area+0.001,individual_id, quantile,.95, na.rm=T)))
  mcp24h.05<-as.numeric(with(x, tapply(x$area+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.mcp24h.days,
                   mcp24h.mean,mcp24h.median,mcp24h.cv,mcp24h.95,mcp24h.05)
  
  return(dats)
}

MCP24h <- FSumMCP1d(mcp.daily)
rownames(MCP24h) <- NULL
```

Function to summarize 7d MCP
```{r function to summarize 7d MCP, class.source = 'fold-hide'}
FSumMCP7d<-function(x){
  individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n.mcp7d.weeks<-as.numeric(with(x, tapply(x$year_week,individual_id, length)))
  
  # 24mcp Displacement
  mcp7d.mean<-as.numeric(with(x, tapply(x$area+0.001,individual_id, mean, na.rm=T)))
  mcp7d.median<-as.numeric(with(x, tapply(x$area+0.001,individual_id, median, na.rm=T)))
  mcp7d.cv<-as.numeric(with(x, tapply(x$area+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  mcp7d.95<-as.numeric(with(x, tapply(x$area+0.001,individual_id, quantile,.95, na.rm=T)))
  mcp7d.05<-as.numeric(with(x, tapply(x$area+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.mcp7d.weeks,
                   mcp7d.mean,mcp7d.median,mcp7d.cv,mcp7d.95,mcp7d.05)
 
  return(dats)
}
MCP7d <- FSumMCP7d(mcp.weekly)
rownames(MCP7d) <- NULL
```

Function to summarize 1m MCP
```{r function to summarize 1m MCP, class.source = 'fold-hide'}
FSumMCP1m<-function(x){
  individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))

  n.mcp1m.months<-as.numeric(with(x, tapply(x$year_month,individual_id, length)))
  
  mcp1m.mean<-as.numeric(with(x, tapply(x$area+0.001,individual_id, mean, na.rm=T)))
  mcp1m.median<-as.numeric(with(x, tapply(x$area+0.001,individual_id, median, na.rm=T)))
  mcp1m.cv<-as.numeric(with(x, tapply(x$area+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  mcp1m.95<-as.numeric(with(x, tapply(x$area+0.001,individual_id, quantile,.95, na.rm=T)))
  mcp1m.05<-as.numeric(with(x, tapply(x$area+0.001,individual_id, quantile,.05, na.rm=T)))
  
  dats<-data.frame(individual_id,n.mcp1m.months,
                   mcp1m.mean,mcp1m.median,mcp1m.cv,mcp1m.95,mcp1m.05)
  
  return(dats)
}
MCP1m <- FSumMCP1m(mcp.monthly)
rownames(MCP1m) <- NULL
```

Function to summarize IoU24h
```{r function to summarize IoU24h, class.source = 'fold-hide'}
FSumIOU24h<-function(x){
 individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n.iou24h.days<-as.numeric(with(x, tapply(x$ymd,individual_id, length)))
  
  # 24hr Displacement
  iou24h.mean<-as.numeric(with(x, tapply(x$iou24h+0.001,individual_id, mean, na.rm=T)))
  iou24h.median<-as.numeric(with(x, tapply(x$iou24h+0.001,individual_id, median, na.rm=T)))
  iou24h.cv<-as.numeric(with(x, tapply(x$iou24h+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  iou24h.95<-as.numeric(with(x, tapply(x$iou24h+0.001,individual_id, quantile,.95, na.rm=T)))
  iou24h.05<-as.numeric(with(x, tapply(x$iou24h+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.iou24h.days,
                   iou24h.mean,iou24h.median,iou24h.cv,iou24h.95,iou24h.05)
  
  return(dats)
}

df.IoU24h <- data.frame(df.IoU24h)
df.IoU24h$individual_id <- factor(df.IoU24h$individual_id)

IOU24h <- FSumIOU24h(df.IoU24h)
rownames(IOU24h) <- NULL
```

Function to summarize IoU1m
```{r function to summarize IoU1m, class.source = 'fold-hide'}
FSumIOU1m<-function(x){
  individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # Get sample size per indivindividual_idual
  n.iou1m.month<-as.numeric(with(x, tapply(x$year_month,individual_id, length)))
  
  # 24hr Displacement
  iou1m.mean<-as.numeric(with(x, tapply(x$iou1m+0.001,individual_id, mean, na.rm=T)))
  iou1m.median<-as.numeric(with(x, tapply(x$iou1m+0.001,individual_id, median, na.rm=T)))
  iou1m.cv<-as.numeric(with(x, tapply(x$iou1m+0.001,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  iou1m.95<-as.numeric(with(x, tapply(x$iou1m+0.001,individual_id, quantile,.95, na.rm=T)))
  iou1m.05<-as.numeric(with(x, tapply(x$iou1m+0.001,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.iou1m.month,
                   iou1m.mean,iou1m.median,iou1m.cv,iou1m.95,iou1m.05)
  
  return(dats)
}

df.IoU1m <- data.frame(df.IoU1m)
df.IoU1m$individual_id <- factor(df.IoU1m$individual_id)

IOU1m <- FSumIOU1m(df.IoU1m)
rownames(IOU1m) <- NULL
```

Function to summarize Diurnality Index
```{r function to summarize IoU1m, class.source = 'fold-hide'}
FSumDI<-function(x){
 individual_id <- with(x, tapply(as.character(x$individual_id),individual_id, unique))
  
  # How many days per indivindividual_idual
  n.di.days<-as.numeric(with(x, tapply(x$ymd,individual_id, length)))
  
  # Diurnality
  di.mean<-as.numeric(with(x, tapply(x$diurnality,individual_id, mean, na.rm=T)))
  di.median<-as.numeric(with(x, tapply(x$diurnality,individual_id, median, na.rm=T)))
  di.cv<-as.numeric(with(x, tapply(x$diurnality,individual_id, function(x) sd(x, na.rm=T) / mean(x, na.rm=T))))
  di.95<-as.numeric(with(x, tapply(x$diurnality,individual_id, quantile,.95, na.rm=T)))
  di.05<-as.numeric(with(x, tapply(x$diurnality,individual_id, quantile,.05, na.rm=T)))
  
  # build dataframe
  dats<-data.frame(individual_id,n.di.days,
                   di.mean,di.median,di.cv,di.95,di.05)
  return(dats)
}

DI <- DI[!is.na(DI$diurnality),]
DI$individual_id <- factor(DI$individual_id)
DI.12m <- FSumDI(DI)
rownames(DI.12m) <- NULL
```

# Build database with summary values

```{r,message=FALSE}
MoveTrait.v0 <- full_join(Displ1h, full_join(Displ24h, full_join(MaxDispl24h, 
             full_join(MaxDispl7d, full_join(MCP24h, full_join(MCP7d, 
             full_join(MCP1m, full_join(IOU24h, full_join(IOU1m, DI.12m)))))))))

movedata2 <- movedata %>% 
  group_by(individual_id) %>% 
  mutate(mean.longitude = mean(Longitude,na.rm=T),
         mean.latitude = mean(Latitude,na.rm=T)) %>% 
  dplyr::select(individual_id,species,animal_mass,contact_person_name,mean.longitude,mean.latitude) %>% 
  dplyr::filter(!duplicated(individual_id))

DBMoveTrait.v0 <- merge(movedata2, MoveTrait.v0, by = "individual_id", all.y=T)
DBMoveTrait.v0 <- DBMoveTrait.v0 %>% droplevels()
```

```{r}
saveRDS(DBMoveTrait.v0,"./DATA/Tucker/MoveTraitsDB.v0.1_Tucker.rds")
```


# Build full database including raw data

```{r,warning=FALSE,message=FALSE}
# Create the full data set with hourly and daily relocation intervals
MoveTrait.v0_spatial <- 
  DBMoveTrait.v0  |> 
  tidyr::nest(data = -individual_id)|> 
  dplyr::select(-data) |>  # just a quick trick to keep things flowing.
  left_join(DBMoveTrait.v0[!duplicated(DBMoveTrait.v0$individual_id),1:66], 
            by = c("individual_id" = "individual_id"))  |>  
  # 1 hourly
  left_join(animlocs.1hourly_sl %>%
              data.frame %>% 
              dplyr::select("individual_id","t_","d1h","x_","y_") %>% 
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(Displ.1h = data) |>  
  
  # 24 hourly
  left_join(animlocs.daily_sl %>%
              dplyr::select("individual_id","t_","d24h","x_","y_") %>% 
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(Displ.24h = data) |>  
  
  # Dmax24
  left_join(dmax24 %>%
              dplyr::select("individual_id","ymd","dmax24h","mean.x", "mean.y") %>% 
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(MaxDispl.24h = data)  |>  
  
  # Dmax7d
  left_join(dmax7d %>%
              dplyr::select("individual_id","week","year_week","dmax7d","mean.x", "mean.y") %>% #,"week","year_week",s"mean.x", "mean.y"
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(MaxDispl.7d = data)  |>  
  
  # mcp.daily
  left_join(mcp.daily %>%
              dplyr::select("individual_id","ymd","area","mean.x", "mean.y") %>% 
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(Mcp.24h = data) |>  
  
  # mcp.weekly
  left_join(mcp.weekly %>%
              dplyr::select("individual_id","week","year_week","area","mean.x", "mean.y") %>% #,"week","year_week",mean.x", "mean.y"
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(Mcp.7d = data)  |>  
  
  # mcp.monthly
  left_join(mcp.monthly %>%
              dplyr::select("individual_id","month","year_month","area","mean.x", "mean.y") %>%  #,"month","year_month"
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(Mcp.1m = data) |> 
  
  # Intensity of Use 24h
  left_join(df.IoU24h %>%
              ungroup() |>
              dplyr::select("individual_id","ymd","iou24h","mean.x", "mean.y") %>% 
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(IoU.24h = data)  |>  
  
  # Intensity of Use 1m
  left_join(df.IoU1m %>%
              ungroup() |>
              dplyr::select("individual_id","month","year","iou1m","mean.x", "mean.y") %>% ## month, year_month
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(IoU.1m = data)  |>  
  
# Diurnality Index
  left_join(DI %>%
              dplyr::select("individual_id","ymd","diurnality","mean.x","mean.y") %>% 
              tidyr::nest(data = -individual_id), 
            by = c("individual_id" = "individual_id")) %>% 
  dplyr::rename(diurnality = data) 

```

```{r}
saveRDS(MoveTrait.v0_spatial,
        "./DATA/Tucker/MoveTraitsDB.v0.1_spatial_Tucker.rds")  

```