M2U2 - Caso: Análisis de Patrones Puntuales de Eventos de Movilidad Waze

1 Introducción

El tramo de la Autopista Norte entre Bogotá, Cajicá y Tocancipá aparece repetido en los nombres de calle de la trama Waze (street) usada en este trabajo. Este documento revisa si los reportes de congestión en ese corredor se concentran en tramos fijos o se reparten al azar, comparando el jueves 26 y el viernes 27 de septiembre de 2024 con técnicas de análisis de patrones puntuales en R.

1.1 Objetivo

Establecer si la congestión vehicular del corredor se concentra en tramos fijos, si esa concentración cambia entre el jueves 26 y el viernes 27 de septiembre de 2024, y qué tramos deberían priorizarse en la gestión del tráfico.

1.2 Fuente de los datos

La trama corresponde a eventos Waze para Colombia, recortada al tramo de la Autopista Norte entre Bogotá, Cajicá y Tocancipá, durante el 26 y 27 de septiembre de 2024.

2 Datos y Limpieza

2.1 Carga del archivo

El archivo trae 17 columnas. Para que ninguna celda con formato de fecha corrompido se pierda al cargar, se fuerza el tipo numeric en las columnas de coordenadas: así, cualquier celda mal formateada como fecha llega a R como su número de serie original en vez de como texto inservible.

# Busca el archivo en varias rutas probables en lugar de asumir una sola.
# Si ninguna existe, el error indica el directorio de trabajo actual y
# todas las rutas revisadas, para localizar el archivo sin adivinar.
ruta_documento_actual <- tryCatch({
  if (requireNamespace("rstudioapi", quietly = TRUE) && rstudioapi::isAvailable()) {
    dirname(rstudioapi::getSourceEditorContext()$path)
  } else {
    NA_character_
  }
}, error = function(e) NA_character_)

buscar_archivo_datos <- function(nombre_archivo = "Trama_Waze.xlsx") {
  candidatos <- unique(c(
    nombre_archivo,
    file.path(getwd(), nombre_archivo),
    if (!is.na(ruta_documento_actual)) file.path(ruta_documento_actual, nombre_archivo),
    "C:/Pontificia_Un_Javeriana/Analisis_Geoespacial/M02_U02_Patrones/Trama_Waze.xlsx",
    file.path(Sys.getenv("USERPROFILE"), "Downloads", nombre_archivo),
    file.path(Sys.getenv("USERPROFILE"), "Documents", nombre_archivo),
    file.path(Sys.getenv("USERPROFILE"), "Desktop", nombre_archivo)
  ))
  candidatos <- candidatos[nzchar(candidatos)]
  encontrados <- candidatos[file.exists(candidatos)]

  if (length(encontrados) == 0) {
    stop(
      "\nNo se encontro '", nombre_archivo, "' en ninguna ruta candidata.\n",
      "Directorio de trabajo actual (getwd()): ", getwd(), "\n",
      "Rutas revisadas:\n  - ", paste(candidatos, collapse = "\n  - "), "\n",
      "Solucion: copia '", nombre_archivo, "' a la carpeta del .Rmd ",
      "(la mas simple), o agrega su ruta real dentro del vector 'candidatos' ",
      "en este chunk."
    )
  }

  message("Archivo de datos encontrado en: ", encontrados[1])
  encontrados[1]
}

ruta_excel <- buscar_archivo_datos("Trama_Waze.xlsx")

tipos_columnas <- c("numeric", "numeric", "text", "numeric", "logical",
                     "numeric", "numeric", "text", "text", "numeric",
                     "numeric", "text", "text", "numeric", "numeric",
                     "numeric", "text")

datos_waze <- read_excel(ruta_excel, sheet = "Hoja 1",
                          col_types = tipos_columnas)

kable(head(datos_waze, 5), caption = "Primeros registros de la trama Waze") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 12) %>%
  scroll_box(width = "100%")

Primeros registros de la trama Waze

