---
title: "Visualizando patrones espaciales y temporales"
output: learnr::tutorial
runtime: shiny_prerendered
---

```{r setup, include=FALSE}
library(learnr)
library(dplyr)
library(ggplot2)
library(ggmap)
library(lubridate)

delitos <- dataviz::delitos
delitos$fecha <- ymd(delitos$fecha)

set.seed("99")
muestra_de_fechas <- delitos %>% 
    sample_n(5) %>% 
    pull(fecha)

bbox <- make_bbox(delitos$longitud, delitos$latitud)
CABA <- get_stamenmap(bbox = bbox, maptype = "toner-lite", zoom = 12)

knitr::opts_chunk$set(message = FALSE, warning = FALSE, comment = "")
```


## Trabajando con fechas

Se hace cada vez más fácil registrar eventos con gran precisión temporal además de espacial y temporal; es decir, datos que incluyen fecha y hora, o coordenadas que los ubican en un sitio determinado. Para entender datasets con datos en gran volumen que poseen atributos de posición y tiempo, es útil visualizar el ritmo en el que ocurren (diario, mensual, anual, etc) y la forma en la que se distribuyen en el espacio.

Para practicar, trabajaremos con un dataset publicado por la Ciudad Autónoma de Buenos Aires, con delitos registrados durante el 2020. Los datos fueron publicados por en el _Mapa del delito_ de la Ciudad (https://mapa.seguridadciudad.gob.ar/). 


```{r echo=FALSE, fig.align="center", out.width="85%"}
knitr::include_graphics('https://live.staticflickr.com/3410/3346781444_48266ac902_b.jpg')
```
<center>
_subtipo: "Siniestro Vial"_

</center>

Las primeras filas del dataset lucen asi:


```{r echo=FALSE}
head(delitos, 50)
```

Teniendo en cuenta que el dataset es del 2020, el Primer Año de la Pandemia, nos podríamos preguntar si se ve algún descenso brusco en robos, hurtos, homicidios o siniestros viales con el inicio de la cuarentena. Para responder eso, vamos a tener que lidiar con datos de tipo "fecha". Esto podría ser bastante engorroso, pero por suerte podemos usar una herramienta que simplifica las cosas.


### La fecha como clase de variable

La fecha es un tipo de dato que puede ser expresado de muchas maneras, dependiendo de que nos interese tener en cuenta: el día de la semana al que corresponde, el mes, el año, etc. El paquete [`lubridate`](https://lubridate.tidyverse.org/) hace fácil tomar una variable que contiene fechas en cualquier formato (por ejemplo "20/07/2020") para extraer el atributo relacionado que deseemos (como su día, "lunes", o su mes, "julio").

Para empezar, convertimos el campo "fecha" al tipo de dato especializado, que se llama... fecha (_date_). Aquí tenemos que prestar atención al formato en que aparecen los datos de la columna, en general algo como "2018-07-21" (mes, día y año) o "2018-07-21 12:14:24" (mes, día, año y hora, minutos, segundos). Con nuestros datos se da el primer caso, por lo cual la función para convertir ese campo en fecha es `ymd()`; para el segundo caso, seria `ymd_hms()`

```{r eval=FALSE}
library(lubridate)

delitos$fecha <- ymd(delitos$fecha)
```


Repasemos algunas de los nuevos trucos que podemos hacer con el tiempo. Tomemos cinco fechas elegidas al azar:

```{r}
muestra_de_fechas
```

Mediante las funciones disponibles en `lubridate`, podemos extraer:

- El día de la semana al que corresponde cada fecha:

```{r fechas1, exercise=TRUE}
wday(muestra_de_fechas)
```

```{r fechas2, exercise=TRUE}
wday(muestra_de_fechas, label = TRUE)
```
- El mes:

```{r fechas3, exercise=TRUE}
month(muestra_de_fechas)
```


```{r fechas4, exercise=TRUE}
month(muestra_de_fechas, label = TRUE)
```

- El año:

```{r fechas5, exercise=TRUE}
year(muestra_de_fechas)
```

Y varias opciones más, que se pueden repasar en https://cran.r-project.org/web/packages/lubridate/vignettes/lubridate.html

Con lo visto hasta aquí, tenemos suficiente para mostrar patrones temporales en los datos.

Empecemos por un gráfico de barras (`geom_bar()`)con la cantidad de eventos registrados por mes. Para que aparezca un conteo por mes del año, asignaremos al eje de las $x$ el _mes_ al que corresponde cada valor de la columna "fecha", o sea `month(fecha, label = TRUE)`:


```{r barras_mes, exercise=TRUE}
ggplot(radios) +
  geom___(___)
```

```{r barras_mes-hint}
ggplot(delitos) + 
    geom_bar(aes(x = month(___)))
```

```{r barras_mes-solution}
ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE)))

```

(Usamos el parámetro `label = TRUE` para obtener el nombre del mes en lugar de su número -"dic" en lugar de 12. Prueben realizar el gráfico sin ese parámetro en la llamada a `month()` para ver que pasa.)

En el gráfico, se ve una reducción drástica a partir de la cuarentena decretada en abril. Y también un incremento gradual y sostenido a partir de allí, que de todos modos no llega a los niveles pre-pandemia de enero y febrero. 

Para ver la composición interna de los conteos mensuales, cuantos casos corresponden a cada categoría, podemos hacer un gráfico de "barras apiladas" como vimos en la clase 3. La sintaxis es igual que antes, pero esta vez asignamos la variable _tipo_"_ al color de relleno de las barras determinado por el parámetro `fill`. 

```{r barras_mes2, exercise=TRUE}
ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE), ___))
```

```{r barras_mes2-hint}
ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE), fill = ___))
```

```{r barras_mes2-solution}
ggplot(delitos) + 
    geom_bar(aes(x = month(fecha, label = TRUE), fill = tipo))
```

Las barras apiladas son prolijas, pero pueden hacer difícil evaluar la evolución de categorías individuales. Recordemos que para mostrar los subconjuntos en barras independientes, una al lado de la otra, podemos usar el parámetro `position = "dodge"`. Inténtenlo agregando el parámetro al ejercicio anterior para ver como queda.

Ahora comparemos la cantidad de eventos registrados, por tipo, para cada día de la semana. Basta con usar la función que extrae el día de la semana, `wday()`, en lugar de `month()`. Lo demás es idéntico, incluyendo el uso de `label = TRUE` para que obtener el nombre del día -"lun"-, en lugar de su posición en la semana -"2", porque para `lubridate` las semanas empiezan el domingo-. 

```{r barras_dia, exercise=TRUE}
ggplot(delitos) + 
    geom_bar(aes(x = ___(fecha, label = TRUE), fill = tipo))
```

```{r barras_dia-solution}
ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo))
```

Como era de esperar, durante los fines de semana se observa una menor cantidad de eventos, aunque quizás no en la categoría Homicidio, que es difícil de discernir por su relativa escasez. Para solucionar el problema podríamos filtrar los datos antes de visualizarlos, como hemos hecho antes con la función `filter()`, mostrando sólo la categoría de interés. 

Otra opción, que probaremos ahora, es una visualización en "facetas". Se trata de generar una visualización que muestra distintos aspectos de los datos en paneles separados. Con un ejemplo vamos a dejarlo mas claro. Para realizar un gráfico en facetas, basta con sumar una línea con la función `facet_wrap(vars(x, y, ...))`, dónde "x", "y", etc son los nombres de las variables cuyas categorías recibirán paneles distintos. Intentémoslo con el último gráfico, agregando "tipo" como variable a facetar:

```{r barras_dia_facet, exercise=TRUE}
ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(___))
```

```{r barras_dia_facet-solution}
ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(tipo))
```
Obtuvimos barras separadas, pero los homicidios siguen difíciles de distinguir. Se debe a que por defecto `facet_wrap()` mantiene a escala todos los paneles, de manera que se puedan comparar cantidades de forma directa. La desventaja es que se pierde legibilidad de categorías con valores ínfimos en comparación con otras. Por eso disponemos del parámetro `scales`, que permite graficar los datos con escala libre, vía `scales = "free"`. Probemos:


```{r barras_dia_facet2, exercise=TRUE}
ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(___), scales = ___)
```

```{r barras_dia_facet2-solution}
ggplot(delitos) + 
    geom_bar(aes(x = wday(fecha, label = TRUE), fill = tipo)) +
    facet_wrap(vars(tipo), scales = "free")
```
Ahora queda más claro que los homicidios siguen un patrón diario distinto al de las demás categorías. También que, en términos relativos, las lesiones son las que mas se reducen durante los fines de semana.

También podemos evaluar el ritmo según la hora del día. ¿Cómo se haría con nuestro dataset? 

## Mirando al espacio

Pasemos ahora al análisis espacial de nuestros datos. Para facilitar la visualización volveremos a usar el paquete `ggmap`, que aporta funciones que facilitan la creación de mapas.

```{r}
library(ggmap)
```

### Obteniendo un mapa base

Como vimos en la clase previa, para obtener un mapa de fondo o "mapa base" necesitamos obtener una _bounding box_ de nuestros datos, que luego pasamos a `get_stamenmap()`.

La diferencia con nuestra experiencia anterior es que antes trabajamos con dataframes espaciales, que contenían geometrías georreferencias: polígonos, líneas, puntos. Pero este dataframe no tiene nada de eso. La buena noticia es que igual podemos hacer mapas, porque contiene dos columnas clave para ello: las de "longitud" y "latitud", que permiten ubicar puntos sobre la faz de la tierra. Y de puntos se trata, porque en general sólo podremos trabajar con columnas de coordenadas cuando los datos espaciales representen puntos, pero no líneas y polígonos. Estas geoemtrías de mayor complejidad se manejan de forma cómodo usando dataframes espaciales como los que exploramos en la clase anterior.

Así como contamos con `st_bbox()` para obtener la _bounding box_ de dataframes espaciales, cuando tenemos las coordenadas en un par de columnas recurrimos a `make_bbox()`,:

```{r}
bbox <- make_bbox(delitos$longitud, delitos$latitud)

bbox
```

En base a la "bounding box" solicitamos nuestro mapa base:

```{r}
CABA <- get_stamenmap(bbox = bbox, maptype = "toner-lite", zoom = 12)
```

Para verlo:

```{r}
ggmap(CABA)
```

### De coordenadas al mapa

De aquí en más podemos suporponer nuestros datos en distintas capas, con la misma sintaxis que conocemos de ggplot. Para mapear las ubicaciones de los delitos en el dataset, usamos `geom_point()` y los campos de longitud y latitud para los ejes $x$ e $y$:

```{r mapa1, exercise = TRUE}
ggmap(CABA) +
    geom____(data = delitos, aes(___))
```

```{r mapa1-hint}
ggmap(CABA) +
    geom_point(data = delitos, aes(x = ___, y = ___))
```


```{r mapa1-solution}
ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud))
```


Aquí nos topamos con un problema habitual al trabajar con grandes volúmenes de datos. Hay tantos puntos proyectados sobre el mapa, que se hace imposible interpretar dónde existen más o menos. Hacemos algunos ajustes: 
- un color que resalte más contra el mapa base, y que no se confunda con él
- un tamaño de punto más pequeño
-y aplicación de una ligera transparencia

Todo ello vía los atributos "color", "size" y "alpha". ¿Cuál es el valor ideal para cada uno? En general, no queda otra que recurrir a la prueba y error para encontrar la receta justa. Probemos con `color = "orange"`, `size = 0.1` y `alpha = 0.1`:

```{r mapa2, exercise = TRUE}
ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud),
               ___)
```

```{r mapa2-hint}
ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud),
               color = ___, size = ___, alpha = ___)
```

```{r mapa2-solution}
ggmap(CABA) +
    geom_point(data = delitos, aes(x = longitud, y = latitud),
               color = "orange", size = 0.1, alpha = 0.1)
```

Ahora si aparecen ciertos patrones, por ejemplo la mayor frecuencia de casos de casos cerca de las principales de circulación de la ciudad. Aún así, se hace difícil identificar de un golpe de vista las "zonas calientes", los puntos de máxima concentración. 

### Mapas de densidad

Una solución práctica para el problema de la cantidad de puntos es una técnica llamada "binning": dividir el espacio en una grilla de celdas, contar cuantos puntos caen dentro de cada una, y visualizar las cantidades agregadas. En el mundo `ggplot` esto se lleva a cabo con `geom_bind2d()`.

```{r binning, exercise = TRUE}
ggmap(CABA) +
    ___(data = delitos, aes(x = longitud, y = latitud))
```

```{r binning-solution}
ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud))
```

Ahora si, resaltan las áreas de mayor concentración de incidentes. Se puede mejorar un poco el gráfico usando:

- una mayor cantidad de celdas para aumentar la resolución
- una escala de colores diseñada para ayudar a detectar diferencias por tonalidad, como Viridis.

la cantidad de celdas se define con el parámetro "bins", por ejemplo `bins = 100`. La escala de color Viridis, como ya habíamos visto, se agrega sumando una llamada a `scale_fill_viridis_c()` -porque aquí la data es continua, si fuera discreta usaríamos `scale_fill_viridis_d()`.

```{r binning2, exercise = TRUE}
ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud), ___) +
    ___
```

```{r binning2-hint}
ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud), bins = ___) +
    scale_fill_viridis_c()
```

```{r binning2-solution}
ggmap(CABA) +
    geom_bin2d(data = delitos, aes(x = longitud, y = latitud), bins = 100) +
    scale_fill_viridis_c()
```


Una alternativa al _binning_ es la llamada _kernel density estimation_, muy utilizada en aplicaciones de análisis espacial para estimar la intensidad de una determinada variable en cualquier punto del área analizada, incluso en aquellos para los cuales no hay observaciones. La idea es asumir que los valores observados corresponden a una distribución continua sobre el espacio, y determinar cual es la más probable en base a los puntos donde existen datos. Podemos visualizar esta distribución estimada `geom_density2d_filled` así:

```{r density, exercise = TRUE}
ggmap(CABA) +
    ___(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5) 
```

```{r density-solution}
ggmap(CABA) +
    geom_density2d_filled(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5) 
```
Nótese que aplicamos transparencia usando el parámetro `alpha = 0.5`. ¿Por qué? ¿Qué pasa si lo quitamos?


### Visualizando multiples categorías

Hasta aquí hemos analizado la distribución espacial de eventos en su totalidad, sin diferenciar su tipo. Veamos ahora las diferencias por categoría. Podemos reintentar el mapa de puntos, esta vez diferenciándolos por color. Recuperamos el código que usamos antes para mostrar puntos, y esta vez asignamos la columna "tipo" al atributo estético `color`:

```{r mapa_categorias, exercise = TRUE}
ggmap(CABA) +
    geom_point(data = delitos, 
               aes(x = longitud, y = latitud, ___),
               size = 0.1, alpha = 0.1)
```

```{r mapa_categorias-hint}
ggmap(CABA) +
    geom_point(data = delitos, 
               aes(x = longitud, y = latitud, color = ___),
               size = 0.1, alpha = 0.1)
```

```{r mapa_categorias-solution}
ggmap(CABA) +
    geom_point(data = delitos, 
               aes(x = longitud, y = latitud, color = tipo),
               size = 0.1, alpha = 0.1)
```


Aquí tenemos dos problemas: 

* La leyenda ("tipo_delito") es difícil de leer, dado que muestra los puntos tal como los definimos: pequeños y con mucha transparencia. Esos atributos son útiles en el mapa, donde tenemos cientos de miles de puntos, pero muy poco prácticos para la leyenda, donde sólo hay un minúsculo punto por etiqueta.

* Los puntos sobre el mapa se superponen en tal medida que es difícil identificar patrones espaciales distintos según su categoría.

El primer problema se resuelve fijando "a mano" los atributos de la leyenda, asi:

```{r}
ggmap(CABA) +
    geom_point(data = delitos,
               aes(x = longitud, y = latitud, color = tipo),
               size = 0.1, alpha = 0.1) +
    guides(color = guide_legend(override.aes = list(size = 1, alpha = 1)))
```

El segundo, usando facetado para mostrar en su propio mapa a cada categoría:

```{r}
ggmap(CABA) +
    geom_point(data = delitos,
               aes(x = longitud, y = latitud, color = tipo),
               size = 0.1, alpha = 0.1) +
    guides(color = guide_legend(override.aes = list(size = 1, alpha = 1))) +
    facet_wrap(vars(tipo))
```

El facetado ayuda a que no se nos mezclen los colores, y hace evidente cuales categorías son mas frecuentes que otras. Pero  con nuestra abundancia de puntos no ayuda encontrar los sitios de alta concentración, y hace que se pierdan de vista los casos de la categoría poco frecuente (homicidios).

Para hacer las diferencias aún mas nítidas, podemos facetar una estimación de densidad en lugar de puntos. ¿Cómo lo haríamos?

```{r facetado_densidad, exercise = TRUE}
ggmap(CABA) +
    ___(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5) +
    facet_wrap(vars(tipo))
```

```{r facetado_densidad-solution}
ggmap(CABA) +
    geom_density2d_filled(data = delitos, aes(x = longitud, y = latitud), alpha = 0.5) +
    facet_wrap(vars(tipo))
```

## Combinando espacio y tiempo


El facetado también nos permite visualizar el cambio de posición a través del tiempo.

Por ejemplo, podemos comparar cierto tipo delito (hurto sin violencia) mostrando dónde ocurre en cada día de la semana.

Primero activamos el paquete `dplyr` para acceder a su función de filtrado de datos,

```{r}
library(dplyr)
```

Y luego mostramos sólo las filas del dataframe donde la columna tipo contiene "Homicidio", 
 - en forma de puntos en el mapa (`geom_point()`), 
 - con el "subtipo" de homicidio representado por el `color` de los puntos
 - y un facetado por día de la semana (`facet_wrap(vars(wday(fecha, label = TRUE)))`)


```{r homicidios_por_dia, exercise = TRUE}
ggmap(CABA) +
    ___(data = filter(delitos, ___),
               aes(x = longitud, y = latitud, ___), alpha = .5) +
    ___(vars(wday(fecha, label = TRUE)))

```

```{r homicidios_por_dia-hint}
ggmap(CABA) +
    geom_point(data = filter(delitos, tipo == ___),
               aes(x = longitud, y = latitud, color = ___), alpha = .5) +
    facet_wrap(vars(wday(fecha, label = TRUE)))

```

```{r homicidios_por_dia-solution}
ggmap(CABA) +
    geom_point(data = filter(delitos, tipo == "Homicidio"),
               aes(x = longitud, y = latitud, color = subtipo), alpha = .5) +
    facet_wrap(vars(wday(fecha, label = TRUE)))

```
Vale aclarar que el poco elegante `facet_wrap(vars(wday(fecha, label = TRUE)))` podría cambiarse por un más legible `facet_wrap(vars(dia_semana))` si como paso previo agregamos al dataframe la columna "dia_semana", guardando allí el valor obtenido con `wday()`. 

Volviendo al tiempo y el espacio, también podemos concentrarnos en un tipo de delito en particular, y evaluar en que zonas se concentra de acuerdo a la hora del día. Nuestra data ya tiene la hora del día declarada en una columna, "franja". Si no fuera así, y tuviéramos la hora como parte de la fecha (estilo "2020-09-18 14:00:00") podríamos obtenerla con ayuda de `hour()` que funciona de forma similar a las ya vistas `month()` y `wday()`.

Entonces, mostremos sólo las filas del dataframe donde la columna tipo contiene "Hurto (sin violencia)", 
 - en forma de mapa de densidad (`geom_density2d_filled()`), 
 - y un facetado por hora del día (`facet_wrap(vars(franja))`)


```{r hurtos_por_hora, exercise = TRUE}
ggmap(CABA) +
    ___(data = filter(delitos, ___),
               aes(x = longitud, y = latitud), alpha = .5) +
    ___
```

```{r hurtos_por_hora-hint}
ggmap(CABA) +
    geom_density2d_filled(data = filter(delitos, tipo == ___),
               aes(x = longitud, y = latitud), alpha = .5) +
    facet_wrap(___)
```

```{r hurtos_por_hora-solution}
ggmap(CABA) +
    geom_density2d_filled(data = filter(delitos, tipo == "Hurto (sin violencia)"),
               aes(x = longitud, y = latitud), alpha = .5) +
    facet_wrap(vars(franja))
```
En el resultado se puede ver como los hurtos se concentran nítidamente en las áreas de mayor actividad comercial durante el día (lo que los porteños llaman "el centro"), sobre todo desde el mediodía hasta 4 o 5 de la tarde, cuando pierde intensidad y se dispersa en la dirección de las principales avenidas de la ciudad. 

Para terminar, pulimos la visualización 
 - filtrando las filas que registran la franja horaria como desconocida (`!is.na(franja)`)
 - retirando la leyenda, ya que nos interesa mostrar como se mueve la densidad a lo largo del día más que las cantidades
 - eligiendo la cantidad de filas en la que se distribuirán las facetas  (`nrow = 4`)
 - agregando título, subtítulo, y nota al pie con fuente
 - eligiendo un tema apropiado

```{r}
ggmap(CABA) +
    geom_density2d_filled(data = filter(delitos, !is.na(franja), tipo == "Hurto (sin violencia)"),
               aes(x = longitud, y = latitud), alpha = .5) +
    guides(fill = FALSE) +
    facet_wrap(vars(franja), nrow = 4) +
    labs(title = "Ciudad de Buenos Aires: concentración espacial de hurtos",
         subtitle = "según hora del día, durante el año 2020",
         caption = "fuente: https://mapa.seguridadciudad.gob.ar") +
    theme_void()
```