---
title: "moveHMM"
author: "Kadlec Ivo"
date: "18 2 2021"
output:
  html_document: default
  word_document: default
  pdf_document: default
---

# moveHMM - použití HMM na samici geparda štíhlého

Data jsou stažena z https://www.movebank.org/. Data pochází ze GPS sledování samice reintrodukovaného geparda žijícího v Jižní Africe v národním parku Pilanesberg. Data se začali sbírat v roce 2014, kdy samice neměla žádné mláďata a ani nebyla březí (v datech jako cheetah_1). V listopadu toho stejného roku zabřezla (Cheetah_2) a v roce 2015 porodila 4 mláďata samce a jedno za chvíli umřelo (Cheetah_3).

## Instalce balíčku, načtení dat a jejich úprava

```{r, warning=FALSE, message=FALSE}
library(moveHMM)
```
Vstupní data musí mít správný formát, aby bylo možné data zpracovat a poté je analyzovat. Data-frame musí obsahovat 2 sloupce:

  - sloupec s názvem "x": Easting nebo zeměpisná délka (longitude)

  - sloupec s názvem "y": Northing nebo zeměpisná šířka (latitude)

Důležitým sloupec je “ID”. Pokud není ve vstupních datech žádný sloupec s názvem “ID”, tak balíček bude předpokládat, že daná data jsou pouze pro jednoho jedince.

```{r}
gepard <- read.table("Gepard_1.csv", header = T, sep = ";", dec=",")
gepard$time = as.POSIXct(gepard$time)

gepard_D <- read.table("Gepard_2.csv", header = T, sep = ";", dec=",")
gepard_D$time = as.POSIXct(gepard_D$time)

gepard_PD <- read.table("Gepard_3.csv", header = T, sep = ";", dec=",")
gepard_PD$time = as.POSIXct(gepard_PD$time)
```
Převod do souřadnicového systému UTM:

```{r, results='hide', message=FALSE, warning=FALSE}
library(sf)

floor((180+gepard$long[1])/6) + 1
gepard_sf <- st_as_sf(gepard, coords = c("long", "lat"), crs = 4326)
gepard_utm <- st_transform(gepard_sf, crs = 32635)
gepard <- cbind(gepard_utm, st_coordinates(gepard_utm))

gepard_D_sf <- st_as_sf(gepard_D, coords = c("long", "lat"), crs = 4326)
gepard_D_utm <- st_transform(gepard_D_sf, crs = 32635)
gepard_D <- cbind(gepard_D_utm, st_coordinates(gepard_D_utm))

gepard_PD_sf <- st_as_sf(gepard_PD, coords = c("long", "lat"), crs = 4326)
gepard_PD_utm <- st_transform(gepard_PD_sf, crs = 32635)
gepard_PD <- cbind(gepard_PD_utm, st_coordinates(gepard_PD_utm))
```

Náhled na vstupní data:

```{r, echo=FALSE}
head(gepard)
```
Data obsahují jak sloupce se souřadnicemi, tak i sloupce s možnými kovariátami:
 
  - TEMP = venkovní teplota (°C)
 
  - Height = nadmořská výška (m n. m.)
  
Navíc je tady i sloupec ID (označení "jedince") a čas sběru dat (time)

V datech lze vidět, že některá měření chybí a je potřeba je doplnit, tak aby 
bylo každé měření po 6 hodinách.

```{r}
gepard_1a=data.frame(time=seq.POSIXt(from = min(gepard$time), by= "6 hour",length.out = 788))

gepard_1b=data.frame(time=seq.POSIXt(from = gepard$time[788],to=max(gepard$time), by= "6 hour"))

geparad_data1=rbind(gepard_1a,gepard_1b)
gepard = merge(geparad_data1,gepard,by="time",all.x = TRUE)
```

```{r}
#gepard_D (březí)
gepard_2a=data.frame(time=seq.POSIXt(from = min(gepard_D$time), to =max(gepard_D$time), by= "6 hour"))

gepard_D = merge(gepard_2a,gepard_D,by="time",all.x = TRUE)

#gepard_PD (s mláďaty)
gepard_3a=data.frame(time=seq.POSIXt(from = min(gepard_PD$time), to =max(gepard_PD$time), by= "6 hour"))

gepard_PD = merge(gepard_3a,gepard_PD,by="time",all.x = TRUE)
```
Hodnoty “x” a “y” jsou v datech vyjádřeny v metrech, proto je potřeba data přetransformovat, aby byla délka v kilometrech.

