Análisis exploratorio del parque automotor vs la calidad del aire

Colombia 2013-2017

======================================================================================================================

Lectura inicial

# Lectura habitual
datos <- read.csv("ParqueAutomotor_Aire/DATOS_DE_CALIDAD_DEL_AIRE_EN_COLOMBIA_2011-2017.csv")

# Con biblioteca data.table
datos <- data.table::fread(file = "DATOS_DE_CALIDAD_DEL_AIRE_EN_COLOMBIA_2011-2017.csv",
                           sep = ",", dec = ".")

Datos de calidad de aire

Conversión a formato .Rdata

# Conversión a formato .Rdata
save(datos, file = "aire.Rdata", compress = "xz")

Cargando datos en formato .Rdata

load("ParqueAutomotor_Aire/aire.Rdata")

Estructura interna

str(datos)

'data.frame':   15657064 obs. of  16 variables:
 $ Fecha                         : Factor w/ 61368 levels "01/01/2011 01:00:00 a. m.",..: 46965 46968 46946 46948 46950 46952 46954 46956 46958 46960 ...
 $ Autoridad.Ambiental           : Factor w/ 27 levels "AMVA","CAM","CAR",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Nombre.de.la.estación         : Factor w/ 337 levels "Aeropuerto","Albania",..: 15 15 15 15 15 15 15 15 15 15 ...
 $ Tecnología                    : Factor w/ 2 levels "Automática","Manual": 1 1 1 1 1 1 1 1 1 1 ...
 $ Latitud                       : num  6.33 6.33 6.33 6.33 6.33 ...
 $ Longitud                      : num  -75.6 -75.6 -75.6 -75.6 -75.6 ...
 $ Código.del.departamento       : int  5 5 5 5 5 5 5 5 5 5 ...
 $ Departamento                  : Factor w/ 23 levels "ANTIOQUIA","ARAUCA",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Código.del.municipio          : int  5088 5088 5088 5088 5088 5088 5088 5088 5088 5088 ...
 $ Nombre.del.municipio          : Factor w/ 150 levels "AGUSTÍN CODAZZI",..: 16 16 16 16 16 16 16 16 16 16 ...
 $ Tipo.de.estación              : Factor w/ 2 levels "Fija","Indicativa": 1 1 1 1 1 1 1 1 1 1 ...
 $ Tiempo.de.exposición          : int  1 1 1 1 1 1 1 1 1 1 ...
 $ Variable                      : Factor w/ 20 levels "CO","Dirección del Viento",..: 8 8 8 8 8 8 8 8 8 8 ...
 $ Unidades                      : Factor w/ 9 levels "%","°","°C","m/s",..: 9 9 9 9 9 9 9 9 9 9 ...
 $ Concentración                 : num  4.88 22.41 48.72 95.86 116.58 ...
 $ Nueva.columna.georreferenciada: Factor w/ 329 levels "(1.216489, -77.282944)",..: 234 234 234 234 234 234 234 234 234 234 ...

Depuración de datos

Datos de aire por departamento

# Bibliotecas
library(dplyr)
library(lubridate)

# Depuración
aire <- datos %>% 
  mutate(Fecha = as.Date(as.character(Fecha), format = "%d/%m/%Y"),
         Año = year(Fecha)) %>% 
  select(Año, Departamento, Nombre.del.municipio, Tiempo.de.exposición,
         Variable, Unidades, Concentración) %>% 
  rename(mpio = Nombre.del.municipio, texpo = Tiempo.de.exposición)

aire <- aire %>% 
  mutate(Año = factor(Año))

# Datos por departamento
filtro <- c("CO", "NO", "NO2", "O3", "PM10", "PM2.5", "PST", "SO2")
depto_aire <- aire %>%
  filter(Variable %in% filtro) %>% 
  group_by(Año, Departamento, Variable, Unidades) %>% 
  summarise(promedio = round(mean(Concentración), digits = 2),
            minimo = round(min(Concentración), digits = 2),
            maximo = round(max(Concentración), digits = 2),
            mediana = round(median(Concentración), digits = 2),
            p98 = round(quantile(Concentración, probs = 0.98), digits = 2))