id	waze_json_trama_id	country	reportRating	reportByMunicipalityUser	confidence	reliability	type	uuid	roadType	magvar	subtype	street	location_x	location_y	pubMillis	creation_Date
16	14	CO	2	FALSE	1	8	HAZARD	74a153fa-6ccd-4d6b-a94b-db01a88b002d	3	96	HAZARD_ON_SHOULDER_CAR_STOPPED	NA	-74016928	4938376	1.727315e+12	2024-09-26 01:53:49.600
17	14	CO	3	FALSE	4	10	HAZARD	b91961a4-e32c-4770-b2a0-551d7add5669	3	153	HAZARD_ON_SHOULDER_CAR_STOPPED	Variante Cajicá / RD45A Ramal A >(S)	-74016985	2733202	1.727313e+12	2024-09-26 01:53:49.600
18	14	CO	0	FALSE	4	10	HAZARD	b9262a77-b38b-46bd-b0b2-cda1b855a549	3	27	HAZARD_ON_SHOULDER_CAR_STOPPED	Bogotá-Tocancipá / RN55-01 >(N)	-73996251	4925537	1.727313e+12	2024-09-26 01:53:49.600
20	15	CO	2	FALSE	1	8	HAZARD	74a153fa-6ccd-4d6b-a94b-db01a88b002d	3	96	HAZARD_ON_SHOULDER_CAR_STOPPED	NA	-74016928	4938376	1.727315e+12	2024-09-26 01:54:02.100
21	15	CO	3	FALSE	4	10	HAZARD	b91961a4-e32c-4770-b2a0-551d7add5669	3	153	HAZARD_ON_SHOULDER_CAR_STOPPED	Variante Cajicá / RD45A Ramal A >(S)	-74016985	2733202	1.727313e+12	2024-09-26 01:54:02.100

2.2 Corrección de coordenadas

location_x y location_y están multiplicadas por 10^6, pero una parte de los registros perdió uno o dos dígitos en la exportación, lo que desplaza la coordenada un orden de magnitud. La función siguiente prueba potencias de diez hasta ubicar el valor en el rango geográfico plausible del corredor (longitud entre 73.9° y 74.1° oeste, latitud entre 4.0° y 5.2° norte), conservando el signo original.

corregir_magnitud <- function(valor, minimo, maximo) {
  base <- abs(valor) / 1e6
  signo <- sign(valor)
  resultado <- rep(NA_real_, length(valor))
  for (potencia in 0:3) {
    candidato <- base * 10^potencia
    pendiente <- is.na(resultado) & candidato >= minimo & candidato <= maximo
    resultado[pendiente] <- signo[pendiente] * candidato[pendiente]
  }
  resultado
}

datos_waze <- datos_waze %>%
  mutate(
    longitud = corregir_magnitud(location_x, 73.9, 74.1),
    latitud  = corregir_magnitud(location_y, 4.0, 5.2)
  )

res_total_filas     <- nrow(datos_waze)
res_filas_corruptas <- sum(is.na(datos_waze$longitud) | is.na(datos_waze$latitud))
res_pct_recuperado  <- formatC(100 * (1 - res_filas_corruptas / res_total_filas),
                                format = "f", digits = 1)

De 5070 registros, 74 no tenían un valor numérico recuperable detrás de la celda corrompida y quedaron fuera. El resto, 98.5%, se reubicó dentro del rango esperado sin necesidad de descartarlo.

Interpretación: perder un dígito en la exportación es un problema sistemático del archivo, no errores sueltos. Corregir por magnitud en vez de eliminar la fila evita perder información sobre tramos que solo aparecen en esos registros.

2.3 Ventana y eventos únicos

La ventana de estudio (ventana) cubre el tramo entre las coordenadas -74.04331° y -73.9929° de longitud y 4.885736° a 4.948562° de latitud, la franja donde aparecen los nombres de calle de Cajicá y Tocancipá en la trama. Cada evento se repite mientras está activo, así que se conserva un solo punto por evento (uuid), tomando su primer reporte de cada día.

ventana <- owin(xrange = c(-74.04331, -73.9929),
                 yrange = c(4.885736, 4.948562))

datos_waze <- datos_waze %>%
  mutate(fecha_hora = ymd_hms(creation_Date, quiet = TRUE),
         fecha = as.Date(fecha_hora))

datos_validos <- datos_waze %>%
  filter(!is.na(longitud), !is.na(latitud),
         longitud >= ventana$xrange[1], longitud <= ventana$xrange[2],
         latitud  >= ventana$yrange[1], latitud  <= ventana$yrange[2])

eventos_unicos <- datos_validos %>%
  filter(fecha %in% as.Date(c("2024-09-26", "2024-09-27"))) %>%
  arrange(fecha_hora) %>%
  distinct(uuid, fecha, .keep_all = TRUE)

tabla_resumen <- eventos_unicos %>%
  count(type, fecha, name = "eventos_unicos") %>%
  arrange(type, fecha)

