#Configuración y librerías:

#Carga de datos

CDMX_incidentes_viales <- read.csv("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/incidentes-viales-ssc.csv",
                            encoding = "UTF-8")
CDMX_alcaldias <- read_sf("/Users/martinsingla/Documents/UTDT - MEU/Ciencia de Datos para Ciudades II/INFO SELECCIONADA/CDMX_Alcaldias/alcaldias.shp")
CDMX_incidentes_viales$FECHA.EVENTO <- as.Date(CDMX_incidentes_viales$FECHA.EVENTO, format = "%Y-%m-%d")

CDMX_incidentes_viales$MES_EVENTO <- factor(CDMX_incidentes_viales$MES_EVENTO , levels= c("Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio", "Julio", "Agosto", "Septiembre", "Octubre", "Noviembre", "Diciembre"))

#Incidentes históricos: mensuales

ggplot() +
  geom_bar(data= CDMX_incidentes_viales, aes(floor_date(FECHA.EVENTO, "month"), fill = CONDICIÓN)) +
  labs(x= "Meses",
       y= "Cantidad de incidentes",
       title= "CDMX: Incidentes viales",
       subtitle= "Cantidad mensual de incidentes viales por condición de accidentado (CDMX, 2018-2019)",
       caption= "Fuente: elaboración propia a partir de datos.cdmx.gob.mx")

#¿¿Escandaloso pensar qué pasó a mediados de 2018 que ascendieron tanto la cantidad de accidentes reportados?? NO. Debido a este punto chequeé la descripción informativa de la fuente , en la cual explica que a partir de JUNIO de 2018 se ampliaron las fuentes por las que se recoletaba esta información teniendo en cuenta las unidades médicas que fueron a atender los casos. Los casos fallecidos están completos, ese dato se ha reportado siempre, sin embargo, los accidentados no estan completos, por lo que no podemos usar la serie pre-Junio 2018.

#Incidentes históricos: Diarios durante un mes (diciembre 2019)

ggplot() +
  geom_bar(data = CDMX_incidentes_viales %>%
              filter(between(FECHA.EVENTO, as.Date("2019-12-01"), as.Date("2019-12-31"))),
            aes(x= floor_date(FECHA.EVENTO, "days"), fill= TIPO.DE.EVENTO)) +
  labs(x= "Días",
       y= "Cantidad de incidentes",
       title= "CDMX: Incidentes viales",
       subtitle = "Cantidad de incidentes viales diarios por tipo (CDMX, Diciembre 2019)",
       caption = "Fuente: elaboración propia a partir de datos.cdmx.gob.mx")
## Warning: position_stack requires non-overlapping x intervals

#Nótese que hay picos viernes y sábados, así como una gran caída en la cantidad de accidentes sucedidos el 25 de diciembre. Este último punto tal vez no se deba a que la gente conduzca mejor, sino que debe conducir mucha menos gente. Sería interesante comparar el ratio de vehículos efectivamente en la calle siendo conducidos en relación a la cantidad de accidentes. Tal vez a modo de hipótesis, si bien en Navidad en términos absolutos se producen menos accidentes, sin embargo en términos proporcionales a la cantidad de vehículos que se conducen, puede que el ratio sea superior. #La proporción de accidentes según tipo de accidente respecto al total pareciera ser constante a lo largo del período

#Accidentes por hora, según tipo (2019)

CDMX_incidentes_viales_porhora <- CDMX_incidentes_viales %>%
  filter(AÑO_EVENTO == 2019) %>%
  count(TIPO.DE.EVENTO, hora_base = hour(hm(HORA_EVENTO))) %>%
  mutate(Particip.porcentual = (n/sum(n))*100)

ggplot() +
  geom_line(data= CDMX_incidentes_viales_porhora,
            aes(x= hora_base, y= n, group = TIPO.DE.EVENTO, color= TIPO.DE.EVENTO)) +
  labs(title= "CDMX: Incidentes viales",
       subtitle= "Distibución horaria de incidentes viales (acumulado 2019)",
       x= "hora",
       y= "frecuencia de incidentes",
       caption= "Fuente: Elaboración propia a partir de datos.cdmx.gob.ar")