```{r}
gepard$X <- gepard$X/1000
gepard$Y <- gepard$Y/1000

gepard_D$X <- gepard_D$X/1000
gepard_D$Y <- gepard_D$Y/1000

gepard_PD$X <- gepard_PD$X/1000
gepard_PD$Y <- gepard_PD$Y/1000
```

## Model pro geparda po vypuštění

Před použití následující funkce „PrebData“ je potřeba pro správné fungovaní baličku spustit následující funkci: 

```{r}
get_closest_id <- function(i, ids=gepard$ID){
  id <- ids[i]
  ii <- i
  while (is.na(id) & ii > 1){
    ii <- ii - 1
    id <- ids[ii]
  }
  ii <- i
  while (is.na(id) & ii < length(ids)){
    ii <- ii + 1
    id <- ids[ii]
  }
  return(id)
}

for (i in which(is.na(gepard$ID))) {
  gepard$ID[i] <- get_closest_id(i)
}
```
Funkce "PrebData" umožňuje vymodelovat délky kroků a jejich úhel otočení a je 
vstupním objektem pro tvorbu modelu.

Duležité je zvolit správné argumenty funkce:

  - type: určuje, zda-li jsou souřadnice v easting/northing (type="UTM"), a nebo
  v zeměpisné šířce/délce (type="LL"), který je výchozí

  - coordNames: pouřívá se v případě, když jsou sloupce se souřadnicemi nejsou
  popsané "x" a "x".

Funkce "PrebData":
```{r}
data <- prepData(gepard,type="UTM",coordNames=c("X","Y"))
```
Varování upozorňuje na data bez kovariáty a říká, že nahradí nejbližší hodnotou.

```{r, echo=FALSE}
head(data)
```
Po použití funkce summary lze vidět v hlavičce, kolik je  sledování pro daného
jedince a ve spodní části je zobrazeno minimum, medián, průměr kovariát 

```{r, message=FALSE, error=FALSE, message=FALSE}
library(dplyr)
data %>% select(-c(geometry)) %>% 
  summary()
```
### Vizualizace "prebData"
"plot("prebData)" zobrazí:

  - mapu trajektorie jedinců

  - časové řady kroků a úhlů

  - histogramy kroků a úhlů

```{r}
plot(data, compact = T)
```

### Model fitHMM

Fit HMM je funkce, která vytvoří model pro skryté Markovovy modely, používající numericko-logarimtickou funkci.

Základní argumenty pro model:
  
  - data: moveData
  
  - nbStates: počet stavů HMM
  
  - stepPar0: počáteční vektor délky
  
  - anglePar0: počáteční vektor úhlů
  
  - formula: regresní formula pro kovariáty (výchozí: 1)
  
  - stepDist: název rozdělení pro délku kroků (výchozí: gama)
  
  - angleDist: název rozdělní pro úhel (výchozí: vm)
  
  - angleMean: vektro úhlů (výchozí: NULL)
  
  - stationary: argument FALSE použije se pokud existují kovariáty (výchozí: FALSE)

Dvoustavový model:
Parametry pro délku kroků:
  
  - m0 = průměrná délka kroku
  
  - sigma0 = odchylka kroku
  
  - zeromass0 - použiju pokud se v datech nachází 0 délka kroku
  
  - StepPa0 = c(mu,sigma0,zeromass0)
  
```{r}
mu0 <- c(0.5,2)
sigma0 <- c(0.5,1)

#kontrola jestli je potřeba zero-mass
whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
#vyšlo číslo 0.02244389, které říká jaký je poměr 0 hodnot v datech
zeromass0 <- c(0.1,0.05) 
stepPar0 <- c(mu0,sigma0,zeromass0)
```
Parametry pro velikost úhlů:
  
  - angleMean0 = průměrná velikost úhlů
  
  - kappa0 = odchylka pro velikost úhlu
  
  - anglePar0 <- c(angleMean0, kappa0)
  
```{r}
angleMean0 <- c(pi,0)
kappa0 <- c(0.7,1.5)
anglePar0 <- c(angleMean0, kappa0)
```
Zvolené rozdělení:
  
  - stepDist - rozdělení pro délku kroku ("gamma" = gamma rozdělení)
  
  - angleDist - rozdělení pro úhly otáčení ("vm" = von Misesovo rozdělení)
  
```{r}
stepDist <- "gamma"
angleDist <- "vm"
```

```{r}
m1 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~1, stepDist=stepDist,
             angleDist =angleDist)#žádná kovariata
m2 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP, stepDist=stepDist,
             angleDist =angleDist)#kovariata -> teplota
m3 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~Height, stepDist=stepDist,
             angleDist =angleDist)#"kovariata -> nadmořská výška
m4 <- fitHMM(data=data, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP+Height, stepDist=stepDist,
             angleDist =angleDist)# kovarita TEMP + Height
```