head(depto_aire)

# Guardando datos de aire por departamento
save(depto_aire, file = "ParqueAutomotor_Aire/Aire_depto.Rdata")

Datos de aire por municipio

# Datos por municipio
mpio_aire <- aire %>%
  filter(Variable %in% filtro) %>% 
  group_by(Año, mpio, Variable, Unidades) %>% 
  summarise(promedio = round(mean(Concentración), digits = 2),
            minimo = round(min(Concentración), digits = 2),
            maximo = round(max(Concentración), digits = 2),
            mediana = round(median(Concentración), digits = 2),
            p98 = round(quantile(Concentración, probs = 0.98), digits = 2))
head(mpio_aire)

# Guardando datos de aire por municipio
save(mpio_aire, file = "ParqueAutomotor_Aire/Aire_mpio.Rdata")

Datos parque automotor

Por departamento

Datos 2013

auto13 <- data.table::fread(file = "ParqueAutomotor_Aire/Parque_Automotor_2013.csv",
                               sep = ",", dec = ".")

# Datos por departamento
auto13_depto <- auto13 %>%
  group_by(Departamento) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2013", 32))

auto13_depto

Datos 2014

auto14 <- data.table::fread(file = "ParqueAutomotor_Aire/Parque_Automotor_2014.csv",
                               sep = ",", dec = ".")

# Datos por departamento
auto14_depto <- auto14 %>%
  group_by(Departamento) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2014", 32))

auto14_depto

Datos 2015

auto15 <- data.table::fread(file = "ParqueAutomotor_Aire/Parque_Automotor_2015.csv",
                               sep = ",", dec = ".")

# Datos por departamento
auto15_depto <- auto15 %>%
  group_by(Departamento) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2015", 32))

auto15_depto

Datos 2016

auto16 <- data.table::fread(file = "ParqueAutomotor_Aire/Parque_Automotor_2016.csv",
                               sep = ",", dec = ".")

# Datos por departamento 
auto16_depto <- auto16 %>%
  group_by(Departamento) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2016", 32))

auto16_depto

Datos 2017

auto17 <- data.table::fread(file = "ParqueAutomotor_Aire/Parque_Automotor_2017.csv",
                               sep = ",", dec = ".")

# Datos por departamento 
auto17_depto <- auto17 %>%
  group_by(Departamento) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2017", 32))

auto17_depto

Parque automotor por departamento 2013-2017

autos_depto <- rbind(auto13_depto, auto14_depto, auto15_depto, auto16_depto,
                     auto17_depto)
save(autos_depto, file = "ParqueAutomotor_Aire/Autos_depto.Rdata")
load("ParqueAutomotor_Aire/Autos_depto.Rdata")
autos_depto

Por municipio

Datos 2013

# Datos por Ciudad
auto13_Ciudad <- auto13 %>%
  group_by(Ciudad) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2013", n()))

auto13_Ciudad

Datos 2014

# Datos por Ciudad en 
auto14_Ciudad <- auto14 %>%
  group_by(Ciudad) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2014", n()))

auto14_Ciudad

Datos 2015

# Datos por Ciudad en 
auto15_Ciudad <- auto15 %>%
  group_by(Ciudad) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2015", n()))

auto15_Ciudad

Datos 2016

# Datos por Ciudad
auto16_Ciudad <- auto16 %>%
  group_by(Ciudad) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2016", n()))

auto16_Ciudad

Datos 2017

auto17_Ciudad <- auto17 %>%
  group_by(Ciudad) %>% 
  summarise(Total = n()) %>% 
  mutate(Año = rep("2017", n()))

auto17_Ciudad