#Interesante, hay un pico enorme de accidentes entre las 7 y las 10 am. , la hora pico de tráfico vial en CDMX.

#Ubicación espacial de los accidentes:

bbox <- make_bbox(CDMX_incidentes_viales$COORDENADA.X, CDMX_incidentes_viales$COORDENADA.Y)

CDMX <- get_stamenmap(bbox = bbox,
                      maptype = "toner-lite",
                      zoom = 9)
## Source : http://tile.stamen.com/toner-lite/9/113/226.png
## Source : http://tile.stamen.com/toner-lite/9/114/226.png
## Source : http://tile.stamen.com/toner-lite/9/115/226.png
## Source : http://tile.stamen.com/toner-lite/9/116/226.png
## Source : http://tile.stamen.com/toner-lite/9/113/227.png
## Source : http://tile.stamen.com/toner-lite/9/114/227.png
## Source : http://tile.stamen.com/toner-lite/9/115/227.png
## Source : http://tile.stamen.com/toner-lite/9/116/227.png
## Source : http://tile.stamen.com/toner-lite/9/113/228.png
## Source : http://tile.stamen.com/toner-lite/9/114/228.png
## Source : http://tile.stamen.com/toner-lite/9/115/228.png
## Source : http://tile.stamen.com/toner-lite/9/116/228.png
ggmap(CDMX) +
  geom_sf(data= CDMX_alcaldias, inherit.aes =  FALSE, fill = NA) +
  geom_point(data = CDMX_incidentes_viales %>%
               filter(AÑO_EVENTO == 2019), 
             aes(x= COORDENADA.X, y= COORDENADA.Y), size = 0.00001, color = "orange", alpha = 0.1)
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#En este mapa tenemos un problema, el dataset está contabilizando también accidentes sucedidos en Cuernavaca, Toluca y Puebla (supuestamente), que no son parte de la Ciudad de México y su área metropolitana. Debemos acotar el dataset a las coordenadas geográficas que engloben la región metropolitana del Valle de México (RMVM), para poder visualizar más en detalle la CDMX.Notamos que de todos modos, en la variable “Alcaldía” se detallan las 16 alcaldías de CDMX. Esto quiere decir que hay coordenadas mal introducidas en el dataset. Corroboramos en la fuente que efectivamente es así, en esta serie hay errores de geo-referenciación, que se incluyeron para no perder la serie completa. #Todos los casos aislados son de 2018. Al filtrar a “2019” ya no se los ve.

#Operación para quitar los puntos que caen fuera de la CDMX. Estos fueron efectivamente accidentes que sucedieron en CDMX, pero no sabemos el punto geo-localizado donde sucedieron por lo que los removeremos.

spdf <- SpatialPointsDataFrame(coords = CDMX_incidentes_viales[, c("COORDENADA.X", "COORDENADA.Y")], data = CDMX_incidentes_viales)
sppoly <- as_Spatial(CDMX_alcaldias)
proj4string(spdf) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
proj4string(sppoly) <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
## Warning in `proj4string<-`(`*tmp*`, value = new("CRS", projargs = "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")): A new CRS was assigned to an object with an existing CRS:
## +proj=longlat +datum=WGS84 +no_defs
## without reprojecting.
## For reprojection, use function spTransform
CDMX_incidentes_viales <- as.data.frame(spdf[!is.na(over(spdf, as(sppoly, "SpatialPolygons"))),])

#De 35,042 casos, el dataframe pasa a tener 34,890 casos. Es decir, eliminamos 152 casos.

#Cramos un nuevo bbox, y volvemos a descarcar el mapa base y crear el mapa de puntos de accidentes viales.

bbox <- c(left = -99.366341, bottom = 19.0461678, right = -98.937752, top = 19.594983)