Tří stavový model:
```{r}
mu0 <- c(0.5,2,4)
sigma0 <- c(0.5,1,1.5)

whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
zeromass0 <- c(0.1,0.05,0.01) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi,pi/2,0)
kappa0 <- c(0.7,1.5,2)
anglePar0 <- c(angleMean0, kappa0)

m1_3 <- fitHMM(data=data, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~1, stepDist=stepDist,
             angleDist =angleDist)# žádná kovariata
m3_3 <- fitHMM(data=data, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~Height, stepDist=stepDist,
             angleDist =angleDist)# kovariata -> nadmořská výška
m4_3 <- fitHMM(data=data, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP+Height, stepDist=stepDist,
             angleDist =angleDist)# kovarita TEMP + Height
```

Model s jedním stavem:
```{r}
mu0 <- c(2)
sigma0 <- c(1)

whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
zeromass0 <- c(0.05) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi/2)
kappa0 <- c(1.5)
anglePar0 <- c(angleMean0, kappa0)

m0 <- fitHMM(data=data, nbStates = 1, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~1, stepDist=stepDist,
             angleDist =angleDist)
```

Test AIC
```{r}
print(AIC(m1,m2,m3,m4,m1_3,m3_3,m4_3,m0))
```


Jako nejlepší model byl vybrán nakonec m2, protože má ze všech modelů 
nejnižší hodnotu AIC

```{r}
m2
```

### Konfidenční intervaly

```{r}
CI(m2)
```

V tabulce CI lze vidět:

  - interval po střední délku kroku -> pro stav 1 = (0.223207117-0.359104588) m; pro stav 2 = (1.51905094-2.041851) m
  
  - interval pro variabliní délku kroku -> pro stav 1 = 0.359104588 m; pro stav 2 = 1.12934926 m
                           
  - střední velikost úhlu -> pro stav 1 = (-5.35580872 do -0.8082987) 
                          -> pro stav 2 = (-0.2698473 -0.2968336)
                          
  - koncetrace úhlů -> pro stav 1 = 0.1948008  
  
                    -> pro stav 2 = 0.61492293 
  
  - pravděpodobnost přechodu z jednoho stavu do druhého a na zpět

### Vizualizace modelu

```{r}
plot(m2, plotCI = T)
```

### Pravděpodobnosti stavu

```{r, message=FALSE}
sp <- stateProbs(m2)
head(sp)
plotStates(m2)
```

### Pravděpodobnost stacionárncího stavu

```{r}
plotStationary(m2, plotCI=TRUE)
```
Stacionární stav bude nejpravděpodobnější mezi 15-20°C. 

### Kontrola modelu

```{r, message=FALSE}
res <- pseudoRes(m2)
plotPR(m2)
```

Test normality

```{r, message=FALSE}
library(tseries)
#test normality kroků
## Jarqueův a Berryho test normality

print(jarque.bera.test(res$step[which(!is.na(res$step))]))
## Shapirro-Wilkův test normality
shapiro.test(res$stepRes)

#test normality pro úhly
print(jarque.bera.test(res$angle[which(!is.na(res$angle))]))
shapiro.test(res$angleRes)
```

Test korelace

```{r}
par(mfrow=c(1,2))
acf(na.omit((data$step)))
acf(na.omit((data$angle)))
```

## Model pro březího geparda
```{r}
get_closest_id <- function(i, ids=gepard_D$ID){
  id <- ids[i]
  ii <- i
  while (is.na(id) & ii > 1){
    ii <- ii - 1
    id <- ids[ii]
  }
  ii <- i
  while (is.na(id) & ii < length(ids)){
    ii <- ii + 1
    id <- ids[ii]
  }
  return(id)
}

for (i in which(is.na(gepard_D$ID))) {
  gepard_D$ID[i] <- get_closest_id(i)
}
```

```{r}
data_D <- prepData(gepard_D,type="UTM",coordNames=c("X","Y"))
```

```{r}
plot(data_D, compact = T)
```