kable(tabla_resumen, caption = "Eventos únicos por tipo y día") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 12)

Eventos únicos por tipo y día

type	fecha	eventos_unicos
ACCIDENT	2024-09-26	4
ACCIDENT	2024-09-27	6
HAZARD	2024-09-26	15
HAZARD	2024-09-27	28
JAM	2024-09-26	76
JAM	2024-09-27	121
ROAD_CLOSED	2024-09-26	3
ROAD_CLOSED	2024-09-27	3

res_pct_jam_total <- formatC(100 * sum(eventos_unicos$type == "JAM") / nrow(eventos_unicos),
                              format = "f", digits = 0)

Interpretación: la congestión (JAM) es el 77% de los 256 eventos únicos de los dos días, más que las otras tres categorías juntas. Cierres de vía y accidentes quedan con tan pocos puntos por día que cualquier prueba estadística sobre ellos debe leerse como exploratoria.

3 Patrón de Congestión

La congestión (JAM) recibe el análisis completo porque es el 77% de los eventos únicos registrados, la única categoría con suficientes puntos para sostener pruebas de aleatoriedad espacial con algo de potencia estadística. Las otras tres se revisan más adelante sin la misma profundidad.

3.1 Construcción del patrón de puntos

jam_26 <- eventos_unicos %>% filter(type == "JAM", fecha == as.Date("2024-09-26"))
jam_27 <- eventos_unicos %>% filter(type == "JAM", fecha == as.Date("2024-09-27"))

jam_26_ppp <- unique(ppp(x = jam_26$longitud, y = jam_26$latitud, window = ventana))
jam_27_ppp <- unique(ppp(x = jam_27$longitud, y = jam_27$latitud, window = ventana))

jam_26_sf <- ppp_a_sf(jam_26_ppp)
jam_27_sf <- ppp_a_sf(jam_27_ppp)

res_n_jam_26 <- npoints(jam_26_ppp)
res_n_jam_27 <- npoints(jam_27_ppp)
res_incremento_jam <- formatC(100 * (res_n_jam_27 / res_n_jam_26 - 1), format = "f", digits = 0)

El jueves se registraron 76 eventos únicos de congestión y el viernes 121, un 59% más. La sección de perfil horario revisa si ese aumento se concentra en alguna franja particular del día, lo que ayudaría a descartar o sostener una explicación de tipo fin de semana.

tm_basemap(mapa_base) +
  tm_shape(jam_26_sf) +
  tm_dots(col = paleta[1], size = 0.05) +
  tm_layout(title = "Congestión, jueves 26 de septiembre")

Figura 1. Patrón espacial de eventos de congestión sobre el mapa de Bogotá, jueves 26 de septiembre de 2024.

tm_basemap(mapa_base) +
  tm_shape(jam_27_sf) +
  tm_dots(col = paleta[2], size = 0.05) +
  tm_layout(title = "Congestión, viernes 27 de septiembre")

Figura 2. Patrón espacial de eventos de congestión sobre el mapa de Bogotá, viernes 27 de septiembre de 2024.

3.2 Perfil horario

Antes de entrar a las pruebas espaciales formales conviene saber a qué hora arrancó cada evento único de congestión, no cuántas veces se reportó mientras seguía activo, que ya se filtró en la deduplicación.

perfil_horario <- eventos_unicos %>%
  filter(type == "JAM") %>%
  mutate(hora = lubridate::hour(fecha_hora)) %>%
  count(fecha, hora)

hora_pico_26 <- perfil_horario %>% filter(fecha == as.Date("2024-09-26")) %>%
  filter(n == max(n)) %>% pull(hora) %>% min()
hora_pico_27 <- perfil_horario %>% filter(fecha == as.Date("2024-09-27")) %>%
  filter(n == max(n)) %>% pull(hora) %>% min()

res_hora_pico_26 <- hora_pico_26
res_hora_pico_27 <- hora_pico_27
res_n_pico_26 <- perfil_horario %>% filter(fecha == as.Date("2024-09-26"), hora == hora_pico_26) %>% pull(n)
res_n_pico_27 <- perfil_horario %>% filter(fecha == as.Date("2024-09-27"), hora == hora_pico_27) %>% pull(n)