Parque automotor por municipio 2013-2017

autos_ciudad <- rbind(auto13_Ciudad, auto14_Ciudad, auto15_Ciudad, auto16_Ciudad,
                      auto17_Ciudad)
save(autos_ciudad, file = "ParqueAutomotor_Aire/Autos_ciudad.Rdata")
load("ParqueAutomotor_Aire/Autos_ciudad.Rdata")
autos_ciudad

Unión de datos

Por departamento

# Lectura de datos de aire por departamento y año
load("ParqueAutomotor_Aire/Aire_depto.Rdata")
depto_aire <- depto_aire %>% 
  ungroup() %>% 
  mutate(Departamento = tolower(Departamento),
         Departamento = gsub("á", "a", Departamento),
         Departamento = gsub("é", "e", Departamento),
         Departamento = gsub("í", "i", Departamento),
         Departamento = gsub("ó", "o", Departamento),
         Departamento = gsub("ú", "u", Departamento))

# Lectura de datos parque automotor por departamento y año
load("ParqueAutomotor_Aire/Autos_depto.Rdata")
autos_depto <- autos_depto %>% 
  mutate(Departamento = tolower(Departamento),
         Departamento = gsub("á", "a", Departamento),
         Departamento = gsub("é", "e", Departamento),
         Departamento = gsub("í", "i", Departamento),
         Departamento = gsub("ó", "o", Departamento),
         Departamento = gsub("ú", "u", Departamento))

# Unión de datos
df_depto <- inner_join(depto_aire, autos_depto, by = c("Departamento", "Año"))

# Guardando y cargando datos df_depto
save(df_depto, file = "ParqueAutomotor_Aire/df_depto.Rdata")
load("ParqueAutomotor_Aire/df_depto.Rdata")

df_depto

Por municipio

# Lectura de datoa de aire por municipio y año
load("ParqueAutomotor_Aire/Aire_mpio.Rdata")
mpio_aire <- mpio_aire %>% 
  ungroup() %>% 
  mutate(mpio = tolower(mpio),
         mpio = gsub("á", "a", mpio),
         mpio = gsub("é", "e", mpio),
         mpio = gsub("í", "i", mpio),
         mpio = gsub("ó", "o", mpio),
         mpio = gsub("ú", "u", mpio))

# Lectura de datos parque automotor por municipio y año
load("ParqueAutomotor_Aire/Autos_ciudad.Rdata")
autos_mpio <- autos_ciudad %>%
  rename(mpio = Ciudad) %>% 
  mutate(mpio = tolower(mpio),
         mpio = gsub("á", "a", mpio),
         mpio = gsub("é", "e", mpio),
         mpio = gsub("í", "i", mpio),
         mpio = gsub("ó", "o", mpio),
         mpio = gsub("ú", "u", mpio))

# Unión de datos
df_mpio <- inner_join(mpio_aire, autos_mpio, by = c("mpio", "Año"))

# Guardando y cargando datos df_mpio
save(df_mpio, file = "ParqueAutomotor_Aire/df_mpio.Rdata")
load("ParqueAutomotor_Aire/df_mpio.Rdata")
df_mpio

Gráficos exploratorios

Por departamento

# Datos necesarios para graficar
load("ParqueAutomotor_Aire/df_depto.Rdata") # datos completos por departamento

# Ajustando datos por departamento
df_depto <- df_depto %>% 
  mutate(Año = factor(Año),
         Departamento = factor(Departamento)) %>% 
  droplevels()

• Estructura interna:

str(df_depto)

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   295 obs. of  10 variables:
 $ Año         : Factor w/ 5 levels "2013","2014",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Departamento: Factor w/ 22 levels "antioquia","arauca",..: 1 1 1 1 1 1 1 5 6 6 ...
 $ Variable    : Factor w/ 8 levels "CO","NO","NO2",..: 1 2 3 4 5 6 7 5 5 6 ...
 $ Unidades    : Factor w/ 1 level "µg/m3": 1 1 1 1 1 1 1 1 1 1 ...
 $ promedio    : num  1397.1 24 31.4 33.6 51.5 ...
 $ minimo      : num  115 0 0 0 1 ...
 $ maximo      : num  8363 412 670 242 472 ...
 $ mediana     : num  1145.6 11 28.2 23.4 46.5 ...
 $ p98         : num  4238.7 122 84.9 118.2 120 ...
 $ Total       : int  130860 130860 130860 130860 130860 130860 130860 12466 18462 18462 ...

Parque automotor vs calidad de aire promedio

# Biblioteca ggplot2, tidyr y RColorBrewer
library(ggplot2)
library(tidyr)

# Identificación de filas con datos atípicos
which(df_depto$Variable == "NO" & df_depto$promedio >= 40)    # fila 231
which(df_depto$Variable == "NO2" & df_depto$promedio >= 100)  # fila 291
which(df_depto$Variable == "O3" & df_depto$promedio >= 100)   # fila 243
which(df_depto$Variable == "PM10" & df_depto$promedio >= 100) # fila 280
which(df_depto$Variable == "PM2.5" & df_depto$promedio >= 40) # fila 239
which(df_depto$Variable == "SO2" & df_depto$promedio >= 500)  # fila 264

# Vector de filas a seleccionar
filas <- c(which(df_depto$Variable == "NO" & df_depto$promedio >= 40),
           which(df_depto$Variable == "NO2" & df_depto$promedio >= 100),
           which(df_depto$Variable == "O3" & df_depto$promedio >= 100),
           which(df_depto$Variable == "PM10" & df_depto$promedio >= 100),
           which(df_depto$Variable == "PM2.5" & df_depto$promedio >= 40),
           which(df_depto$Variable == "SO2" & df_depto$promedio >= 500))

# Relación del parque automotor vs calidad de aire
df_depto %>%
  slice(-filas)%>% 
  ggplot(data = ., mapping = aes(x = Total, y = promedio)) +
  facet_wrap(facets = ~Variable, scales = "free_y") +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, aes(color = "blue"), lwd = 0.7) +
  geom_smooth(method = "loess", se = FALSE,  aes(color = "red"), lwd = 0.7) +
  scale_color_identity(guide = "legend",
                       name = "Modelo:",
                       breaks = c("blue", "red"),
                       labels = c("Lineal", "Loess")) +
  labs(x = "Número de carros (n)", y = "Promedio µg/m3",
       title = "Relación del parque automotor vs calidad de aire",
       subtitle = "Colombia 2013-2017", 
       caption = "Cada punto representa un departamento.") +
  theme_linedraw() +
  theme(legend.position = "bottom")

Parque automotor vs calidad de aire Percentil 98

df_depto %>%
  slice(-filas)%>% 
  ggplot(data = ., mapping = aes(x = Total, y = p98)) +
  facet_wrap(facets = ~Variable, scales = "free_y") +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, aes(color = "blue"), lwd = 0.7) +
  geom_smooth(method = "loess", se = FALSE,  aes(color = "red"), lwd = 0.7) +
  scale_color_identity(guide = "legend",
                       name = "Modelo:",
                       breaks = c("blue", "red"),
                       labels = c("Lineal", "Loess")) +
  labs(x = "Número de carros (n)", y = "Percentil 98 µg/m3",
       title = "Relación del parque automotor vs calidad de aire",
       subtitle = "Colombia 2013-2017", 
       caption = "Cada punto representa un departamento.") +
  theme_linedraw() +
  theme(legend.position = "bottom")
  

Por municipio

# Datos necesarios para graficar
load("ParqueAutomotor_Aire/df_mpio.Rdata")  # datos completos por municipio

# Ajustando datos por municipio
df_mpio <- df_mpio %>% 
  mutate(Año = factor(Año),
         mpio = factor(mpio)) %>% 
  droplevels()