Model:
```{r}
# Parametry 2 stavy
mu0 <- c(0.2,1)
sigma0 <- c(0.5,1)

whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
zeromass0 <- c(0.1,0.05) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi,0)
kappa0 <- c(0.7,1.5)
anglePar0 <- c(angleMean0, kappa0)

stepDist <- "gamma"
angleDist <- "vm"

# Model 2 stavy
m2_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~TEMP, stepDist=stepDist,
               angleDist =angleDist)
m1_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~Height, stepDist=stepDist,
               angleDist =angleDist)
m3_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP+Height, stepDist=stepDist,
             angleDist =angleDist)# kovarita TEMP + Height
m4_D <- fitHMM(data=data_D, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~1, stepDist=stepDist,
               angleDist =angleDist)

# Parametry 3 stavy
mu0 <- c(0.5,2,4)
sigma0 <- c(0.5,1,1.5)

whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
zeromass0 <- c(0.1,0.05,0.01) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi,pi/2,0)
kappa0 <- c(0.7,1.5,2)
anglePar0 <- c(angleMean0, kappa0)

m1_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~1, stepDist=stepDist,
             angleDist =angleDist)#žádná kovariata
m3_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~Height, stepDist=stepDist,
             angleDist =angleDist)#"kovariata -> nadmořská výška
m2_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP, stepDist=stepDist,
             angleDist =angleDist)#"kovariata -> nadmořská výška
m4_3_D <- fitHMM(data=data_D, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP+Height, stepDist=stepDist,
             angleDist =angleDist)# kovarita TEMP + Height

# Parametry 1 stav
mu0 <- c(1)
sigma0 <- c(0.5)

whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
zeromass0 <- c(0.1) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi)
kappa0 <- c(0.7)
anglePar0 <- c(angleMean0, kappa0)

stepDist <- "gamma"
angleDist <- "vm"

# Model 1 stav
m0_D <- fitHMM(data=data_D, nbStates = 1, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~TEMP, stepDist=stepDist,
               angleDist =angleDist)
```
AIC

```{r}
print(AIC(m2_D,m0_D,m1_D,m4_D, m1_3_D,m2_3_D,m3_3_D, m4_3_D))
```

```{r}
m4_D
```

### Konfidenční intervaly

```{r}
CI(m4_D)
```

### Vizualizace modelu

```{r}
plot(m4_D, plotCI = T)
```

### Pravděpodobnosti stavu

```{r}
sp <- stateProbs(m4_D)
head(sp)
plotStates(m4_D)
```

### Kontrola modelu

```{r}
res <- pseudoRes(m4_D)
plotPR(m4_D)
```

Test normality

```{r, message=FALSE}
library(tseries)
#test normality kroků
## Jarqueův a Berryho test normality
print(jarque.bera.test(res$step[which(!is.na(res$step))]))
## Shapirro-Wilkův test normality
shapiro.test(res$stepRes)

#test normality pro úhly
print(jarque.bera.test(res$angle[which(!is.na(res$angle))]))
shapiro.test(res$angleRes)
```

Test korelace

```{r}
par(mfrow=c(1,2))
acf(na.omit((data_D$step)))
acf(na.omit((data_D$angle)))

```

## Model pro geparda s mláďaty v poporodní fázi

```{r}
get_closest_id <- function(i, ids=gepard_PD$ID){
  id <- ids[i]
  ii <- i
  while (is.na(id) & ii > 1){
    ii <- ii - 1
    id <- ids[ii]
  }
  ii <- i
  while (is.na(id) & ii < length(ids)){
    ii <- ii + 1
    id <- ids[ii]
  }
  return(id)
}

for (i in which(is.na(gepard_PD$ID))) {
  gepard_PD$ID[i] <- get_closest_id(i)
}
```

```{r}
data_PD <- prepData(gepard_PD,type="UTM",coordNames=c("X","Y"))
```

```{r}
plot(data_PD, compact = T)
```
Model:

```{r}
# Parametry 2 stavy
mu0 <- c(1,2)
sigma0 <- c(0.5,1)

whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
zeromass0 <- c(0.1,0.05) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi,0)
kappa0 <- c(0.7,1.5)
anglePar0 <- c(angleMean0, kappa0)

stepDist <- "gamma"
angleDist <- "vm"

# Model 2 stavy
m2_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP, stepDist=stepDist,
             angleDist =angleDist)
m1_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~Height, stepDist=stepDist,
               angleDist =angleDist)
m3_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~Height+TEMP, stepDist=stepDist,
               angleDist =angleDist)
m4_PD <- fitHMM(data=data_PD, nbStates = 2, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~1, stepDist=stepDist,
               angleDist =angleDist)

# Parametry 3 stavy
mu0 <- c(0.5,2,4)
sigma0 <- c(0.5,1,1.5)

whichzero <- which(data$step == 0)
length(whichzero)/nrow(data)
zeromass0 <- c(0.1,0.05,0.01) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi,pi/2,0)
kappa0 <- c(0.7,1.5,2)
anglePar0 <- c(angleMean0, kappa0)

# Model 3 stavy
m1_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~1, stepDist=stepDist,
             angleDist =angleDist)#žádná kovariata
m3_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~Height, stepDist=stepDist,
             angleDist =angleDist)#"kovariata -> nadmořská výška
m2_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP, stepDist=stepDist,
             angleDist =angleDist)#"kovariata -> teplota
m4_3_PD <- fitHMM(data=data_PD, nbStates = 3, stepPar0 = stepPar0, anglePar0 = anglePar0,
             formula = ~TEMP+Height, stepDist=stepDist,
             angleDist =angleDist)# kovarita TEMP + Height

# Parametry 1 stav
mu0 <- c(0.5)
sigma0 <- c(0.5)

whichzero <- which(data_D$step == 0)
length(whichzero)/nrow(data_D)
zeromass0 <- c(0.1) 
stepPar0 <- c(mu0,sigma0,zeromass0)

angleMean0 <- c(pi)
kappa0 <- c(0.7)
anglePar0 <- c(angleMean0, kappa0)

stepDist <- "gamma"
angleDist <- "vm"

# Model 1 stav
m0_PD <- fitHMM(data=data_PD, nbStates = 1, stepPar0 = stepPar0, anglePar0 = anglePar0,
               formula = ~TEMP, stepDist=stepDist,
               angleDist =angleDist)
```
AIC

```{r}
print(AIC(m2_PD,m1_PD,m4_PD,m0_PD,m1_3_PD,m2_3_PD,m3_3_PD,m4_3_PD))
```

```{r}
m2_PD
```

### Konfidenční intervaly

```{r}
CI(m2_PD)
```

### Vizualizace modelu

```{r}
plot(m2_PD, plotCI = T)
```

### Pravděpodobnosti stavu

```{r}
sp <- stateProbs(m2_PD)
head(sp)
plotStates(m2_PD)
```

### Pravděpodobnost stacionárního stavu

```{r}
plotStationary(m2_PD, plotCI=TRUE)
```

### Kontrola modelu

```{r}
res <- pseudoRes(m2_PD)
plotPR(m2_PD)
```

Test normality

```{r, message=FALSE}
library(tseries)
#test normality kroků
print(jarque.bera.test(res$step[which(!is.na(res$step))]))
## Shapirro-Wilkův test normality
shapiro.test(res$stepRes)

#test normality pro úhly
print(jarque.bera.test(res$angle[which(!is.na(res$angle))]))
shapiro.test(res$angleRes)
```

Test korelace

```{r}
par(mfrow=c(1,2))
acf(na.omit((data_PD$step)))
acf(na.omit((data_PD$angle)))
```

## Použití balíčku trajr na vyhlazení trajektorie a její nasledné zobrazení pomocí leaflet

```{r, warning=FALSE}
library(trajr)
trj1<-read.table("Cheetah_1.csv", dec=",",sep = ";", head=T)
trj1$time=as.Date(trj1$time)
trj2<-read.table("Cheetah_2.csv", dec=",",sep = ";", head=T)
trj2$time=as.Date(trj2$time)
trj3<-read.table("Cheetah_3.csv", dec=",",sep = ";", head=T)
trj3$time=as.Date(trj3$time)

Ch_1 <- TrajFromCoords(trj1, xCol = "x", yCol = "y", timeCol = "time")
Ch_2 <- TrajFromCoords(trj2, xCol = "x", yCol = "y", timeCol = "time")
Ch_3 <- TrajFromCoords(trj3, xCol = "x", yCol = "y", timeCol = "time")

p=9
n=p + 3 - p%%2
Cheetah_1 <- TrajSmoothSG(Ch_1, p, n)
Cheetah_2 <- TrajSmoothSG(Ch_2, p, n)
Cheetah_3 <- TrajSmoothSG(Ch_3, p, n)

library(leaflet)

map = leaflet() %>%
  addTiles() %>%
  addPolylines(data = Cheetah_1, lat = ~y, lng = ~x, color = "orange",weight =2, group = )%>%
  addPolylines(data = Cheetah_2, lat = ~y, lng = ~x, color = "blue", weight= 2, group = )%>%
  addPolylines(data = Cheetah_3, lat = ~y, lng = ~x, color = "green", weight = 2, group = )%>%
  addLegend(position = "topright", labels = c("Samice po vypuštění", "Těhotná samice", "Samice s mláďaty"),
            colors = c("orange","blue","green"))
map  
```