ggplot(perfil_horario, aes(x = hora, y = n, color = factor(fecha))) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  scale_x_continuous(breaks = seq(0, 23, 2)) +
  scale_color_manual(values = paleta[1:2], labels = c("26 sept", "27 sept"),
                      name = "Día") +
  labs(title = "Eventos únicos de congestión por hora de inicio",
       x = "Hora del día", y = "Eventos únicos") +
  theme_minimal()

Figura 3. Hora de inicio de los eventos únicos de congestión, jueves y viernes.

Interpretación: el jueves la mayor cantidad de eventos arrancó a las 11:00, con 17 inicios. El viernes el pico se corrió a las 20:00, con 21 inicios. Ese corrimiento hacia la noche es compatible con un patrón de salida de fin de semana, aunque esta trama no incluye datos de movilidad inter-municipal para confirmarlo de forma directa.

4 Prueba CSR

4.1 Ji-cuadrado por cuadrantes

Se divide el corredor en una malla de 4 por 3 cuadrantes, dimensionada para que cada celda tenga en promedio más de cinco eventos esperados bajo aleatoriedad, condición que necesita la aproximación ji-cuadrado para ser válida.

prueba_26 <- quadrat.test(jam_26_ppp, nx = 4, ny = 3)
prueba_27 <- quadrat.test(jam_27_ppp, nx = 4, ny = 3)

res_chi_26 <- formatC(prueba_26$statistic, format = "f", digits = 2)
res_p_26   <- formatC(prueba_26$p.value, format = "f", digits = 4)
res_chi_27 <- formatC(prueba_27$statistic, format = "f", digits = 2)
res_p_27   <- formatC(prueba_27$p.value, format = "f", digits = 4)

# Disparidad real entre cuadrantes, solo contando los que tienen al menos
# un evento (los cuadrantes en cero suelen caer fuera del trazado vial y
# no aportan una comparación útil de intensidad).
conteo_26 <- as.vector(quadratcount(jam_26_ppp, nx = 4, ny = 3))
conteo_26_activos <- conteo_26[conteo_26 > 0]
res_cuadrante_max_26 <- max(conteo_26_activos)
res_cuadrante_min_26 <- min(conteo_26_activos)
res_cuadrante_razon_26 <- formatC(res_cuadrante_max_26 / res_cuadrante_min_26,
                                   format = "f", digits = 1)

plot(quadratcount(jam_26_ppp, nx = 4, ny = 3),
     main = "Conteo por cuadrantes, 26 de septiembre")
points(jam_26_ppp, col = paleta[1], pch = 16, cex = 0.6)

Figura 4. Conteo de eventos de congestión por cuadrante, jueves 26 de septiembre.

El estadístico ji-cuadrado es 65.47 el jueves y 99.07 el viernes, ambos con valor p de 0.0000 y 0.0000 respectivamente.

Interpretación: los dos valores p quedan muy por debajo de 0.05, así que se rechaza la aleatoriedad espacial completa los dos días. El cuadrante con más eventos el jueves tiene 18, frente a 2 en el que menos tiene entre los que registran al menos un evento, una razón de 9.0 a 1 que ya anticipa que la intensidad no es constante y que conviene usar la función K inhomogénea más adelante.

4.2 Vecino más cercano (función G)

La prueba de cuadrantes mira la intensidad por zonas, pero no dice nada sobre qué tan cerca está cada evento de su vecino más próximo. La función G y la distancia media al vecino más cercano completan ese diagnóstico a una escala más fina.

nn_26 <- nndist(jam_26_ppp) * m_por_grado
nn_27 <- nndist(jam_27_ppp) * m_por_grado

res_nn_media_26 <- formatC(mean(nn_26), format = "f", digits = 0)
res_nn_media_27 <- formatC(mean(nn_27), format = "f", digits = 0)
res_nn_mediana_26 <- formatC(median(nn_26), format = "f", digits = 0)
res_nn_mediana_27 <- formatC(median(nn_27), format = "f", digits = 0)

area_m2 <- area.owin(ventana) * m_por_grado^2
lambda_26 <- res_n_jam_26 / area_m2
lambda_27 <- res_n_jam_27 / area_m2
nn_csr_26 <- 1 / (2 * sqrt(lambda_26))
nn_csr_27 <- 1 / (2 * sqrt(lambda_27))

res_nn_csr_26 <- formatC(nn_csr_26, format = "f", digits = 0)
res_nn_csr_27 <- formatC(nn_csr_27, format = "f", digits = 0)
res_nn_ratio_26 <- formatC(100 * mean(nn_26) / nn_csr_26, format = "f", digits = 0)
res_nn_ratio_27 <- formatC(100 * mean(nn_27) / nn_csr_27, format = "f", digits = 0)