• Estructura interna:

# Estructura interna df_depto
str(df_mpio)

Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   632 obs. of  10 variables:
 $ Año     : Factor w/ 5 levels "2013","2014",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ mpio    : Factor w/ 73 levels "agustin codazzi",..: 2 2 6 7 7 9 9 9 9 10 ...
 $ Variable: Factor w/ 8 levels "CO","NO","NO2",..: 5 7 5 4 5 2 3 4 5 1 ...
 $ Unidades: Factor w/ 1 level "µg/m3": 1 1 1 1 1 1 1 1 1 1 ...
 $ promedio: num  29.2 59.9 20.6 33.6 29.1 ...
 $ minimo  : num  5.24 19.74 10 0.2 10.13 ...
 $ maximo  : num  77.9 110.7 40 188.8 54.3 ...
 $ mediana : num  27.1 56.4 20 29.1 28.5 ...
 $ p98     : num  69.4 105 30.8 83.4 52 ...
 $ Total   : int  197 197 6470 780 780 1713 1713 1713 1713 11000 ...

Parque automotor vs calidad de aire promedio

# Identificación de filas con datos atípicos
which(df_mpio$Variable == "NO" & df_mpio$promedio >= 90)    
which(df_depto$Variable == "NO2" & df_mpio$promedio >= 250)
which(df_depto$Variable == "O3" & df_mpio$promedio >= 500)
which(df_depto$Variable == "PM10" & df_mpio$promedio >= 100)
which(df_depto$Variable == "SO2" & df_mpio$promedio >= 500)

# Vector de filas a seleccionar
filas2 <- c(which(df_mpio$Variable == "NO" & df_mpio$promedio >= 90),
           which(df_depto$Variable == "NO2" & df_mpio$promedio >= 250),
           which(df_depto$Variable == "O3" & df_mpio$promedio >= 500),
           which(df_depto$Variable == "PM10" & df_mpio$promedio >= 100),
           which(df_depto$Variable == "SO2" & df_mpio$promedio >= 500))

# Relación del parque automotor vs calidad de aire
df_mpio %>%
  slice(-filas2) %>%   
  ggplot(data = ., mapping = aes(x = Total, y = promedio)) +
  facet_wrap(facets = ~Variable, scales = "free_y") +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, aes(color = "blue"), lwd = 0.7) +
  geom_smooth(method = "loess", se = FALSE,  aes(color = "red"), lwd = 0.7) +
  scale_color_identity(guide = "legend",
                       name = "Modelo:",
                       breaks = c("blue", "red"),
                       labels = c("Lineal", "Loess")) +
  labs(x = "Número de carros (n)", y = "Promedio µg/m3",
       title = "Relación del parque automotor vs calidad de aire",
       subtitle = "Colombia 2013-2017", 
       caption = "Cada punto representa un municipio.") +
  theme_linedraw() +
  theme(legend.position = "bottom")

Parque automotor vs calidad de aire Percentil 98

df_mpio %>%
  slice(-filas2)%>% 
  ggplot(data = ., mapping = aes(x = Total, y = p98)) +
  facet_wrap(facets = ~Variable, scales = "free_y") +
  geom_point() +
  geom_smooth(method = lm, se = FALSE, aes(color = "blue"), lwd = 0.7) +
  geom_smooth(method = "loess", se = FALSE,  aes(color = "red"), lwd = 0.7) +
  scale_color_identity(guide = "legend",
                       name = "Modelo:",
                       breaks = c("blue", "red"),
                       labels = c("Lineal", "Loess")) +
  labs(x = "Número de carros (n)", y = "Percentil 98 µg/m3",
       title = "Relación del parque automotor vs calidad de aire",
       subtitle = "Colombia 2013-2017", 
       caption = "Cada punto representa un municipio.") +
  theme_linedraw() +
  theme(legend.position = "bottom")s