datos_telemetricos.qmd

---
format:
  html:
    toc: true
    toc-title: "Contenidos"
    toc-expand: 1
    code-overflow: wrap
editor_options:
  markdown:
    wrap: sentence
    canonical: true
lang: es
---
```{r}
#| include: false
source("_common.R")
```

# Datos de telemetría satelital {#sec-datos-telemetricos}

Como he hecho énfasis hasta ahora, antes de iniciar con análisis, es necesario asegurarse que la información ingresada sea la apropiada para evitar cualquier tipo de errores, desde mensajes de advertencia, a impedimentos para correr funciones, y errores en los resultados finales.

## Importación, limpieza y visualización de datos telemétricos

Es de suma importancia que estos pasos sean realizados de manera correcta, ya que cualquier error en la importación y limpieza de los datos nos puede traer grandes errores.

A manera de ejercicio, he tomado material del curso [Landscape and Analysis Modeling](https://mltconsecol.github.io/TU_LandscapeAnalysis_Documents/index.html) de The University of Tulsa, al que pueden acceder gratuitamente para profundizar más en los temas que trataremos en estas clases. Además, modifiqué el formato del archivo original `MigratoryZebra.csv`, por un motivo que será obvio a continuación.

Es siempre aconsejable tener una mirada general de los datos antes de leerlos como objetos en R, ya que se debe utilizar la función adecuada para su lectura. [Descarga](https://mega.nz/folder/VwVQEbJK#lxYXI88_90s6hUihhwNXgg) y abre el archivo que usarás en este módulo `zebras.csv`, en un archivo de texto y observa algunas características importantes del archivo como los separadores de columnas y los indicadores de decimales.

:::::: {.callout-caution}
## **IMPORTANTE**
Evita abrir cualquier archivo que será importado a R, utilizando Excel como primera opción. Utiliza editores de texto como Notepad o Notepad++.
:::

![Captura de pantalla de zebras.csv](docs/Figures/zebras_csv_screenshot.png){#fig-captura .column-page-right .border fig-alt="Captura de pantalla del archivo zebras.csv"}

Puedes ver que las columnas están separadas por `;`, y los decimales mediante `,`.

```{r eval = FALSE}
# Intenta importar el archivo y mira las características de estos objetos mediante summary() y str()

data1 <- read.csv("zebras.csv")
data2 <- read.csv2("zebras.csv")
data3 <- read.table("zebras.csv", header = TRUE, na.strings = "NA",  sep=";")
```

Ahora activa el paquete _readr_ (también parte de _tidyverse_) e intenta importar los datos mediante la función `read_csv2`.

```{r include=FALSE}
library(readr)
data <- read_csv2("docs/datos/zebras.csv")
```
```{r eval = FALSE}
data <- read_csv2("zebras.csv")
summary(data)
str(data)
```

Ya que tienes el archivo debidamente importado, extrae algunas columnas de interés y asigna un nombre más intuitivo, fácil de escribir y recordar. Esta es una preferencia personal, así que puedes omitir esta parte según tu preferencia.

```{r}
zebras <- data[, c("event-id", "individual-local-identifier",
                "location-long", "location-lat", "study-local-timestamp")]

zebras <- zebras %>%
              rename(id = "event-id", identifier = "individual-local-identifier",
              long = "location-long",
              lat = "location-lat",
              timestamp = "study-local-timestamp")

# Inspecciona las primeras filas y los nombres de las columnas
head(zebras)
```

Ahora puedes observar la distribución de estas zebras de una manera sencilla mediante funciones base.

::: {#fig-todas-zebras-plot fig-cap="Distribución espacial de todos los individuos en el data frame `zebras`"}
```{r}
# Cambia el valor de pch y mira que sucede con el gráfico
plot(zebras[, c("long", "lat")], pch = 20)
```
:::

Los individuos muestreados se encuentran en una extensión de terreno enorme, por lo que es preferible para los fines de esta primera lección, que analices animales que comparten un mismo espacio geográfico. Para esto, puedes hacer un subset del data frame original, y además puedes practicar tus habilidades en manejar _dyplyr_.

```{r}
dos_zebras <- zebras %>%
          filter(identifier == "Z3864" | identifier == "Z6405") %>%
          dplyr::select(identifier, long, lat, timestamp)
```

Realiza nuevamente un gráfico asignando colores a cada individuo.

::: {#fig-dos-zebras-plot fig-cap="Distribución espacial de dos zebras extraídas del data frame `zebras`"}
```{r}
#cambia el valor de pch y mira que sucede con el gráfico
plot(dos_zebras[, c("long", "lat")], pch = 20, col = c("#440D54", "#3CBB75"))
legend("topright", legend = c("Z3864", "Z6405"),
      fill = c("#440D54", "#3CBB75"),
      pch = 20, box.lty = 0)
```
:::

:::::: {.callout-note}
## **Tarea**
Utiliza tus habilidades y conocimiento en `ggplot2` para realizar un gráfico similar al que acabas de realizar con `plot`.
:::

Como has aprendido en lecciones anteriores, no es conveniente embarcarse en algúm tipo de análisis sin antes realizar una limpieza de tus datos. Ahora puedes proceder a remover los NA’s de las columnas señaladas por `summary`.

```{r}
# Investiga si el set de datos posee NAs
summary(dos_zebras)
which(is.na(dos_zebras$lat))
which(is.na(dos_zebras$long))
```

:::::: {.callout-tip}
## Ejercicio
Utiliza cualquier método para eliminar los datos de longitud y latitud faltantes en la fila 1802, y verifica que el proceso de limpieza de datos sea correcto.
:::
```{r include=FALSE}
# Investiga si el set de datos posee NAs
dos_zebras <- dos_zebras %>% na.omit
```

!Haz limpiado los datos con éxito! Ahora crea un objeto que será transformado y proyectado posteriormente.

```{r}
zebras.proj <- dos_zebras
```

## Transformación y proyección cartográfica {#sec-transf-proy}

Para continuar y entender cual es el propósito de los siguientes pasos, primero debes entender que es una proyección cartográfica. En breve, la posición en el espacio tridimensional de un animal tomada mediante coordenadas decimales deben ser transformadas a un objeto bidimensional. Esta transformación viene acompañada de un gran problema como lo es la distorsión, y esta será más o menos evidente según el tipo de proyección que elijas.

El canal de YouTube Vox tiene un excelente video que muestra gráficamente este problema.
{{< video https://www.youtube.com/watch?v=kIID5FDi2JQ&ab_channel=Vox >}}

El paquete _sp_ es el indicado para esta tarea ya que ofrece clases y métodos para manejar datos espaciales. Así, empezarás por transformar el objeto `zebras.proj` a un objeto de clase `SpatialPointsDataFrame`.

```{r warnings = FALSE}
library(sp)
library(sf)

zebras.proj <- SpatialPointsDataFrame(coords = as.data.frame(cbind(zebras.proj$long, zebras.proj$lat)),
                data = zebras.proj,
                proj4string = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs"))
```

Para entender el código que antecede a la función `proj4string` primero te debes preguntar, ¿qué significa CRS, y qué argumentos contiene esta función? *proj4string* no es nada más que una librería que dispone de métodos que permiten la transformación entre diferentes sistemas de coordenadas de referencia **(CRS)**.

`CRS` permite determinar la ubicación de un punto de una manera estandarizada, utilizando uno o más números que representan tanto su posición vertical como horizontal. Para lograr esto, utiliza como referencia un elipsoide que es una representación matemática de la Tierra.
Una vez definido el sistema de coordenadas, `datum` define la posición y orientación del elipsoide de referencia en relación con el centro de la Tierra y el meridiano usado como longitud cero, el meridiano principal. Investiga sobre el significado de **WGS84**.

Como ya deberia ser habitual para ti, analiza tu objeto antes de proseguir. ¿Qué sucedió con `zebras.proj`? Aplica funciones que has aprendido para analizar este nueva clase de objeto.

:::::: {.callout-caution}
He decidido mostrar el código antiguo de definir la proyección para mejorar la comprensión de lo que sucede tras escenas del código. Sin embargo `proj4string` ahora debe definirse como `CRS(SRS_string = "EPSG:4326")`
:::


```{r}
class(zebras.proj)
summary(zebras.proj)
```

Ya has verificado la transformación de las CRS, sin embargo, estos siguen siendo puntos en un mapa representados por medidas angulares. ¿De qué manera podemos medir la distancia entre ellos para poder calcular el área descrita por los puntos posteriormente? Transformando esta proyección a un sistema basado en metros.

```{r}
zebras.transf <- spTransform(zebras.proj,
                  CRS("+proj=lcc +lat_1=20 +lat_2=-23 +lat_0=0 +lon_0=25 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs"))
```

Esta CRS es una proyección cónica y representa a la forma más precisamente que a un área. Esta transformación nos da como resultado un mapa matemático de coordenadas esféricas en un plano bidimensional, cuyas unidades pueden son expresadas en metros. Con todo esto en mente, ya podemos plotear este objeto proyectado de clase `SpatialPointsDataFrame`.

::: {#fig-zebras-proyectadas fig-cap="Zebras extraídas del data frame `zebra` observadas mediante una proyección cónica y transformada a un sistema basado en metros"}

```{r}
# ¿Qué diferencia tiene este grafico con el anterior?
plot(zebras.transf[, c("long", "lat")],
        pch = 20, col = c("#440D54", "#3CBB75"), axes=TRUE)
```
:::

¿Cuál es la diferencia entre la @fig-zebras-proyectadas y @fig-dos-zebras-plot ?

## Análisis de trayectorias

Debes tener en cuenta que la limpieza de los sets de datos no es el único paso a seguir en el análisis espacial, también debes explorar métodos de segmentación o partición para poder discernir cambios en la conducta o trayectoria de un animal (@gurarie_what_2016). Según la pregunta de investigación que tengas debes indagar más en otros métodos como First Passage Time (FPT), step selection function, etc.

Aquí, verás brevemente como puedes obtener información sobre la trayectoria de las zebras que extrajimos previamente.

## _move2_

Este paquete es la continuación de _move_, y asi como _sp_, realizó una transición para incorporar funciones de _sf_. Para iniciar necesitarás convertir el data frame `zebras` a un objeto de clase `move2`, y además deberás ocuparte de otro problema, la transformación del tiempo.

```{r}
dos_zebras$timestamp <- as.POSIXct(dos_zebras$timestamp, "%Y-%m-%d %H:%M:%S", tz="UTC")
```

:::::: {.callout-important}
Establecimos la zona horaria como UTC, que es la que usualmente usan como predeterminada los collares satelitales. Sin embargo deber ser cauteloso al momento de definir la zona horaria, y transformarla a tu zona de estudio, si es que fuese necesario. Esto lo puedes hacer con el paquete `lubridate`.
:::

Mediante `as.POSIXct`, te has asegurado de transformar una columna que fue importada como caracter, y ahora la tienes en formato `POSIXct`, con el cual R maneja el tiempo. Una vez tengas tu set de datos con el formato adecuado, puedes empezar.

```{r warnings = FALSE}
library(move2)
library(units)

dos_zebras_move <- mt_as_move2(dos_zebras, coords = c("long", "lat"),
            time_column = "timestamp", track_id_column = "identifier") %>%
            # en Mercator porque funciones no corren en Lambert
            sf::st_set_crs(4326)

# Estadísticas del objeto
mt_n_tracks(dos_zebras_move)

# ¿Cada cuantas horas se tomaron posiciones satelitales?
timeLags <- set_units(mt_time_lags(dos_zebras_move), hours)

# Distribucion de time lags
timeLags_h <- units::drop_units(timeLags)
hist(timeLags_h, breaks = 50, main=NA, xlab = "Time lag en horas")
```

Puedes graficar la distribucion de puntos por individuo.
```{r warnings = FALSE}
# Grafico por individuos con ggplot2
ggplot() +
  geom_sf(data = dos_zebras_move) +
  geom_sf(data = mt_track_lines(dos_zebras_move), aes(color = identifier)) +
  theme_light()
```

Ahora puedes anotar algunos parametros de movimiento como:

1)  Acimut o dirección angular: como sugiere su nombre estima la dirección engrados angulares en los cuales el animal se mueve.
1)  Desplazamiento cuadrado neto (NSD): mide el movimiento del animal como una distancia cuadrada neta y permite predecir si el animal mantiene un movimiento estacionario, migatorio. Mira la @fig-NSD.
1)  Velocidad: ten en cuenta que este paquete mide velocidad lineal entre un punto A y B, en un tiempo t0 y t1.

![Comportamiento teórico basado en la estimación del NSD. Tomado de @pretorius_movement_2020](docs/Figures/NSD1.png){#fig-NSD fig-alt=""}

```{r}
dos_zebras_annot <- dos_zebras_move %>%
                      mutate(azimuth = mt_azimuth(dos_zebras_move),
                      speed = mt_speed(dos_zebras_move),
                      turnangle = mt_turnangle(dos_zebras_move),
                      NSD = (st_distance(dos_zebras_move,
                      dos_zebras_move[1,]))^2)
```

Debes tener en cuenta que es necesario agregar un límite en el cual tu infieras que el animal está en movimiento o no. Esto dependerá de tu conocimiento sobre el comportamiento de tu organismo de estudio, el terreno, el tipo de collar que utilizaste para tomar los puntos, entre otros factores. Así puedes establecer, por ahora, un limite de 5 m/min como indicador de movimiento en las zebras.

```{r warnings = FALSE}
# Cuando se mueven
zebras_movimiento <- dos_zebras_annot %>% filter(speed > set_units(5, "m/min"))
```


```{r warnings = FALSE}
zebras_movimiento$NSD <- set_units(zebras_movimiento$NSD, "km^2")
head(zebras_movimiento$NSD)

ggplot(zebras_movimiento, aes(x = mt_time(zebras_movimiento), y = NSD,
  group = identifier)) +
  geom_line() +
  facet_grid(~identifier, scales = "free_x") +
  labs(x = "", y ="Desplazamiento cuadrático neto") +
  theme_bw()
```


```{r warnings = FALSE}
ggplot(zebras_movimiento, aes(x = azimuth, y = speed)) +
  geom_point(color= alpha("black", 0.3)) +
  scale_x_units(unit = "degrees", breaks = c(-2:2) * 90) +
  scale_y_units(unit = "km/h") +
  theme_linedraw()
```

:::::: {.callout-note}
## **TAREA**
Grafica la distribución de estas zebras sobre un mapa.
:::

## Fuentes