g_26 <- Gest(jam_26_ppp, correction = "km")
g_27 <- Gest(jam_27_ppp, correction = "km")

par(mfrow = c(1, 2))
plot(g_26, main = "G, 26 sept", col = c(paleta[1], "grey50"))
plot(g_27, main = "G, 27 sept", col = c(paleta[2], "grey50"))

par(mfrow = c(1, 1))

Figura 5. Función G observada contra la teórica de Poisson, jueves y viernes.

La distancia media al vecino más cercano es de 121 metros el jueves y 54 metros de mediana; bajo aleatoriedad completa, con 76 eventos en esta ventana, esa distancia esperada sería de 358 metros. El viernes la media observada es de 84 metros frente a 284 metros esperados.

Interpretación: la distancia observada es apenas el 34% de la esperada el jueves y el 30% el viernes. Los eventos están mucho más cerca entre sí de lo que predice el azar, y esa proporción se mantiene similar entre los dos días aunque el viernes tenga más eventos, lo que sugiere que el grado de agrupamiento relativo no cambia, solo su volumen.

5 Función K

5.1 Distancia de referencia

El radio máximo de análisis se calcula desde las dimensiones de la ventana, no desde las distancias entre puntos, para no calcular pares de distancias de forma cuadrática sobre cientos de eventos.

ancho_ventana <- diff(ventana$xrange)
alto_ventana  <- diff(ventana$yrange)
r_max <- 0.25 * sqrt(ancho_ventana^2 + alto_ventana^2)
r_seq <- seq(0, r_max, length.out = 128)

res_r_max_m <- formatC(r_max * m_por_grado, format = "f", digits = 0)

El radio máximo equivale a 2235 metros, una cuarta parte de la diagonal del corredor.

5.2 K homogénea

k_homog_26 <- Kest(jam_26_ppp, r = r_seq, correction = "Ripley")
k_homog_27 <- Kest(jam_27_ppp, r = r_seq, correction = "Ripley")

par(mfrow = c(1, 2))
plot(k_homog_26, main = "K homogénea, 26 sept", col = c(paleta[1], "grey50"))
plot(k_homog_27, main = "K homogénea, 27 sept", col = c(paleta[2], "grey50"))

par(mfrow = c(1, 1))

Figura 6. Función K homogénea de Ripley, jueves y viernes.

Interpretación: la curva observada queda por encima de la de Poisson en los dos días, lo que de nuevo apunta a agrupamiento. Pero como la prueba de cuadrantes ya mostró que la intensidad cambia por zona, parte de esta separación puede ser solo el efecto de tener más eventos en unas zonas que en otras, no agrupamiento adicional entre eventos cercanos.

5.3 K inhomogénea

Como la intensidad no es constante, se usa la función K inhomogénea, que compara cada evento contra una intensidad local en lugar de un promedio único para todo el corredor.

k_inhom_26 <- Kinhom(jam_26_ppp, r = r_seq, correction = "Ripley")
k_inhom_27 <- Kinhom(jam_27_ppp, r = r_seq, correction = "Ripley")

par(mfrow = c(1, 2))
plot(k_inhom_26, main = "K inhomogénea, 26 sept", col = c(paleta[1], "grey50"))
plot(k_inhom_27, main = "K inhomogénea, 27 sept", col = c(paleta[2], "grey50"))

par(mfrow = c(1, 1))

Figura 7. Función K inhomogénea de Ripley, jueves y viernes.

Interpretación: al controlar por la intensidad local, la curva se acerca más a la de Poisson que en la versión homogénea. Buena parte de lo que parecía agrupamiento es variación de la densidad de fondo del corredor, lo que ya habían anticipado tanto el ji-cuadrado como la función G.

6 Densidad Kernel

6.1 Selección del ancho de banda

bw_diggle_26 <- bw.diggle(jam_26_ppp)
bw_diggle_27 <- bw.diggle(jam_27_ppp)
bw_scott_26  <- bw.scott(jam_26_ppp)
bw_scott_27  <- bw.scott(jam_27_ppp)

