3 Visualizaciones

Cargamos paquetes

library(tidyverse)
library(gganimate)
library(gifski)
library(bookdown)
library(rmarkdown)
library(purrr)
library(lubridate)
library(stringr)
library(knitr)
#library(xts)
#library(zoo)
library(gridExtra)
#library(fpp2)
#library(RcppRoll)
#library(kableExtra)
options(knitr.table.format = "html")

Cargamos los datos

air_data_2 <- readRDS("data_rds/air_data_2.rds")

Echamos un vistazo a las variables

glimpse(air_data_2)
## Observations: 775,328
## Variables: 31
## $ station       <fct> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ station_name  <fct> Estación Avenida Constitución, Estación Avenida ...
## $ latitude      <dbl> 43.52981, 43.52981, 43.52981, 43.52981, 43.52981...
## $ longitude     <dbl> -5.673428, -5.673428, -5.673428, -5.673428, -5.6...
## $ date_time_utc <dttm> 2000-01-01 00:00:00, 2000-01-01 01:00:00, 2000-...
## $ SO2           <dbl> 23, 29, 40, 50, 39, 39, 40, 43, 38, 39, 48, 66, ...
## $ NO            <dbl> 89, 73, 53, 46, 35, 26, 27, 37, 30, 32, 33, 35, ...
## $ NO2           <dbl> 65, 60, 57, 53, 50, 49, 51, 48, 46, 39, 36, 39, ...
## $ CO            <dbl> 1.97, 1.61, 1.13, 1.06, 0.95, 0.82, 0.83, 0.88, ...
## $ PM10          <dbl> 53, 63, 56, 58, 50, 50, 57, 53, 40, 53, 54, 52, ...
## $ O3            <dbl> 9, 8, 7, 5, 6, 7, 7, 4, 5, 6, 10, 9, 11, 9, 18, ...
## $ wd            <dbl> 245, 222, 228, 239, 244, 218, 240, 248, 244, 235...
## $ ws            <dbl> 0.34, 1.06, 0.71, 0.84, 0.89, 0.71, 0.80, 1.11, ...
## $ TMP           <dbl> 5.7, 5.4, 5.3, 5.1, 4.6, 4.6, 4.5, 4.1, 3.8, 3.8...
## $ HR            <dbl> 76, 73, 72, 71, 72, 69, 68, 69, 70, 70, 64, 59, ...
## $ PRB           <dbl> 1026, 1025, 1025, 1025, 1024, 1024, 1024, 1026, ...
## $ RS            <dbl> 33, 33, 33, 33, 33, 33, 33, 33, 35, 112, 171, 36...
## $ LL            <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## $ BEN           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ TOL           <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ MXIL          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ PM25          <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ station_alias <fct> Constitucion, Constitucion, Constitucion, Consti...
## $ year          <dbl> 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, ...
## $ month         <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...
## $ date          <date> 2000-01-01, 2000-01-01, 2000-01-01, 2000-01-01,...
## $ week_day      <dbl> 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, ...
## $ hour          <int> 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14...
## $ holiday_type  <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ no_lab_days   <chr> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, ...
## $ wd_code       <fct> WS, SW, SW, WS, WS, SW, WS, WS, WS, SW, SW, SW, ...

Obtenemos las medias anuales de todos los contaminantes por estación y las plasmamos en un gráfico para observar sus evoluciones

# Calculamos las medias anuales
year_avgs <- air_data_2 %>% select(station_alias, date_time_utc, PM10, PM25, SO2, NO2, NO, O3, BEN, CO, MXIL, TOL) %>%
  group_by(station_alias, year = year(date_time_utc)) %>%
  summarise_all(funs(mean(., na.rm = TRUE))) %>% 
  select(-date_time_utc) # We drop this variable


# Convertimos los resultados a formato largo

year_avgs_long <- gather(year_avgs, contaminante, value, 3:length(year_avgs)) %>% 
                    filter(!(grepl('Constit', station_alias) 
                             & year == '2006' & 
                               contaminante %in% c('BEN', 'MXIL', 'TOL'))) %>%
                                # We filter this data because is only completed in 0.01%
                    filter(!(grepl('Constit', station_alias) &
                              year == '2008' & contaminante == 'PM25'))
                                # We filter this data because is only completed in 0.02%

# Finalmente representamos los indicadores generados en una tabla de gráficos

ggplot(year_avgs_long, aes(x = year, y = value)) + 
  geom_line() + 
  facet_grid(contaminante~station_alias,scales="free_y") +
   theme(axis.text = element_text(size = 6))

Nos quedamos sólo con el contaminante PM10

year_avgs_long_pm10 <- year_avgs_long %>% filter(contaminante == "PM10",
                                                 station_alias != "Montevil")

ggplot(year_avgs_long_pm10, aes(x = year, y = value)) + 
  geom_line() + 
  facet_grid(contaminante~station_alias,scales="free_y") +
   theme(axis.text = element_text(size = 6)) + ylim(0, 100)

ggplot(year_avgs_long_pm10, aes(x = year, y = value, group = station_alias, color = station_alias)) + 
  geom_line()  +
  ylim(0, 100)

