Indicadores
Cargar librerías
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.0 ✔ readr 2.1.5
## ✔ ggplot2 3.5.1 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(lubridate)
library(ggthemes)
library(showtext)
## Loading required package: sysfonts
## Loading required package: showtextdb
library(sysfonts)
Cargar archivos
registros_2022 <- read.csv("data_2022.csv")
registros_2023 <- read.csv("data_2023.csv")
registros_2024 <- read.csv("data_2024.csv")
registros <- rbind(registros_2022, registros_2023, registros_2024)
limpiar la maestra
registros <- registros %>%
mutate(
estacion = case_when(
Sensor_id == "ANL1" ~ "Guadalupe",
Sensor_id == "ANL2" ~ "San Nicolás",
Sensor_id == "ANL3" ~ "Santa Catarina",
Sensor_id == "ANL4" ~ "San Pedro",
Sensor_id == "ANL5" ~ "Universidad",
Sensor_id == "ANL6" ~ "García",
Sensor_id == "ANL7" ~ "San Bernabé",
Sensor_id == "ANL8" ~ "Cadereyta",
Sensor_id == "ANL9" ~ "Escobedo",
Sensor_id == "ANL10" ~ "Apodaca",
Sensor_id == "ANL11" ~ "Monterrey (sur)",
Sensor_id == "ANL12" ~ "Obispado",
Sensor_id == "ANL13" ~ "Juárez",
Sensor_id == "ANL15" ~ "Pesquería",
Sensor_id == "ANL16" ~ "San Juan",
TRUE ~ NA_character_
)
)
Base de prueba
prueba <- registros %>%
select(estacion, Dia, PM10, PM25, O3, CO, NO1, NO2, NOx, SO2) %>%
pivot_longer(cols = PM10:SO2,
names_to = "contaminante",
values_to = "valor") %>%
mutate(
Dia = ymd_hms(Dia), # Use ymd_hms instead of dmy_hm
Dia = Dia + hours(6),
fecha = as.Date(Dia), # Extract only the date
hora = format(Dia, "%H:%M:%S"), # Extract the time in HH:MM:SS
valor = as.numeric(valor) # Convert 'valor' to numeric
) %>%
select(estacion, fecha, hora, contaminante, valor)
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `valor = as.numeric(valor)`.
## Caused by warning:
## ! NAs introduced by coercion
# Contar registros válidos (no NA) por fecha y contaminante
conteo_validos <- prueba %>%
filter(!is.na(valor)) %>% # Filtra valores no nulos
group_by(estacion, fecha, contaminante) %>%
summarise(n_registros = n(), .groups = "drop") # Cuenta registros válidos
# Unir conteo con la base original y asignar "suficiencia"
prueba <- prueba %>%
left_join(conteo_validos, by = c("estacion", "fecha", "contaminante")) %>%
mutate(suficiencia = ifelse(n_registros >= 18, TRUE, FALSE)) %>%
mutate(cumple_norma = ifelse(
suficiencia == TRUE, TRUE, TRUE
)) %>%
select(estacion, fecha, hora, contaminante, valor, suficiencia) # Mantener orden
pm2_5 <- prueba %>%
filter(contaminante == "PM25") %>%
mutate(valor = if_else(valor < 0, NA_real_, valor))
pm10_norma <- prueba %>%
filter(contaminante == "PM10") %>% # Trabajamos solo con PM10
group_by(estacion, fecha) %>% # Agrupamos por estación y fecha
summarise(
media_pm10 = mean(valor, na.rm = TRUE),
# Suponemos que dentro de cada grupo (día y estación para PM10)
# el valor de 'suficiencia' es el mismo en todas las filas, así que tomamos el primero.
suficiencia = first(suficiencia),
.groups = "drop"
) %>%
mutate(
# Si no se cumplen los 18 registros, se asigna "ND" sin evaluar el promedio
cumple_norma = if_else(
!suficiencia,
"ND",
# Si suficiencia es TRUE, se evalúa el promedio:
if_else(media_pm10 <= 60, "CUMPLE", if_else(media_pm10 > 60, "NO CUMPLE", NA_character_))
)
)
# Se une la clasificación calculada para PM10 a la base original 'prueba'
prueba <- prueba %>%
left_join(
pm10_norma %>% select(estacion, fecha, cumple_norma_pm10 = cumple_norma),
by = c("estacion", "fecha")
) %>%
mutate(
# Asignamos la clasificación solo para PM10; para otros contaminantes se deja NA
cumple_norma = if_else(contaminante == "PM10", cumple_norma_pm10, NA_character_)
) %>%
select(estacion, fecha, hora, contaminante, valor, suficiencia, cumple_norma)
PM 2.5
pm25_norma <- pm2_5 %>%
group_by(estacion, fecha) %>%
summarise(
media_pm25 = mean(valor, na.rm = TRUE),
maximo_pm25 = max(valor, na.rm = TRUE),
n_registros = sum(!is.na(valor)), # Cuenta solo registros válidos
suficiencia = first(suficiencia),
.groups = "drop"
) %>%
mutate(
suficiencia = replace_na(suficiencia, FALSE),
cumple_norma = if_else(
!suficiencia,
"ND",
if_else(media_pm25 <= 25, "CUMPLE", "NO CUMPLE")
)
) %>%
mutate(
AQI_cat = case_when(
media_pm25 < 25 ~ "Buena",
media_pm25 < 45 ~ "Aceptable",
media_pm25 < 79 ~ "Mala",
media_pm25 < 147 ~ "Muy Mala",
TRUE ~ "Extremadamente Mala"
)
)
## Warning: There were 2836 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `maximo_pm25 = max(valor, na.rm = TRUE)`.
## ℹ In group 283: `estacion = "Apodaca"` and `fecha = 2022-10-25`.
## Caused by warning in `max()`:
## ! no non-missing arguments to max; returning -Inf
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 2835 remaining warnings.
write.csv(pm25_norma, "pm25_max_norma.csv", row.names = FALSE)
Visualizaciones
Visualización 1
pm25_norma %>%
filter(year(fecha) == 2024, cumple_norma == "NO CUMPLE") %>%
group_by(estacion) %>%
summarise(dias_sobre_norma = n()) %>%
ggplot(aes(x = estacion, y = dias_sobre_norma)) +
# Gráfica de barras en color verde oscuro
geom_col(fill = "darkgreen") +
labs(
title = "Total de días sobre la norma de PM2.5 por estación (2024)",
x = "Estación",
y = "Número de días NO CUMPLE"
) +
theme_minimal() +
# Rotar etiquetas del eje X para que no se sobrepongan
theme(
axis.text.x = element_text(angle = 45, hjust = 1)
)
Visualización 2
pm25_norma %>%
filter(year(fecha) == 2024, cumple_norma == "CUMPLE") %>%
group_by(estacion) %>%
summarise(dias_cumplen_norma = n()) %>%
ggplot(aes(x = estacion, y = dias_cumplen_norma)) +
geom_col(fill = "darkgreen") +
labs(
title = "Total de días que SÍ cumplen la norma de PM25 por estación (2024)",
x = "Estación",
y = "Número de días CUMPLE"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1)
)
Visualización HeatMap versión 1
pm25_norma %>%
# 1. Filtrar sólo los días que NO cumplen la norma en 2024
filter(year(fecha) == 2024, cumple_norma == "NO CUMPLE") %>%
# 2. Crear una columna con el mes (abreviado)
mutate(mes = month(fecha, label = TRUE, abbr = TRUE)) %>%
# 3. Agrupar y contar días
group_by(estacion, mes) %>%
summarise(dias_no_cumple = n(), .groups = "drop") %>%
# 4. Convertir 'estacion' a factor ordenado alfabéticamente
# y luego se usa rev() para que la primera quede arriba.
mutate(estacion = factor(estacion, levels = rev(sort(unique(estacion))))) %>%
# 5. (Opcional) Ordenar los meses de Jan a Dec, si quieres un orden específico
mutate(mes = factor(mes,
levels = c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"))) %>%
# 6. Crear el heatmap
ggplot(aes(x = mes, y = estacion, fill = dias_no_cumple)) +
geom_tile() +
geom_text(aes(label = dias_no_cumple), color = "black") +
scale_fill_gradient(low = "white", high = "tomato") +
labs(
x = "Mes",
y = "Estación",
fill = "Días NO cumple",
title = "Total de días sobre la norma para el contaminante PM2.5"
) +
theme_minimal() +
theme(legend.position = "none")
#scale_fill_viridis_c(direction = -1, option = "rocket") +
Visualización HeatMap versión 3
# 1. Calculamos los días NO CUMPLE por estacion y mes
df <- pm25_norma %>%
filter(year(fecha) == 2024, cumple_norma == "NO CUMPLE") %>%
mutate(mes = month(fecha, label = TRUE, abbr = TRUE)) %>%
group_by(estacion, mes) %>%
summarise(dias_no_cumple = n(), .groups = "drop")
# 2. Calculamos el total por mes (suma de todas las estaciones)
df_totales_mes <- df %>%
group_by(estacion) %>%
summarise(dias_no_cumple = sum(dias_no_cumple)) %>%
mutate(mes = "Total") # Etiquetamos la "categoría" de estación como TOTAL
# 3. Unimos
df_final <- bind_rows(df, df_totales_mes)
# 4. Aseguramos que "TOTAL" sea un nivel al final de la variable estacion
df_final$mes <- factor(
df_final$mes,
levels = c(levels(df$mes), "Total")
)
df_final$estacion <- factor(
df_final$estacion,
levels = rev(sort(unique(df_final$estacion)))
)
# 5. Graficamos
ggplot(df_final, aes(x = mes, y = estacion, fill = log(dias_no_cumple))) + # transformción logarítmica para mantener
# 1. Capa de "tiles" coloreados
geom_tile() +
# 2. Capa de texto con la cuenta de días
geom_text(aes(label = dias_no_cumple), color = "black") +
# Escala de color (puedes ajustarla a tu gusto)
scale_fill_viridis_c(direction = -1, option = "rocket", alpha = 0.80) +
# scale_fill_gradient(low = "white", high = "red") +
labs(
x = "Mes",
y = "Estación",
fill = "Días NO cumple"
) +
theme_minimal() +
theme(legend.position = "none")
Función Gráfica
plot_estacion <- function(estacion_param) {
pm25_norma %>%
filter(year(fecha) == 2024) %>%
filter(suficiencia == TRUE) %>%
filter(estacion == estacion_param) %>%
ggplot(aes(x = fecha, y = media_pm25, fill = AQI_cat)) +
geom_col() +
scale_fill_manual(
values = c(
"Buena" = "green",
"Aceptable" = "yellow",
"Mala" = "orange",
"Muy Mala" = "red",
"Extremadamente Mala"= "purple"
)
) +
scale_x_date(
limits = c(as.Date("2024-01-01"), as.Date("2024-12-31")),
date_breaks = "1 month",
date_labels = "%b"
) +
labs(
title = paste("Evolución diaria de la calidad del aire (PM2.5) - 2024 -", estacion_param),
subtitle = "Cada barra representa el valor promedio diario de PM2.5",
x = "Fecha",
y = expression(paste("PM"[2.5], " (µg/m³)")),
fill = "Categoría"
) +
theme_minimal() +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none"
)
}
Gráfica por estación
Apodaca <- plot_estacion("Apodaca")
Cadereyta <- plot_estacion("Cadereyta")
Escobedo <- plot_estacion("Escobedo")
Guadalupe <- plot_estacion("Guadalupe")
Juarez <- plot_estacion("Juárez") # Asegúrate de que coincida con tu base de datos
Monterrey_sur <- plot_estacion("Monterrey (sur)")
Obispado <- plot_estacion("Obispado")
SanBernabe <- plot_estacion("San Bernabé")
SanNicolas <- plot_estacion("San Nicolás")
SanPedro <- plot_estacion("San Pedro")
SantaCatarina <- plot_estacion("Santa Catarina")
Universidad <- plot_estacion("Universidad")
Matriz de gráficas
# Con patchwork, concatenas usando + y luego indicas layout
#matriz_graficas <-
# Apodaca + Cadereyta + Escobedo + Guadalupe +
# Juarez + Monterrey_sur + Obispado + SanBernabe +
# SanNicolas + SanPedro + SantaCatarina + Universidad +
# plot_layout(ncol = 2) # 4 columnas => 3 filas
# Visualizas la matriz:
#matriz_graficas
# Si usaste patchwork y guardaste en 'matriz_graficas'
#ggsave("matriz_estaciones.png", matriz_graficas, width = 16, height = 14, units = "in")
Estadísticas
1) ¿Cuántos días no cumplen con la suficiencia (por estación)?
dias_no_suficiencia_por_estacion <- pm25_norma %>%
filter(!suficiencia) %>%
filter(year(fecha) == 2024) %>%
group_by(estacion) %>%
summarise(
dias_no_suficiencia = n(),
.groups = "drop"
)
dias_no_suficiencia_por_estacion
## # A tibble: 15 × 2
## estacion dias_no_suficiencia
## <chr> <int>
## 1 Apodaca 164
## 2 Cadereyta 51
## 3 Escobedo 26
## 4 García 169
## 5 Guadalupe 209
## 6 Juárez 29
## 7 Monterrey (sur) 87
## 8 Obispado 156
## 9 Pesquería 364
## 10 San Bernabé 326
## 11 San Juan 364
## 12 San Nicolás 61
## 13 San Pedro 55
## 14 Santa Catarina 135
## 15 Universidad 67
2) Promedio de datos recopilados por día (por estación y
general)
promedio_datos_dia_estacion <- pm25_norma %>%
group_by(estacion) %>%
filter(year(fecha) == 2024) %>%
summarise(
promedio_datos = mean(n_registros, na.rm = TRUE),
.groups = "drop"
)
promedio_datos_dia_estacion
## # A tibble: 15 × 2
## estacion promedio_datos
## <chr> <dbl>
## 1 Apodaca 13.7
## 2 Cadereyta 21.0
## 3 Escobedo 21.6
## 4 García 14.0
## 5 Guadalupe 12.7
## 6 Juárez 21.9
## 7 Monterrey (sur) 18.7
## 8 Obispado 14.8
## 9 Pesquería 0
## 10 San Bernabé 3.91
## 11 San Juan 0.239
## 12 San Nicolás 20.4
## 13 San Pedro 20.9
## 14 Santa Catarina 15.9
## 15 Universidad 20.3
promedio_datos_dia_general <- pm25_norma %>%
summarise(
promedio_datos = mean(n_registros, na.rm = TRUE)
)
promedio_datos_dia_general
## # A tibble: 1 × 1
## promedio_datos
## <dbl>
## 1 16.1
3) Estación con el mayor porcentaje de insuficiencia
estacion_mayor_pct_insuficiencia <- pm25_norma %>%
group_by(estacion) %>%
filter(year(fecha) == 2024) %>%
summarise(
total_dias = n(),
dias_no_suficiencia = sum(!suficiencia),
pct_insuficiencia = 100 * dias_no_suficiencia / total_dias,
.groups = "drop"
) %>%
arrange(desc(pct_insuficiencia)) %>%
slice(1)
estacion_mayor_pct_insuficiencia
## # A tibble: 1 × 4
## estacion total_dias dias_no_suficiencia pct_insuficiencia
## <chr> <int> <int> <dbl>
## 1 Pesquería 364 364 100
top5_insuficiencia <- pm25_norma %>%
group_by(estacion) %>%
filter(year(fecha) == 2024) %>%
summarise(
total_dias = n(),
dias_no_suficiencia = sum(!suficiencia),
pct_insuficiencia = 100 * dias_no_suficiencia / total_dias,
.groups = "drop"
) %>%
arrange(desc(pct_insuficiencia)) %>%
slice(1:5)
top5_insuficiencia
## # A tibble: 5 × 4
## estacion total_dias dias_no_suficiencia pct_insuficiencia
## <chr> <int> <int> <dbl>
## 1 Pesquería 364 364 100
## 2 San Juan 364 364 100
## 3 San Bernabé 364 326 89.6
## 4 Guadalupe 364 209 57.4
## 5 García 364 169 46.4