tabla_bandwidth <- data.frame(
  dia = c("26 sept", "27 sept"),
  bw_diggle_m = c(formatC(bw_diggle_26 * m_por_grado, format = "f", digits = 0),
                   formatC(bw_diggle_27 * m_por_grado, format = "f", digits = 0)),
  bw_scott_x_m = c(formatC(bw_scott_26[1] * m_por_grado, format = "f", digits = 0),
                    formatC(bw_scott_27[1] * m_por_grado, format = "f", digits = 0)),
  bw_scott_y_m = c(formatC(bw_scott_26[2] * m_por_grado, format = "f", digits = 0),
                    formatC(bw_scott_27[2] * m_por_grado, format = "f", digits = 0))
)

kable(tabla_bandwidth,
      col.names = c("Día", "bw.diggle (m)", "bw.scott eje X (m)", "bw.scott eje Y (m)"),
      caption = "Comparación de anchos de banda") %>%
  kable_styling(bootstrap_options = c("striped", "hover"), font_size = 12)

Comparación de anchos de banda

Día	bw.diggle (m)	bw.scott eje X (m)	bw.scott eje Y (m)
26 sept	67	575	805
27 sept	37	555	751

res_bw_26_m <- formatC(bw_diggle_26 * m_por_grado, format = "f", digits = 0)
res_bw_27_m <- formatC(bw_diggle_27 * m_por_grado, format = "f", digits = 0)

res_bw_scott_26_m <- formatC(bw_scott_26[1] * m_por_grado, format = "f", digits = 0)
res_bw_diferencia_pct_26 <- formatC(100 * (bw_scott_26[1] - bw_diggle_26) / bw_diggle_26,
                                     format = "f", digits = 0)

Interpretación: bw.scott da 575 m el jueves, un 757% más ancho que 67 m de bw.diggle. Esa diferencia es consistente con que bw.scott asume una dispersión cercana a la normal bivariada, mientras bw.diggle ajusta por validación cruzada y responde mejor a un patrón agrupado sobre una vía angosta como este corredor, por lo que se usa para los mapas de intensidad (67 m el jueves, 37 m el viernes).

6.2 Mapas de intensidad

# Definidas aquí porque a partir de este punto cada mapa de intensidad
# se necesita tanto en su propia figura como, más adelante, en la
# diferencia entre los dos días.
convertir_a_raster <- function(imagen) {
  matriz <- as.matrix(imagen)
  capa <- raster::raster(matriz)
  raster::extent(capa) <- raster::extent(imagen$xrange[1], imagen$xrange[2],
                                          imagen$yrange[1], imagen$yrange[2])
  capa <- raster::flip(capa, direction = "y")
  raster::crs(capa) <- "+proj=longlat +datum=WGS84"
  capa
}

normalizar <- function(capa) {
  (capa - raster::cellStats(capa, "min")) /
    (raster::cellStats(capa, "max") - raster::cellStats(capa, "min"))
}

intensidad_26 <- density.ppp(jam_26_ppp, sigma = bw_diggle_26, edge = TRUE,
                              dimyx = c(128, 128))
raster_26 <- normalizar(convertir_a_raster(intensidad_26))
raster_26_vis <- enmascarar_bajo_umbral(raster_26)

tm_basemap(mapa_base) +
  tm_shape(raster_26_vis) +
  tm_raster(palette = c(paleta[5], paleta[1]), style = "cont", alpha = 0.85,
            title = "Intensidad normalizada") +
  tm_shape(jam_26_sf) +
  tm_dots(col = "black", size = 0.03) +
  tm_layout(title = "Densidad de congestión, 26 de septiembre")

Figura 8. Densidad kernel de eventos de congestión sobre el mapa de Bogotá, jueves 26 de septiembre.

intensidad_27 <- density.ppp(jam_27_ppp, sigma = bw_diggle_27, edge = TRUE,
                              dimyx = c(128, 128))
raster_27 <- normalizar(convertir_a_raster(intensidad_27))
raster_27_vis <- enmascarar_bajo_umbral(raster_27)

tm_basemap(mapa_base) +
  tm_shape(raster_27_vis) +
  tm_raster(palette = c(paleta[5], paleta[2]), style = "cont", alpha = 0.85,
            title = "Intensidad normalizada") +
  tm_shape(jam_27_sf) +
  tm_dots(col = "black", size = 0.03) +
  tm_layout(title = "Densidad de congestión, 27 de septiembre")

Figura 9. Densidad kernel de eventos de congestión sobre el mapa de Bogotá, viernes 27 de septiembre.