La estación Argentina es la que presenta peores valores de PM10 a lo largo del tiempo. Echamos un vistazo a su evolución a través de la función ‘summary’. Para ello vamos a comparar dos histogramas, uno con las frecuencias promedio horarios de PM10 de los primeros 5 años de la serie y otro con las correspondientes a los últimos 5 años.

pm10_argentina_2000_2004 <- air_data_2 %>% 
                           select(date, 
                                   station_alias,
                                    PM10) %>%
                           filter(year(date) < "2005",
                                  station_alias == "Argentina")

pm10_argentina_2014_2018 <- air_data_2 %>% 
                           select(date, 
                                   station_alias,
                                    PM10) %>%
                           filter(year(date) > "2013",
                                  station_alias == "Argentina")


summary(pm10_argentina_2000_2004$PM10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   39.00   53.00   58.13   72.00  986.00     922
summary(pm10_argentina_2014_2018$PM10)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    0.00   16.00   24.00   28.84   35.00  522.00     330

Como vemos en los resultados, la evolución ha sido muy importante. Durante el periodo 2000-2004 el 75% de los promedios horarios se situaban por encima de 39 µg/m³. Mientras, en los últimos 5 años, el 75% de los promedios están por debajo de 35 µg/m³.

Vamos a comparar los histogramas de frecuencias de promedios horarios de PM10 en la estación Argentina en dos periodos, los primeros 5 años de la serie 2000-2004 contra los 5 últimos años.

graph_1 <- ggplot(pm10_argentina_2000_2004, aes(x = `PM10`)) +
                    geom_histogram() + 
                           xlim(0, 300) +
                           labs(title = "Distribución promedios horarios PM10 estación Argentina - (2000-2004)")


graph_2 <- ggplot(pm10_argentina_2014_2018, aes(x = `PM10`)) +
                    geom_histogram() + 
                           xlim(0, 300) +
                          labs(title = "Distribución promedios horarios PM10 estación Argentina - (2014-2018)")

grid.arrange(graph_1, 
             graph_2, 
             ncol = 1) 

pm10_argentina <- air_data_2 %>% 
                           select(date_time_utc, 
                                  station_alias,
                                  PM10) %>%
                           filter(station_alias == "Argentina",
                                  PM10 != "NA")

pm10_argentina_ranges <- pm10_argentina %>% 
                         mutate(rangos = if_else(PM10 <= 20, "<=20 µg/m³", 
                                         if_else(PM10 > 20 & PM10 <= 40, "(20-40] µg/m³", ">40 µg/m³")))

pm10_argentina_ranges$rangos <- factor(pm10_argentina_ranges$rangos, levels = c("<=20 µg/m³", "(20-40] µg/m³", ">40 µg/m³"))
 
pm10_argentina_ranges_conteo <- pm10_argentina_ranges %>%
                                group_by(year = year(date_time_utc), rangos) %>%
                                summarise(n = n()) %>%
                                ungroup() %>%
                                group_by(year) %>%
                                mutate(n_total_year = sum(n),
                                       prop_year = n / n_total_year)

graph_3 <- ggplot(pm10_argentina_ranges_conteo, aes(x = rangos, y = prop_year, fill = rangos)) +
                  geom_col() +
                  facet_grid(. ~year) +
                  scale_color_manual(values = c("<=20 µg/m³" = "green", "(20-40] µg/m³" = "orange", ">40 µg/m³" = "red"),
                                     aesthetics = c("colour", "fill"))
graph_4  <- ggplot(pm10_argentina_ranges_conteo, aes(x = rangos, y = prop_year, fill = rangos, label = rangos)) +
                  geom_col(show.legend = FALSE) +
                  geom_text(vjust=-0.5, size = 7, face = "bold", color = "grey")+
                  scale_color_manual(values = c("<=20 µg/m³" = "green", "(20-40] µg/m³" = "orange", ">40 µg/m³" = "red"),
                                     aesthetics = c("colour", "fill")) +
                  ggtitle("Evolución PM10 - Estación Argentina (Gijón)") +
                  ylab("% rangos PM10 - Promedios horarios") +
                  xlab("Rangos de niveles de PM10-Promedios horarios") +
                  theme(axis.text=element_text(size=12),
                  axis.title=element_text(size=16,face="bold", color = "grey"),
                  plot.title=element_text(size = 20, face="bold", color = "grey"),
                  plot.subtitle = element_text(size = 40, face="bold", color = "grey"),
                  axis.text.x=element_blank(),
                  axis.ticks.x = element_blank(),
                  panel.grid.major = element_blank(),
                  panel.grid.minor = element_blank(),
                  panel.background = element_blank()) +
                  scale_y_continuous(labels = scales::percent)

graph_4

Los porcentajes en el eje Y sale mal, pero por lo que sea al animarlo aparece correctamente. Así que lo dejamos así.

pm10_Argentina_2000_2018_anim <- graph_4 + transition_time(as.integer(year)) +
  labs(title = "Evolución PM10 estación Argentina (Gijón)",
        subtitle = "Año: {frame_time}")

pm10_Argentina_2000_2018_anim

anim_save("pm10_Argentina_2000_2018.gif", pm10_Argentina_2000_2018_anim)