Mapas Interactivos con Leaflet: Ecuador y Siniestros de Tránsito
Autor: Nelson Abad (nelson_abadz@outlook.com)
Fecha: 2023-02-17
Introducción
Desde el año 2017 hasta el año 2022, Ecuador tuvo aproximadamente ciento cuarenta mil accidentes de tránsito, en los cuales once mil veintisiete hubo al menos un fallecido. Según Our World in Data en el año 2019 la séptima causa de muerte más frecuente en Ecuador era por accidentes de tránsito o lesiones provocadas por los mismos, llegando a ser cuatro mil setecientos ochenta y ocho, y para el año 2022 la causa de muerte número uno en Ecuador son los accidentes de tránsito, siendo el único país junto con Paraguay en tener esta causa como la principal según Latinometrics.
Conociendo entonces cuál es la situación de los siniestros de tránsito en Ecuador, lo que se busca es aprovechar estos datos para así poder explicar como se construyen ciertos gráficos que nos ayudan a entender geográficamente en que partes del país hay mayor número de accidentes de tránsito y qué información valiosa podemos rescatar de cada uno de estos accidentes.
Librerías a Utilizar
Para crear los mapas se utilizará la librería \(\texttt{leaflet}\), la librería \(\texttt{geojsonio}\) para poder manipular e importar datos tipo Geo JSON, la librería \(\texttt{knitr}\) para crear tablas y la librería \(\texttt{tidyverse}\) para poder realizar análisis exploratorio y manipulación de la base de datos y poder utilizar el operador \(\texttt{pipe}\) para concatenar las funciones.
library(tidyverse)
library(leaflet)
library(geojsonio)
library(knitr)Datos a Analizar
Los datos a analizar provienen de la Agencia Nacional de Tránsito del Ecuador, el cual es un registro que va (hasta la fecha de elaborado el presente trabajo) desde enero del 2017 hasta diciembre del 2022.
siniestros = read_csv("BDD_DICIEMBRE.csv")Seleccionamos ahora las variables de interés.
siniestros = siniestros %>%
select(ANIO,
LESIONADOS,
FALLECIDOS,
LATITUD_Y,
LONGITUD_X,
DPA_1,
DPA_2,
PROVINCIA,
CANTON,
ZONA_PLANIFICACION,
ZONA,
FECHA,
HORA,
FERIADO,
TIPO_DE_SINIESTRO,
TIPO_DE_VEHICULO_1,
)Análisis Exploratorio
Observemos primero cuantos accidentes se produjeron por año.
siniestros %>%
group_by(ANIO) %>%
summarise(accidentes = n()) %>%
kable(col.names = c("Año", "Numero de Accidentes"),
align = rep("c", 6))| Año | Numero de Accidentes |
|---|---|
| 2017 | 28967 |
| 2018 | 25530 |
| 2019 | 24595 |
| 2020 | 16972 |
| 2021 | 21352 |
| 2022 | 21739 |
El año en donde más accidentes se produjeron fue el 2017. Tomaremos este año de referencia para nuestros mapas.
siniestros = siniestros %>% filter(ANIO == 2017)Observemos ahora en un gráfico cuáles fueron las provincias con mayor número de accidentes.
siniestros %>%
mutate(PROVINCIA = ifelse(PROVINCIA == "SANTO DOMINGO DE LOS TSACHILAS",
"SANTO DOMINGO",
PROVINCIA)) %>%
mutate(PROVINCIA = PROVINCIA %>% factor() %>% fct_infreq() %>% fct_rev()) %>%
ggplot(mapping = aes(x = PROVINCIA, y = (..count..)/sum(..count..))) +
geom_bar(color = "black", fill = "#00005a") +
xlab("") + ylab("") +
coord_flip() +
theme_minimal()
Se puede observar que las provincias con mayor proporción de accidentes de tránsito son Pichincha y Guayas, la provincia con menor proporción de accidentes de tránsito fue Galápagos.
Ahora, una variable que podría ser interesante de analizar sería la zona, esto indica si el accidente se produjo en zona urbana o rural, en el siguiente gráfico podemos observar que tanta es la proporción de accidentes que se dieron en zona urbana o rural y por provincia.
siniestros = siniestros %>%
mutate(ZONA = ifelse(ZONA == "Rural", "RURAL", ZONA))
siniestros %>%
mutate(PROVINCIA = ifelse(PROVINCIA == "SANTO DOMINGO DE LOS TSACHILAS",
"SANTO DOMINGO",
PROVINCIA)) %>%
mutate(PROVINCIA = PROVINCIA %>% factor() %>% fct_infreq() %>% fct_rev()) %>%
ggplot(mapping = aes(x = PROVINCIA,
y = (..count..)/sum(..count..),
fill = ZONA)) +
geom_bar(color = "black", position = "dodge") +
xlab("") + ylab("") +
scale_fill_manual(values = c("#3f75f2", "#00005a")) +
guides(fill = guide_legend(title = "Zona:")) +
coord_flip() +
theme_minimal()
Se puede notar que en casi todas las provincias hay mayor número de accidentes en zonas urbanas que en zonas rurales. En la siguiente tabla podemos ver en términos absolutos y de todo el país el número de accidentes por zona.
siniestros %>%
group_by(ZONA) %>%
summarise(accidentes = n()) %>%
kable(col.names = c("Zona", "Numero de Accidentes"),
align = rep("c", 2))| Zona | Numero de Accidentes |
|---|---|
| RURAL | 8948 |
| URBANA | 20019 |
Para los mapas creados a continuación no tendremos en cuenta las islas Galápagos, debido a que, por su condición geográfica, la cantidad de accidentes que se dan es muy escasa, además que podría dificultar en análisis del mapa debido a su lejanía con las costas del Ecuador continental
Mapa de Puntos
Para poder realizar este mapa necesitamos las variables que nos den las coordenadas geográficas para poder graficar cada punto, las cuales son latitud y longitud. Lo primero que debemos hacer es crear un mapa y fijar la vista en Ecuador. Para no saturar el mapa solo graficaremos los accidentes en los cuales solo haya habido al menos un fallecido o un lesionado, además que podemos modificar el diseño original.
siniestros %>%
filter(FALLECIDOS > 0 | LESIONADOS > 0) %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.DE)Una vez creado el mapa y los datos están listos para ser usados lo que tenemos que hacer es graficar los puntos, cada uno de estos puntos va a representar un accidente y serán ubicados geográficamente utilizando las variables latitud y longitud.
siniestros %>%
filter(FALLECIDOS > 0 | LESIONADOS > 0) %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addCircles(lng = ~LONGITUD_X,
lat = ~LATITUD_Y)Prácticamente, el mapa ya está creado, lo que faltaría por hacer es agregar ciertos detalles que puedan hacerlo más comunicativo e interactivo. Por ejemplo, hacer los puntos más pequeños para que el mapa no se vea tan saturado, cambiar el color de los puntos y agregar información relevante del accidente que se va a desplegar si damos clic en cualquiera de los puntos.
siniestros %>%
filter(FALLECIDOS > 0 | LESIONADOS > 0) %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addCircles(lng = ~LONGITUD_X,
lat = ~LATITUD_Y,
radius = 0.05,
color = "#3f75f2",
popup = ~paste("Provincia:", PROVINCIA, "<br/>",
"Cantón:", CANTON, "<br/>",
"Lesionados:", LESIONADOS, "<br/>",
"Fallecidos:", FALLECIDOS, "<br/>",
"Zona:", ZONA, "<br/>",
"Siniestro:", TIPO_DE_SINIESTRO, "<br/>",
"Vehículo:", TIPO_DE_VEHICULO_1))Se pueden realizar acercamientos en el mapa para tener mayor detalle en que provincia, cantón, parroquia o dirección se dio el accidente.
Para culminar la construcción de este mapa, podemos agregar un detalle extra, vamos a definir de un color los accidentes que se hayan dado en zona urbana y de otro color los accidentes que se hayan dado en zona rural, para esto también debemos agregar una leyenda.
siniestros %>%
filter(FALLECIDOS > 0 | LESIONADOS > 0) %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addCircles(lng = ~LONGITUD_X,
lat = ~LATITUD_Y,
color = ~ifelse(ZONA == "RURAL", "#3f75f2", "#00005a"),
radius = 0.05,
popup = ~paste("Provincia:", PROVINCIA, "<br/>",
"Cantón:", CANTON, "<br/>",
"Lesionados:", LESIONADOS, "<br/>",
"Fallecidos:", FALLECIDOS, "<br/>",
"Zona:", ZONA, "<br/>",
"Siniestro:", TIPO_DE_SINIESTRO, "<br/>",
"Vehículo:", TIPO_DE_VEHICULO_1)) %>%
addLegend(title = "Zona", labels = c("Rural", "Urbana"),
colors = c("#3f75f2", "#00005a"), opacity = 1.0)Aparentemente, son más los accidentes que se dan en zonas rurales que en zonas urbanas; sin embargo, sabemos por un análisis previo que el número de accidentes que se dio en zonas rurales fue ocho mil novecientos cuarenta y ocho; en cambio, en las zonas urbanas se dieron veinte mil diecinueve. Lo que si podemos observar es que los accidentes en las zonas rurales son más dispersos que los accidentes en las zonas urbanas, esto debido a que en las zonas urbanas hay una mayor densidad poblacional.
Si quisiéramos observar únicamente los accidentes que se dieron en las provincias con mayor número de accidentes (en este caso Pichincha y Guayas) simplemente agregamos un par de sentencias adicionales a la función \(\texttt{filter()}\).
siniestros %>%
filter((PROVINCIA == "PICHINCHA" | PROVINCIA == "GUAYAS") &
(FALLECIDOS > 0 | LESIONADOS > 0)) %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addCircles(lng = ~LONGITUD_X,
lat = ~LATITUD_Y,
color = ~ifelse(ZONA == "RURAL", "#3f75f2", "#00005a"),
radius = 0.05,
popup = ~paste("Provincia:", PROVINCIA, "<br/>",
"Cantón:", CANTON, "<br/>",
"Lesionados:", LESIONADOS, "<br/>",
"Fallecidos:", FALLECIDOS, "<br/>",
"Zona:", ZONA, "<br/>",
"Siniestro:", TIPO_DE_SINIESTRO, "<br/>",
"Vehículo:", TIPO_DE_VEHICULO_1)) %>%
addLegend(title = "Zona", labels = c("Rural", "Urbana"),
colors = c("#3f75f2", "#00005a"), opacity = 1.0)Mapa Coroplético
Un mapa coroplético se construye sombreando áreas con mayor o menor intensidad de una misma gama cromática, en donde la intensidad de esa gama cromática representa los distintos valores de una variable estadística característica de esa región geográfica. En este caso particular, las áreas de interés serían las provincias y la variable estadística sería el número de accidentes.
Para realizar este mapa necesitamos delimitar cada provincia geográficamente, para esto necesitamos importar datos Geo JSON de las provincias del Ecuador y manipular estos datos Geo JSON para poder asignar a cada provincia la información deseada. Para esto utilizamos la librería \(\texttt{geojsonio}\).
provincias = geojson_read(
"https://raw.githubusercontent.com/Handerica/COVID-19_Ecuador/master/provs-ec.json",
what = "sp")Podemos ver que el objeto \(\texttt{provincias}\) es un dataframe del tipo espacial-poligonal.
class(provincias)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"Podemos comprobar si los datos delimitan correctamente las provincias del Ecuador.
provincias %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.HOT) %>%
addPolygons(
weight = 2,
opacity = 1,
color = "black")Observemos ahora el número de accidentes, fallecidos y lesionados en cada provincia. Para este caso usaremos los departamentos.
provincia_data = siniestros %>%
select(DPA_1, FALLECIDOS, LESIONADOS) %>%
group_by(DPA_1) %>%
summarise(TOTAL = n(), FALLECIDOS = sum(FALLECIDOS), LESIONADOS = sum(LESIONADOS)) %>%
rename(dpa_provin = DPA_1)
provincia_data %>%
kable(align = rep("c", 4),
col.names = c("Provincia", "Accidentes", "Fallecidos", "Lesionados"))| Provincia | Accidentes | Fallecidos | Lesionados |
|---|---|---|---|
| 1 | 1497 | 92 | 1127 |
| 2 | 167 | 23 | 245 |
| 3 | 176 | 38 | 232 |
| 4 | 186 | 17 | 149 |
| 5 | 500 | 131 | 272 |
| 6 | 790 | 88 | 295 |
| 7 | 687 | 86 | 577 |
| 8 | 250 | 74 | 245 |
| 9 | 8414 | 455 | 8073 |
| 10 | 1323 | 81 | 487 |
| 11 | 537 | 25 | 254 |
| 12 | 910 | 181 | 957 |
| 13 | 1305 | 150 | 1141 |
| 14 | 176 | 31 | 189 |
| 15 | 95 | 25 | 144 |
| 16 | 54 | 10 | 81 |
| 17 | 9363 | 375 | 5345 |
| 18 | 1120 | 92 | 760 |
| 19 | 82 | 17 | 84 |
| 20 | 2 | 0 | 1 |
| 21 | 76 | 35 | 65 |
| 22 | 34 | 15 | 25 |
| 23 | 582 | 72 | 606 |
| 24 | 641 | 40 | 664 |
Ahora, como muchas veces los datos GEO JSON tienen ciertas imprecisiones, hay zonas dentro del Ecuador que las define como ZONA NO DELIMITADA.
Como estos datos son relacionales en donde la llave es \(\texttt{dpa_provin}\), por lo tanto, podemos unir los datos mediante la función \(\texttt{full_join()}\). Recordar que el objeto \(\texttt{provincias}\) es un dataframe del tipo espacial-poligonal, para hacer deferencias a los datos internos usamos \(\texttt{@data}\).
provincias@data = provincias@data %>%
select(dpa_provin, dpa_despro) %>%
mutate(dpa_provin = as.numeric(dpa_provin)) %>%
full_join(provincia_data, by = "dpa_provin")Tenemos listos los datos y las delimitaciones de las provincias, solo falta definir la gama cromática.
pal <- colorNumeric("Blues", NULL, na.color = "black")Creamos entonces el mapa coroplético.
provincias %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addTiles() %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addPolygons(fillColor = ~pal(TOTAL),
weight = 2,
opacity = 1,
color = "black",
dashArray = "3",
fillOpacity = 0.85,
popup = ~paste("Provincia:", dpa_despro, "<br/>",
"Total:", TOTAL, "<br/>",
"Fallecidos:", FALLECIDOS, "<br/>",
"Lesionados:", LESIONADOS)) %>%
addLegend(title = "Número de <br/> Accidentes", pal = pal, values = ~TOTAL, opacity = 1.0)Dando clic en cualquiera de las provincias se va a desplegar el número de accidentes, de fallecidos y lesionados en el año 2017. Las zonas no delimitadas fueron pintadas de color negro y al no tener información como tal tienen únicamente \(\texttt{NA}\).
Si quisiéramos hacer un mapa coroplético de los cantones del Ecuador, tendríamos que hacer un proceso muy similar, evidentemente podría ser un poco más desafiante debido a que son doscientos veintiún cantones; sin embargo, podemos aprovechar la ventaja de que son datos relacionales y su manipulación es sencilla.
Importamos los datos GEO JSON para delimitar los cantones.
cantones = geojson_read(
"https://data.humdata.org/dataset/e66dbc70-17fe-4230-b9d6-855d192fc05c/resource/6fa37b41-ad28-40a6-9641-3b4efd4dbe13/download/ecuador.geojson",
what = "sp")Comprobamos si los datos delimitan correctamente los cantones del Ecuador.
cantones %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addProviderTiles(providers$OpenStreetMap.HOT) %>%
addPolygons(
weight = 2,
opacity = 1,
color = "black")Ahora, en el objeto \(\texttt{siniestros}\) puede que no estén algunos cantones, esto debido a que quizás no se han registrado accidentes de tránsito en ese sector; sin embargo, en el objeto \(\texttt{cantones}\) si están todos los cantones con sus respectivas delimitaciones, aquí es muchísimo más útil el uso de sentencias \(\texttt{SQL}\) como lo son \(\texttt{full_join()}\) debido a que aunque exista un cantón que no haya tenido un solo accidente de tránsito este seguirá apareciendo pero con datos \(\texttt{NA}\). Esto es especialmente útil para que el código mantenga su correcto funcionamiento cuando se están realizando manipulación de datos.
cantones_data = siniestros %>%
select(DPA_2, FALLECIDOS, LESIONADOS) %>%
group_by(DPA_2) %>%
summarise(TOTAL = n(), FALLECIDOS = sum(FALLECIDOS), LESIONADOS = sum(LESIONADOS)) %>%
rename(DPA_CANTON = DPA_2)
cantones@data = cantones@data %>%
select(DPA_CANTON, DPA_DESCAN, DPA_DESPRO) %>%
mutate(DPA_CANTON = as.numeric(DPA_CANTON)) %>%
full_join(cantones_data, by = "DPA_CANTON")Podríamos definir una gama cromática distinta.
pal <- colorNumeric("Reds", NULL, na.color = "black")Finalmente, creamos el mapa coroplético.
cantones %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addTiles() %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addPolygons(fillColor = ~pal(TOTAL),
weight = 2,
opacity = 1,
color = "black",
dashArray = "3",
fillOpacity = 0.85,
popup = ~paste("Provincia:", DPA_DESPRO, "<br/>",
"Canton:", DPA_DESCAN, "<br/>",
"Total:", TOTAL, "<br/>",
"Fallecidos:", FALLECIDOS, "<br/>",
"Lesionados:", LESIONADOS)) %>%
addLegend(title = "Total:", pal = pal, values = ~TOTAL, opacity = 1.0)Realizamos un último filtrado, en este caso para observar las provincias con mayor número de accidentes.
cantones_data = siniestros %>%
select(DPA_2, FALLECIDOS, LESIONADOS, PROVINCIA) %>%
filter(PROVINCIA == "GUAYAS" | PROVINCIA == "PICHINCHA") %>%
group_by(DPA_2) %>%
summarise(TOTAL = n(), FALLECIDOS = sum(FALLECIDOS), LESIONADOS = sum(LESIONADOS)) %>%
rename(DPA_CANTON = DPA_2)
cantones@data = cantones@data %>%
select(DPA_CANTON, DPA_DESCAN, DPA_DESPRO) %>%
mutate(DPA_CANTON = as.numeric(DPA_CANTON)) %>%
full_join(cantones_data, by = "DPA_CANTON")
cantones %>%
leaflet() %>%
setView(-78, -1.8, 6) %>%
addTiles() %>%
addProviderTiles(providers$OpenStreetMap.DE) %>%
addPolygons(fillColor = ~pal(TOTAL),
weight = 2,
opacity = 1,
color = "black",
dashArray = "3",
fillOpacity = 0.85,
popup = ~paste("Provincia:", DPA_DESPRO, "<br/>",
"Canton:", DPA_DESCAN, "<br/>",
"Total:", TOTAL, "<br/>",
"Fallecidos:", FALLECIDOS, "<br/>",
"Lesionados:", LESIONADOS)) %>%
addLegend(title = "Total:", pal = pal, values = ~TOTAL, opacity = 1.0)