CDMX <- get_stamenmap(bbox= bbox, maptype = "toner-lite", zoom= 11)
## Source : http://tile.stamen.com/toner-lite/11/458/910.png
## Source : http://tile.stamen.com/toner-lite/11/459/910.png
## Source : http://tile.stamen.com/toner-lite/11/460/910.png
## Source : http://tile.stamen.com/toner-lite/11/461/910.png
## Source : http://tile.stamen.com/toner-lite/11/458/911.png
## Source : http://tile.stamen.com/toner-lite/11/459/911.png
## Source : http://tile.stamen.com/toner-lite/11/460/911.png
## Source : http://tile.stamen.com/toner-lite/11/461/911.png
## Source : http://tile.stamen.com/toner-lite/11/458/912.png
## Source : http://tile.stamen.com/toner-lite/11/459/912.png
## Source : http://tile.stamen.com/toner-lite/11/460/912.png
## Source : http://tile.stamen.com/toner-lite/11/461/912.png
## Source : http://tile.stamen.com/toner-lite/11/458/913.png
## Source : http://tile.stamen.com/toner-lite/11/459/913.png
## Source : http://tile.stamen.com/toner-lite/11/460/913.png
## Source : http://tile.stamen.com/toner-lite/11/461/913.png
ggmap(CDMX) +
  geom_sf(data= CDMX_alcaldias, inherit.aes = FALSE, fill = NA) +
  geom_point(data= CDMX_incidentes_viales %>%
               filter(AÑO_EVENTO == "2019"), 
             aes(x= COORDENADA.X, y= COORDENADA.Y), color= "salmon", alpha= 0.5, size= 0.0000000001) +
  labs(x= "",
       y= "",
       title= "CDMX: Incidentes viales",
       subtitle= "Ubicación de indidentes viales acaecidos en 2019",
       caption= "Fuente: Elaboración propia a partir de datos.cdmx.gob.mx")
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#Los datos están muy compactados por toda la ciudad, será mejor crear un mapa de densidad.

ggmap(CDMX) +
  geom_bin2d(data = CDMX_incidentes_viales %>%
               filter(AÑO_EVENTO == "2019"),
             aes(x= COORDENADA.X, y= COORDENADA.Y), bins= 100) +
  scale_fill_viridis_b() +
  labs(x= "",
       y= "",
       title= "CDMX: incidentes viales",
       subtitle= "Mapa de densidad de accidentalidad vial (total 2019)",
       caption= "Fuente: elaboración propia a partir de datos.cdmx.gob.mx")

#Vemos muy claramente que las alcaldías centrales de la ciudad concentran la mayor cantidad de accidentes viales, lo cual es deducible ya que así mismo concentran la mayor parte de la afluencia de vehículos.

#Hacemos un facetado por tipo de accidentes

ggmap(CDMX) +
  geom_density_2d(data= CDMX_incidentes_viales %>%
               filter(AÑO_EVENTO == "2019"), 
             aes(x= COORDENADA.X, y= COORDENADA.Y, color= stat(level)))  +
  geom_sf(data= CDMX_alcaldias, inherit.aes = FALSE, fill= NA, color= "grey", size= 0.3) +
  facet_wrap(~TIPO.DE.EVENTO) +
  labs(subtitle= "CDMX: Localización espacial de accidentes viales por tipo de incidente (2019)",
       caption= "Fuente: elaboración propia a partir de datos.cdmx.gob.mx")
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#Se hicieron dos series de mapas, con puntos y con densidades. Se exhibe el mapa de densidades ya que es visualmente más claro, si bien del de puntos también se extrajeron conclusiones: #Los atropellos están extremadamente concentrados en el centro de la ciudad, casi de manera exclusiva en la alcaldía CUAUHTÉMOC. #Las caídas de ciclistas se concentran en la CUAUHTÉMOC, BENITO JUAREZ y la zona de MIGUEL HIDALGO correspondiente a Polanco/Chapultepec. Estos son los alcaldías mejor preparadas para el tránsito de ciclistas, donde se concentra todo el sistema Ecobici y mayor densidad de infraestructura de ciclovías. #Si bien con patrones de densidad distintos, los choques, derrapados y caídas de pasajeros están concentrados por toda la ciudad urbanizada. #Las volcaduras no se distinguen bien en el mapa de densidad, pero en un mapa de puntos se observa muy claramente que suceden a lo largo de toda la red de autopistas urbanas de CDMX. Es en estos ejes viales donde los vehículos circulan a velocidades que permiten este tipo de siniestros. Las volcaduras desaparecen allí donde no hay más autopistas.

#Visualización de la distribución espacial de “volcaduras”:

ggmap(CDMX) +
  geom_sf(data= CDMX_alcaldias, inherit.aes = FALSE, fill = NA) +
  geom_point(data= CDMX_incidentes_viales %>%
               filter(TIPO.DE.EVENTO == "VOLCADURA"), 
             aes(x= COORDENADA.X, y= COORDENADA.Y), color= "red", size= 0.1) +
  labs(x= "", y= "", title= "CDMX: Incidentes viales",
       subtitle= "Ubicación de volcaduras",
       caption= "Fuente: Elaboración propia a partir de datos.cdmx.gob.mx")
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

#Se observa que las volcaduras tienden a alinearse con los ejes viales principales y autopistas que figuran en el mapa. Muchas de estas autopistas son las demarcaciones de las alcaldías, y los puntos se tienden a alinear ahí.

#Veamos la serie de tiempo de accidentes según día de la semana:

CDMX_incidentes_viales$DIA_EVENTO <- factor(CDMX_incidentes_viales$DIA_EVENTO , levels= c("lunes", "martes", "miércoles", "jueves", "viernes", "sábado", "domingo"))

ggmap(CDMX) +
  geom_density_2d(data= CDMX_incidentes_viales %>%
                   filter(AÑO_EVENTO == "2019"),
                 aes(x= COORDENADA.X, y= COORDENADA.Y, color= stat(level))) +
  scale_colour_viridis_c()+
  facet_wrap(~DIA_EVENTO)+
  labs(subtitle= "CDMX: distribución espacial de totalidad incidentes viales (2019)")

Accidentes totales: Relativamente homogénea distribución de desnidad según día de la semana.

ggmap(CDMX) +
  geom_density_2d(data= CDMX_incidentes_viales %>%
                   filter(AÑO_EVENTO == "2019",
                          TIPO.DE.EVENTO == "ATROPELLADO"),
                 aes(x= COORDENADA.X, y= COORDENADA.Y, color= stat(level))) +
  facet_wrap(~DIA_EVENTO)+
  labs(subtitle= "CDMX: distribución espacial de choques vehiculares (2019)")

#Atropellados, choques y derrapados, sigue un patrón similar todos los días de la semana.

ggmap(CDMX) +
  geom_density_2d(data= CDMX_incidentes_viales %>%
                   filter(AÑO_EVENTO == "2019",
                          TIPO.DE.EVENTO == "CAIDA DE CICLISTA"),
                 aes(x= COORDENADA.X, y= COORDENADA.Y, color= stat(level))) +
  facet_wrap(~DIA_EVENTO) +
  scale_colour_viridis_c()+
  labs(subtitle= "CDMX: distribución espacial incidentes c/ ciclistas (2019)")

#Accidentes con ciclistas: esta es una serie que observa patrones distintos según el día de la semana. Claramente hay más densidad de accidentes los martes, miércoles y jueves, en los que las personas salen más en bicicleta por las alcaldías céntricas como medio de transporte. Los lunes y viernes hay poca movilización en bicicleta por lo que hay poca densidad de accidentes. En los fines de semana, el uso de la bicicleta cambia, pasa a ser más bien recreativo, y la densidad de accidentes se expande por todo el corredor que va desde Ciudad Universitaria en el sur de la ciudad, hasta la Basilica de Guadalupe en el norte.

#Veamos estos patrones por horarios de los choques vehículares:

ggmap(CDMX) +
    geom_density2d(data = CDMX_incidentes_viales %>%
                     filter(TIPO.DE.EVENTO == "CHOQUE",
                            AÑO_EVENTO == "2019",
                            DIA_EVENTO == "viernes" | DIA_EVENTO == "sábado" | DIA_EVENTO == "domingo"),
                   aes(x = COORDENADA.X, 
                       y = COORDENADA.Y, 
                       color = stat(level))) +
    scale_color_viridis_c() +
    facet_wrap(~hour(hm(HORA_EVENTO)), nrow = 4) +
  labs(subtitle= "Densidad espacial de choques, según hora del día (CDMX, 2019)")

#Vemos que la densidad de choques va cambiando y se va concentrando en ciertos rangos horarios en concreto, si bien la ciudad en su conjunto pareciera bastante afectada por este problema.