Interpretación: en las dos figuras la franja de mayor intensidad corre sobre la vía que cruza el centro de Cajicá y sigue hacia Calahorra, no se reparte por el resto del rectángulo de la ventana. Eso es justamente lo que el ji-cuadrado y la función K ya habían señalado con números: el problema no es la zona en general, es ese tramo vial puntual, y por eso al enmascarar las celdas de baja densidad el mapa base queda visible en casi toda la ventana salvo ahí.

Nota de verificación: la orientación del raster depende de cómo spatstat almacena la matriz de intensidad. Si al ver la Figura 8 el área de mayor intensidad no coincide con la nube de puntos negros, retirar la línea flip() dentro de convertir_a_raster() en lugar de aplicarla.

7 Comparación 26-27

7.1 Normalización y diferencia

Cada mapa de intensidad ya se normalizó entre 0 y 1 de forma independiente al construir las Figuras 8 y 9, para que la comparación no se vea afectada por tener 76 eventos un día y 121 el otro.

7.2 Diferencia celda a celda

diferencia_raster <- raster_27 - raster_26

res_max_expansion   <- formatC(raster::cellStats(diferencia_raster, "max"), format = "f", digits = 2)
res_max_contraccion <- formatC(raster::cellStats(diferencia_raster, "min"), format = "f", digits = 2)

# Para el mapa, se oculta el cambio casi nulo (|diferencia| < 0.05) y se
# deja ver solo dónde la congestión relativa subió o bajó de forma
# apreciable entre los dos días.
diferencia_raster_vis <- diferencia_raster
diferencia_raster_vis[abs(diferencia_raster) < 0.05] <- NA

El cambio celda a celda va de -0.92 (mayor contracción relativa) a 0.66 (mayor expansión relativa) en la escala normalizada.

Interpretación: los valores positivos marcan tramos donde la congestión relativa creció el viernes frente al jueves, y los negativos donde bajó. Los puntos que aparecen en ambos días bajo una y otra intensidad, sin depender del cambio de calendario, son los que más sugieren un cuello de botella estructural; los que solo aparecen el viernes están más ligados al efecto de fin de semana que al diseño de la vía.

7.3 Mapas finales

tm_basemap(mapa_base) +
  tm_shape(raster_26_vis) +
  tm_raster(palette = c(paleta[5], paleta[1]), style = "cont", alpha = 0.85,
            title = "Intensidad normalizada") +
  tm_layout(title = "Congestión, jueves 26")

Figura 10. Mapa 1, intensidad normalizada de congestión sobre el mapa de Bogotá, jueves 26 de septiembre.

tm_basemap(mapa_base) +
  tm_shape(raster_27_vis) +
  tm_raster(palette = c(paleta[5], paleta[2]), style = "cont", alpha = 0.85,
            title = "Intensidad normalizada") +
  tm_layout(title = "Congestión, viernes 27")

Figura 11. Mapa 2, intensidad normalizada de congestión sobre el mapa de Bogotá, viernes 27 de septiembre.

tm_basemap(mapa_base) +
  tm_shape(diferencia_raster_vis) +
  tm_raster(palette = c(paleta[1], "white", paleta[2]), style = "cont",
            midpoint = 0, alpha = 0.85, title = "Cambio relativo") +
  tm_layout(title = "Diferencia 27 menos 26 de septiembre")

Figura 12. Mapa 3, diferencia celda a celda sobre el mapa de Bogotá, entre el viernes y el jueves.

8 Otros Eventos

Cierres de vía, peligros y accidentes se revisan de forma exploratoria. Con cierres de vía y accidentes el número de eventos únicos por día es tan bajo que solo se ubican en el mapa, sin prueba estadística.

construir_ppp_tipo <- function(tipo, dia) {
  sub <- eventos_unicos %>% filter(type == tipo, fecha == as.Date(dia))
  if (nrow(sub) == 0) return(NULL)
  unique(ppp(x = sub$longitud, y = sub$latitud, window = ventana))
}

hazard_26 <- construir_ppp_tipo("HAZARD", "2024-09-26")
hazard_27 <- construir_ppp_tipo("HAZARD", "2024-09-27")
cerrado_26 <- construir_ppp_tipo("ROAD_CLOSED", "2024-09-26")
cerrado_27 <- construir_ppp_tipo("ROAD_CLOSED", "2024-09-27")
accidente_26 <- construir_ppp_tipo("ACCIDENT", "2024-09-26")
accidente_27 <- construir_ppp_tipo("ACCIDENT", "2024-09-27")

res_n_hazard_26 <- npoints(hazard_26)
res_n_hazard_27 <- npoints(hazard_27)
res_n_cerrado_26 <- npoints(cerrado_26)
res_n_cerrado_27 <- npoints(cerrado_27)
res_n_accidente_26 <- npoints(accidente_26)
res_n_accidente_27 <- npoints(accidente_27)

Peligros en vía pasaron de 15 eventos el jueves a 28 el viernes, una muestra suficiente para una prueba CSR orientativa. Cierres de vía (3 y 3) y accidentes (4 y 6) son demasiado escasos para cualquier prueba formal.

prueba_hazard_26 <- quadrat.test(hazard_26, nx = 2, ny = 2)
prueba_hazard_27 <- quadrat.test(hazard_27, nx = 2, ny = 2)

res_p_hazard_26 <- formatC(prueba_hazard_26$p.value, format = "f", digits = 4)
res_p_hazard_27 <- formatC(prueba_hazard_27$p.value, format = "f", digits = 4)

El valor p de la prueba ji-cuadrado para peligros en vía es 0.0034 el jueves y 0.0000 el viernes.

Interpretación: con muestras de 15 y 28 eventos, esta prueba es orientativa. Sirve para no descartar la categoría de entrada, no para sostener una conclusión estadística firme sobre peligros en vía.

construir_sf_etiquetado <- function(objeto_ppp, etiqueta) {
  if (is.null(objeto_ppp) || npoints(objeto_ppp) == 0) return(NULL)
  sf_obj <- ppp_a_sf(objeto_ppp)
  sf_obj$evento <- etiqueta
  sf_obj
}

otros_eventos_lista <- list(
  construir_sf_etiquetado(cerrado_26, "Cierre de vía, 26 sept"),
  construir_sf_etiquetado(cerrado_27, "Cierre de vía, 27 sept"),
  construir_sf_etiquetado(accidente_26, "Accidente, 26 sept"),
  construir_sf_etiquetado(accidente_27, "Accidente, 27 sept")
)
otros_eventos_lista <- Filter(Negate(is.null), otros_eventos_lista)
otros_eventos_sf <- do.call(rbind, otros_eventos_lista)

# El modo "view" de tmap dibuja siempre marcadores circulares y no
# soporta combinar color y forma en una sola capa de puntos; por eso
# tipo y día se combinan en una sola categoría de color en vez de
# usar dos estéticas (color y forma) al mismo tiempo.
tm_basemap(mapa_base) +
  tm_shape(otros_eventos_sf) +
  tm_dots(col = "evento",
          palette = c("Cierre de vía, 26 sept" = paleta[3],
                      "Cierre de vía, 27 sept" = paleta[4],
                      "Accidente, 26 sept" = paleta[1],
                      "Accidente, 27 sept" = paleta[2]),
          size = 0.08, title = "Evento y día") +
  tm_layout(title = "Cierres de vía y accidentes, 26-27 de septiembre")

Figura 13. Cierres de vía y accidentes sobre el mapa de Bogotá, jueves y viernes, sin prueba estadística por tamaño de muestra.

9 Conclusiones

9.1 Hallazgos principales

El ji-cuadrado, la distancia al vecino más cercano y la función K coinciden en lo mismo desde tres ángulos distintos: el corredor no reparte la congestión al azar, se concentra sobre la vía que cruza Cajicá y sigue hacia Calahorra, en tramos que se repiten el jueves y el viernes. El viernes el número de eventos sube 59% y el pico horario se corre de las 11:00 a las 20:00, un patrón compatible con la salida hacia Cajicá y Tocancipá antes del fin de semana.

9.2 Recomendaciones para planificación urbana

El aumento de 59% en eventos del viernes, concentrado en la franja de la noche, justifica medidas de tráfico acotadas a esa ventana horaria, como agentes o semaforización adaptativa entre las 20:00 y el cierre de la noche, sin necesidad de mantenerlas el resto de la semana. Los tramos que se repiten en ambas jornadas, visibles donde coinciden las Figuras 10 y 11 antes de calcular la diferencia, son los que justifican una intervención permanente de infraestructura, porque su congestión no depende del día.

9.3 Limitaciones

El análisis cubre dos días y un corredor específico, por lo que no debe extenderse a toda Bogotá sin más evidencia. Cierres de vía y accidentes necesitan una ventana de tiempo más larga antes de poder analizarse con el mismo rigor estadístico que la congestión.