<jupyter_start><jupyter_text>INFLUENCIA DE LAS PRECIPITACIONES Y EL VIENTO EN LA CONTAMINACIÓN  1. Introducción La comunidad de Madrid ha estado incumpliendo durante años los niveles establecidos por los parámetros de la calidad del aire. Hasta la fecha existen pocos datos sobre como influye la meteorología en relación a la contaminación. Las teorías actuales sugieren que las emisiones contaminantes de los vehículos son la causa principal, reduciendo o limitando la circulación de los mismos. Una de las medidas que el Ayuntamiento de Madrid tomo dentrode su plan de calidad del aire es la aplicación de restricciones de tráfico de paso al centro de la capital(Madrid central).La contaminación supone un problema grave para la salud, según la Organización Mundial de la Salud (OMS) "la contaminación atmosférica urbana aumenta el riesgo de padecer enfermedades respiratorias agudas, como la neumonía y crónicas, como el cáncer de pulmón y las enfermedades cardiovasculares".Nuestro objetivo en el presente trabajo es comprobar si otras causas como la meteorología pueden influir en la calidad del aire en la ciudad de Madrid.  2. Objetivos  2.1 Hipótesis La contaminación en la comunidad de Madrid se ve reducida cuando aumenta el nivel de precipitaciones y la velocidad del viento.   2.2 Restricciones - La precisión del estudio es a nivel diario ya que no hay disponibles fuentes de datos para la precipitación a nivel horario.- Unicamente se estudia la relación de la precipitación y velocidad del viento con la contaminación atmosférica   y no se incluyen otros tipos como la contaminación, acústica, lumínica, etc.  2.3 Alcance El estudio se centra en la ciudad de Madrid (España)  3. Metodología  3.1 Explicación de la metodología La metodología elegida es SEMMA, es  el  acrónimo  a  las  cinco  fases:  (Sample,  Explore,  Modify, Model, Assess). A continuación se explica cada una de las fases:  **Sample**: Se toma una muestra de los datos y se crea con ellos un dataset (más) reducido. Las muestras deben ser lo suficientemente grandes como para contener toda la información significativa pero también lo suficientemente pequeñas como para poder procesarlas rápidamente.  **Explore**: Se analizan (exploran) los datos para comprenderlos, buscando relaciones previstas, tendencias imprevistas y/o anomalías, así como ideas para su posterior procesamiento.  **Modify**: Se modifican los datos (creando, seleccionando y transformando las variables adecuadas) para facilitar el proceso siguiente (modelización).  **Model**: Se entrenan y aplican distintos modelos para el procesamiento de los datos, seleccionando aquel cuya tecnología y/o parámetros produce unos resultados aparentemente más adecuados.  **Assess**: Se evalúan los resultados del proceso de minería de datos, midiendo su utilidad y fiabilidad.   ![semma.png](attachment:semma.png)  3.2 Planificación A continuación se visualiza la planificación llevada a cabo por los integrantes del proyecto con un Diagrama de Gantt:  [>>Enlace al Diagrama Gant<<](https://drive.google.com/file/d/192o2H7ss1-ZGJ13-0Reet4vwuaRe5dlM/view)  4. Implementación Durante la elaboración de este estudio se utilizará la metodología SEMMA como ya se ha citado (Sample,  Explore,  Modify, Model, Assess), por ello la implementación seguirá las fases de dicha metodología.  Instalamos Paquetes y cargamos librerias<jupyter_code>#Limpiamos memoria
rm(list=ls());
is.installed <- function(mypkg) {
  is.element(mypkg, installed.packages()[,1]) 
}

#Instalamos paquetes para agrupar csv
if (!is.installed("data.table")){
  install.packages("data.table",repos = "http://cran.us.r-project.org")
};
#Intslamos paquete para manipular ts
if (!is.installed("xts")){
  install.packages("xts",repos = "http://cran.us.r-project.org")
};
#Instalamos paquete para agrupar datos
if (!is.installed("dplyr")){
  install.packages("dplyr",repos = "http://cran.us.r-project.org")
};
# Instalamos paquete para leer exccels
if (!is.installed("readxl")){
  install.packages("readxl",repos = "http://cran.us.r-project.org")
};
#Instalamos paquetes para modelo predictivo de arbol
if (!is.installed("party")){
  install.packages("party",repos = "http://cran.us.r-project.org")
};
if (!is.installed("rpart")){
  install.packages("rpart",repos = "http://cran.us.r-project.org")
};
if (!is.installed("randomForest")){
  install.packages("randomForest",repos = "http://cran.us.r-project.org")
};
if (!is.installed("rpart.plot")){
  install.packages("rpart.plot",repos = "http://cran.us.r-project.org")
};
if (!is.installed("rattle")){
  install.packages("rattle",repos = "http://cran.us.r-project.org")
};
if (!is.installed("RColorBrewer")){
  install.packages("RColorBrewer",repos = "http://cran.us.r-project.org")
};
if (!is.installed("tidyr")){
  install.packages("tidyr",repos = "http://cran.us.r-project.org")
};
#Instalamos paquetes para graficas de correlaciones
if (!is.installed("corrplot")){
  install.packages("corrplot",repos = "http://cran.us.r-project.org")
};
# Cargamos las librerias
library(data.table)
library(xts)
library(dplyr)    # loads %>%
library(readxl)
library(ggplot2)
library(party)
library(rpart)
library(randomForest)
library(rpart.plot)
library(rattle)
library(RColorBrewer)
library(tidyr)
library(corrplot)<jupyter_output><empty_output><jupyter_text>1º FASE: Sample Se toma una muestra de los datos y se crea con ellos un dataset (más) reducido. Las muestras deben ser lo suficientemente grandes como para contener toda la información significativa pero también lo suficientemente pequeñas como para poder procesarlas rápidamente.  Cargamos datos calidad del aire 2001-2018 Madrid Cargamos las mediciones de los contaminantes de Madrid entre 2001-2018 y creamos un dataframe con todas ellas<jupyter_code>#Cargamos los datos de los anos en un dataset (utilziamos read.csv ya que el separador es,)
data2001<- read.csv(file = "data/madrid_2001.csv",header = TRUE);
data2002<- read.csv(file = "data/madrid_2002.csv",header = TRUE);
data2003<- read.csv(file = "data/madrid_2003.csv",header = TRUE);
data2004<- read.csv(file = "data/madrid_2004.csv",header = TRUE);
data2005<- read.csv(file = "data/madrid_2005.csv",header = TRUE);
data2006<- read.csv(file = "data/madrid_2006.csv",header = TRUE);
data2007<- read.csv(file = "data/madrid_2007.csv",header = TRUE);
data2008<- read.csv(file = "data/madrid_2008.csv",header = TRUE);
data2009<- read.csv(file = "data/madrid_2009.csv",header = TRUE);
data2010<- read.csv(file = "data/madrid_2010.csv",header = TRUE);
data2011<- read.csv(file = "data/madrid_2011.csv",header = TRUE);
data2012<- read.csv(file = "data/madrid_2012.csv",header = TRUE);
data2013<- read.csv(file = "data/madrid_2013.csv",header = TRUE);
data2014<- read.csv(file = "data/madrid_2014.csv",header = TRUE);
data2015<- read.csv(file = "data/madrid_2015.csv",header = TRUE);
data2016<- read.csv(file = "data/madrid_2016.csv",header = TRUE);
data2017<- read.csv(file = "data/madrid_2017.csv",header = TRUE);
data2018<- read.csv(file = "data/madrid_2018.csv",header = TRUE);
# unimos los datos 
data <- data.table::rbindlist(list(data2001,data2002,data2003,data2004,data2005,data2006,
                                   data2007,data2008,data2009,data2010,data2011,data2012,
                                   data2013,data2014,data2015,data2016,data2017,data2018),
                                   fill = TRUE);

# convertimos la lista en un dataset
airMadrid <- as.data.frame(data); 

# liberamos memoria
rm(list=ls(pattern="data"))<jupyter_output><empty_output><jupyter_text>Carga datos Precipitaciones 2001-2018 Cargamos las mediciones de las precipitaciones de Madrid entre 2001-2018 y creamos un dataframe con ellas<jupyter_code>#Cargamos los datos de las precipitaciones
rainfall_2001<- readxl::read_xlsx("data/rainfall_2001.xlsx",col_names = TRUE);
rainfall_2002<- readxl::read_xlsx("data/rainfall_2002.xlsx",col_names = TRUE);
rainfall_2003<- readxl::read_xlsx("data/rainfall_2003.xlsx",col_names = TRUE);
rainfall_2004<- readxl::read_xlsx("data/rainfall_2004.xlsx",col_names = TRUE);
rainfall_2005<- readxl::read_xlsx("data/rainfall_2005.xlsx",col_names = TRUE);
rainfall_2006<- readxl::read_xlsx("data/rainfall_2006.xlsx",col_names = TRUE);
rainfall_2007<- readxl::read_xlsx("data/rainfall_2007.xlsx",col_names = TRUE);
rainfall_2008<- readxl::read_xlsx("data/rainfall_2008.xlsx",col_names = TRUE);
rainfall_2009<- readxl::read_xlsx("data/rainfall_2009.xlsx",col_names = TRUE);
rainfall_2010<- readxl::read_xlsx("data/rainfall_2010.xlsx",col_names = TRUE);
rainfall_2011<- readxl::read_xlsx("data/rainfall_2011.xlsx",col_names = TRUE);
rainfall_2012<- readxl::read_xlsx("data/rainfall_2012.xlsx",col_names = TRUE);
rainfall_2013<- readxl::read_xlsx("data/rainfall_2013.xlsx",col_names = TRUE);
rainfall_2014<- readxl::read_xlsx("data/rainfall_2014.xlsx",col_names = TRUE);
rainfall_2015<- readxl::read_xlsx("data/rainfall_2015.xlsx",col_names = TRUE);
rainfall_2016<- readxl::read_xlsx("data/rainfall_2016.xlsx",col_names = TRUE);
rainfall_2017<- readxl::read_xlsx("data/rainfall_2017.xlsx",col_names = TRUE);
rainfall_2018<- readxl::read_xlsx("data/rainfall_2018.xlsx",col_names = TRUE);
# unimos los datos
rainfall <- data.table::rbindlist(list(rainfall_2001,rainfall_2002,rainfall_2003,rainfall_2004,
              rainfall_2005,rainfall_2006,rainfall_2007,rainfall_2008,rainfall_2009,rainfall_2010,
              rainfall_2011,rainfall_2012,rainfall_2013,rainfall_2014,rainfall_2015,rainfall_2016,
rainfall_2017,rainfall_2018), fill = TRUE)

# convertimos la lista en un dataset
dataFrameRainfall <- as.data.frame(rainfall);

# liberamos memoria
rm(list=ls(pattern="rainfall"))<jupyter_output><empty_output><jupyter_text>Carga datos Viento 2001-2018 Cargamos  las mediciones del viento de Madrid entre 2001-2018 y creamos un dataframe con ellas<jupyter_code>#Cargamos los datos de las precipitaciones
wind_2001<- readxl::read_xlsx("data/wind_2001.xlsx",col_names = TRUE);
wind_2002<- readxl::read_xlsx("data/wind_2002.xlsx",col_names = TRUE);
wind_2003<- readxl::read_xlsx("data/wind_2003.xlsx",col_names = TRUE);
wind_2004<- readxl::read_xlsx("data/wind_2004.xlsx",col_names = TRUE);
wind_2005<- readxl::read_xlsx("data/wind_2005.xlsx",col_names = TRUE);
wind_2006<- readxl::read_xlsx("data/wind_2006.xlsx",col_names = TRUE);
wind_2007<- readxl::read_xlsx("data/wind_2007.xlsx",col_names = TRUE);
wind_2008<- readxl::read_xlsx("data/wind_2008.xlsx",col_names = TRUE);
wind_2009<- readxl::read_xlsx("data/wind_2009.xlsx",col_names = TRUE);
wind_2010<- readxl::read_xlsx("data/wind_2010.xlsx",col_names = TRUE);
wind_2011<- readxl::read_xlsx("data/wind_2011.xlsx",col_names = TRUE);
wind_2012<- readxl::read_xlsx("data/wind_2012.xlsx",col_names = TRUE);
wind_2013<- readxl::read_xlsx("data/wind_2013.xlsx",col_names = TRUE);
wind_2014<- readxl::read_xlsx("data/wind_2014.xlsx",col_names = TRUE);
wind_2015<- readxl::read_xlsx("data/wind_2015.xlsx",col_names = TRUE);
wind_2016<- readxl::read_xlsx("data/wind_2016.xlsx",col_names = TRUE);
wind_2017<- readxl::read_xlsx("data/wind_2017.xlsx",col_names = TRUE);
wind_2018<- readxl::read_xlsx("data/wind_2018.xlsx",col_names = TRUE);
# unimos los datos
wind <- data.table::rbindlist(list(wind_2001,wind_2002,wind_2003,wind_2004,
              wind_2005,wind_2006,wind_2007,wind_2008,wind_2009,wind_2010,
              wind_2011,wind_2012,wind_2013,wind_2014,wind_2015,wind_2016,
wind_2017,wind_2018), fill = TRUE)

dataFrameWind <- as.data.frame(wind);

rm(list=ls(pattern="wind"))

wind <-dataFrameWind
rm(dataFrameWind)<jupyter_output><empty_output><jupyter_text>2º FASE: Explore Se analizan (exploran) los datos para comprenderlos, buscando relaciones previstas, tendencias imprevistas y/o anomalías, así como ideas para su posterior procesamiento.Para cada uno de los dataframe creados se visualizará el principo y el final de los mismos para observar si el intervalo en el que se encuentran abarca 2001-2018  Visualización inicial datos calidad del aire<jupyter_code>head(airMadrid)
tail(airMadrid)<jupyter_output><empty_output><jupyter_text>Visualización inicial datos precipitaciones<jupyter_code>head(dataFrameRainfall)
tail(dataFrameRainfall)<jupyter_output><empty_output><jupyter_text>Visualización inicial datos viento<jupyter_code>head(wind)
tail(wind)<jupyter_output><empty_output><jupyter_text>3º FASE: Modify Se modifican los datos (creando, seleccionando y transformando las variables adecuadas) para facilitar el proceso siguiente (modelización).  Modificaciones dataframe precipitaciones Se cambiarán los nombres de las columnas para una mejor utilización de las mismas. También se modificará el tipo de la columna "*Precipitación: l/m2*" a númerico pra poder tratar los datos.<jupyter_code># cambiamos nombre de las columnas
data.table::setnames(dataFrameRainfall, "FECHA", "date")
data.table::setnames(dataFrameRainfall, "Precipitación: l/m2", "quantity")
#Convertimos tipos
dataFrameRainfall <- as.data.frame(dataFrameRainfall)
dataFrameRainfall$quantity <- as.numeric(dataFrameRainfall$quantity)<jupyter_output><empty_output><jupyter_text>Modificaciones dataframe viento Se cambiarán los nombres de las columnas para una mejor utilización de las mismas. Se eliminarán las columnas que no se utilizarán en el estudio y se comprobará el tipo de las columnas restantes, transformando la velocidad del viento a m/s para su posterior discretización.<jupyter_code># cambiamos nombre de las columnas
data.table::setnames(wind, "FECHA", "date")
data.table::setnames(wind, "Veloc. Media (Km/h)", "Speed")
data.table::setnames(wind, "Racha Max (Km/h)", "maxSpeed")
data.table::setnames(wind, "Hora Racha", "time")

#Eliminamos columnas que no nos interesan
wind$time = NULL
wind$maxSpeed= NULL

#Preparamos la fecha en formato Date para cruzar los datos
wind$date <- as.Date(wind$date)

#Ponemos la velocidad del viento en m/s

wind$Speed <- as.numeric(wind$Speed)
wind$Speed = wind$Speed * (1000/3600)<jupyter_output><empty_output><jupyter_text>Limpieza de datos Tras las modificaciones en los dataset anteriores procedemos a sustitutir los NaN como NA para que puedan ser utilizados en el estudio.Se realizará esta operacion en los dataframes:- Muestras contaminación Madrid- Precipitaciones- Viento<jupyter_code>is.nan.data.frame <- function(x)
    do.call(cbind,lapply(x,is.nan))
airMadrid <- as.data.frame(airMadrid)
airMadrid[is.nan(airMadrid)] <- NA
tail(airMadrid)
is.nan.data.frame <- function(x)
    do.call(cbind,lapply(x,is.nan))
dataFrameRainfall <- as.data.frame(dataFrameRainfall)
dataFrameRainfall[is.nan(dataFrameRainfall)] <- NA
tail(dataFrameRainfall)
is.nan.data.frame <- function(x)
    do.call(cbind,lapply(x,is.nan))
wind <- as.data.frame(wind)
wind[is.nan(wind)] <- NA
tail(wind)<jupyter_output><empty_output><jupyter_text>A continuación se eliminarán del dataframe de contaminación aquellos gases de los cuales se tenga poco datos, para ello se establecerá la línea de corte en 35%, ya que como se ha observado en diversas fuentes los gases que más afectan son NO_2,O_3,PM10,SO_2,CO y son consecuentemente, los gases que cuentan con mayor número de mediciones<jupyter_code>#1) Calculamos el porcenje de NA para descartar aquellas columnas cuyo porcentaje sea muy alto
naPerColumn <- colSums(is.na(airMadrid))

#2) Todos los gases tienen al menos un NA por, tanto, estudiamos gases a descartar, para ello
#quitamos los demás parámetros
percentageNA <- round(naPerColumn/nrow(airMadrid) *100, digits = 2)
percentageNA <- percentageNA[percentageNA > 0]
percentageNA <- percentageNA[order(percentageNA)]
percentageNA
#2)Pintamos graficos de barras para observar volumen NA en los gases contaminantes
barplot(percentageNA, main ="Porcentaje de NA",xlab = "Gas contaminante",ylab = "%",
        ylim = c(0,100), las=3, cex.names=.5, col = gray.colors(17))
#Nos quedamos con aquellos gases que tengan menos del 35% de NA
percentageNA[percentageNA < 35]
#Procedemos a quitar el resto de columnas
airMadrid$BEN = NULL;
airMadrid$EBE = NULL;
airMadrid$MXY = NULL;
airMadrid$NMHC = NULL;
airMadrid$NOx = NULL;
airMadrid$OXY = NULL;
airMadrid$PXY = NULL;
airMadrid$TCH = NULL;
airMadrid$TOL = NULL;
airMadrid$PM25 = NULL;
airMadrid$NO = NULL;
airMadrid$CH4 = NULL;

#Mostramos dataframe para observar cambios
head(airMadrid)
tail(airMadrid)
# Liberamos memoria
rm(percentageNA)
rm(naPerColumn)<jupyter_output><empty_output><jupyter_text>Agrupacion de datos Dado que las muestras de la contaminación son por hora, se agruparan las mismas realizando la media de cada uno de los gases diaria. Para ello se seguirán los siguientes pasos:- Mostrar dataframe ordenado por estación- Cambiar el formato de la fecha- Agrupar por fecha y mostrar resultado  Mostrar dataframe ordenado por estación Con esta visualización podemos comprobar que efectivamente al ordenar por estación, el siguiente campo ordenado es la fecha de la medición, por lo que mostramos le principio y el final para comprobar el rango de fechas y estaicones.<jupyter_code>#Ordenamos por estacion para comprobar que tras las 24 mediciones, pasamos al dia siguiente
airMadrid<- airMadrid[order(airMadrid$station),]
head(airMadrid, n=24)
tail(airMadrid)<jupyter_output><empty_output><jupyter_text>Cambiar el formato de la fecha Se cambia el formato de la fecha para realizar la agrupación de los datos por dicha columa.<jupyter_code>#Formato fecha
airMadrid$date <- as.Date(airMadrid$date, format = "%Y-%m-%d")
#Mostramos dataframe
head(airMadrid, n = 24)
tail(airMadrid)
#Formateamos fecha en dataframe precipitaciones
dataFrameRainfall$date <- as.Date(dataFrameRainfall$date)
rainfall<-dataFrameRainfall
rm(dataFrameRainfall)
rainfall[is.na(rainfall)] <- 0 # convertir los NAS a ceros para que pinte bien la grafica
head(rainfall, n = 10)<jupyter_output><empty_output><jupyter_text>Agrupar por fecha obteniendo la media diaria de cada gas contaminante Para los conjuntos de 24 mediciones existentes para cada día y estación, se realiza una media de las mismas, agrupando por estación y por fecha, obteniendo los siguientes resultados:<jupyter_code>madrid <- airMadrid %>%
  group_by(station,date) %>%
  summarise(CO = mean(CO, na.rm = TRUE),
            NO_2 = mean(NO_2, na.rm = TRUE),
            O_3 = mean(O_3, na.rm = TRUE),
            PM10 = mean(PM10, na.rm = TRUE),
            SO_2 = mean(SO_2, na.rm = TRUE))
head(madrid)
tail(madrid)
rm(airMadrid)<jupyter_output><empty_output><jupyter_text>Observamos que la estacion 28079099 no tiene muestras desde 2011, se debe a que fue dada de baja ese año por lo que se procede a eliminar del estudio. A continuación mostramos el final del data.frame para observar que ya no se encuentra dicha estación. Además sustituimos los NA en las mediciones de los gases por cero para poder utilizarlos.<jupyter_code>#Eliminamos estación 28079099
madrid <-  madrid[!madrid$station == 28079099,]
madrid[is.na(madrid)] <- 0
tail(madrid)<jupyter_output><empty_output><jupyter_text>Visualización de los datos a lo largo de los años A continuación se procede a visualizar los datos colocados en un timeseries tanto a nivel diario como mensual ya que apoyarán al estudio.  Gráfica con la evolución de la contaminación diaria Conviertimos el dataframe en un timeseries para visualizar la evolución diaria de la contaminacion<jupyter_code>madrid_ts_day <- rowSums(madrid[,3:7]) # agrego los gases contaminantes que corresponden a las columnas 3 a 7
madrid_ts_day <- madrid_ts_day[madrid_ts_day != 0] # elimino los ceros del timeseries porque corresponden a Nas
madrid_ts_day <- ts(madrid_ts_day, start=c(2001, 1), end=c(2018, 1), frequency=365)
rainfall_ts <- ts(rainfall$quantity, start=c(2001, 1), end=c(2018, 1), frequency=365)
plot(madrid_ts_day,main ="Contaminación diaria", xlab="Year", ylab="Contaminacion",
       col="navy")<jupyter_output><empty_output><jupyter_text>Se puede observar como la contaminación ha ido disminuyendo a lo largo de los años, así como los valores diarios maximos de la misma  Gráficos con la evolución de la contaminación y precipitación mensuales<jupyter_code># divido por año, mes y dia
madrid_month_pol <- madrid %>% separate(date, sep="-", into = c("year", "month", "day"))
pollutionMonth <- madrid %>% separate(date, sep="-", into = c("year", "month", "day"))

# sumo todos los dias del año para cada gas
madrid_month_pol <- madrid_month_pol %>%
  group_by(month) %>%
  summarise(CO = mean(CO, na.rm = TRUE),
            NO_2 = mean(NO_2, na.rm = TRUE),
            O_3 = mean(O_3, na.rm = TRUE),
            PM10 = mean(PM10, na.rm = TRUE),
            SO_2 = mean(SO_2, na.rm = TRUE))

madrid_month_pol <- madrid_month_pol[-c(18),] # quito 2018 porque solo llega hasta mayo

madrid_month_pol$totalPollution <- rowSums(madrid_month_pol[,2:6])

head(madrid_month_pol, n = 24)
# divido por año, mes y dia
rainfall_month <-rainfall %>% separate(date, sep="-", into = c("year", "month", "day"))
rainfallMonth <-rainfall %>% separate(date, sep="-", into = c("year", "month", "day"))

rainfall_month <- rainfall_month %>%
  group_by(month) %>%
  summarise(quantity = mean(quantity, na.rm = TRUE))

rainfall_month <- rainfall_month[-c(18),] # quito 2018 como en la tabla de contaminacion


head(rainfall_month, n = 12)
#Cruzamos por mes
cont_lluvia <- merge ( x = madrid_month_pol, y = rainfall_month, by = "month")
head(cont_lluvia, n = 12)
# sumo todos los dias del año para cada gas
pollutionMonth <- pollutionMonth %>%
  group_by(month,year) %>%
  summarise(CO = mean(CO, na.rm = TRUE),
            NO_2 = mean(NO_2, na.rm = TRUE),
            O_3 = mean(O_3, na.rm = TRUE),
            PM10 = mean(PM10, na.rm = TRUE),
            SO_2 = mean(SO_2, na.rm = TRUE))

pollutionMonth <- pollutionMonth[-c(18),] # quito 2018 porque solo llega hasta mayo
pollutionMonth$totalPollution <- rowSums(pollutionMonth[,3:7])
head(pollutionMonth)
rainfallMonth <- rainfallMonth %>%
  group_by(month,year) %>%
  summarise(quantity = mean(quantity, na.rm = TRUE))

rainfallMonth <- rainfallMonth[-c(18),] # quito 2018 como en la tabla de contaminacion


head(rainfallMonth, n = 12)
barplot(madrid_month_pol$totalPollution, main ="Contaminación mensual",xlab = "mes",ylab = "Contaminacion", 
        names.arg= madrid_month_pol$month, cex.names=.5)<jupyter_output><empty_output><jupyter_text>Se puede observar que en los meses mas cálidos como junio, julio, agosto y septiembre la contaminación es mayor, esto nos puede llevar a pensar que aspectos meteorológicos como la lluvia, la temperatura y el viento pueden afectar a la calidad del aire<jupyter_code>barplot(rainfall_month$quantity, main ="Precipitación mensual",xlab = "mes",ylab = "Precipitacion", 
        names.arg= rainfall_month$month, cex.names=.5)<jupyter_output><empty_output><jupyter_text>La escasez de lluvias en junio, julio, agosto y septiembre podria ser una de las causas del incremento de la contaminacion en dichos meses  Gráficas con la evolución de la contaminación y precipitaciones anuales<jupyter_code># divido por año, mes y dia
madrid_year_pol <- madrid %>% separate(date, sep="-", into = c("year", "month", "day"))

# sumo todos los dias del año para cada gas
madrid_year_pol <- madrid_year_pol %>%
  group_by(year) %>%
  summarise(CO = mean(CO, na.rm = TRUE),
            NO_2 = mean(NO_2, na.rm = TRUE),
            O_3 = mean(O_3, na.rm = TRUE),
            PM10 = mean(PM10, na.rm = TRUE),
            SO_2 = mean(SO_2, na.rm = TRUE))

madrid_year_pol <- madrid_year_pol[-c(18),] # quito 2018 porque solo llega hasta mayo

madrid_year_pol$totalPollution <- rowSums(madrid_year_pol[,2:6]) # agrego los gases contaminantes que corresponden a las columnas 2 a 6
head(madrid_year_pol, n = 18)

# divido por año, mes y dia
rainfall_year <-rainfall %>% separate(date, sep="-", into = c("year", "month", "day"))

rainfall_year <- rainfall_year %>%
  group_by(year) %>%
  summarise(quantity = mean(quantity, na.rm = TRUE))

rainfall_year <- rainfall_year[-c(18),] # quito 2018 como en la tabla de contaminacion

head(rainfall_year, n = 18)
barplot(madrid_year_pol$totalPollution, main ="Contaminación anual",xlab = "año",ylab = "Contaminación", 
        names.arg= madrid_year_pol$year, cex.names=.5)

barplot(rainfall_year$quantity, main ="Precipitación anual",xlab = "año",ylab = "Precipitación", 
        names.arg= rainfall_year$year, cex.names=.5)<jupyter_output><empty_output><jupyter_text>No podemos sacar conclusiones claras cuando evaluamos la contaminación y precipitaciones anuales, lo que si se observa es un descenso notable de la contaminacion a partir de 2009, y entre los años 2015 a 2017 el aumento de las precipitaciones en el año 2016 parece provocar un descenso en el nivel de contaminacion  Correlación de variables<jupyter_code>madrid_cor <- madrid
head(madrid_cor)
#Cruzamos por fecha y ordenamos por estacion
madrid_cor <- merge ( x = madrid_cor, y = rainfall, by = "date")
madrid_cor <- madrid_cor[order(madrid_cor$station),]
head(madrid_cor)
#Eliminamos columnas
madrid_cor$date = NULL
madrid_cor$station = NULL
head(madrid_cor)
# Añadimos columnas gases contaminantes dia siguiente
madrid_cor$CO_NextDay <- c(madrid_cor$CO[-1], NA)
madrid_cor$NO2_NextDay <- c(madrid_cor$NO_2[-1], NA)
madrid_cor$O_3NextDay <- c(madrid_cor$O_3[-1], NA)
madrid_cor$PM10_NextDay <- c(madrid_cor$PM10[-1], NA)
madrid_cor$SO_2NextDay <- c(madrid_cor$SO_2[-1], NA)

# Quitamos columnas gases contaminantes dia actual
madrid_cor$CO = NULL
madrid_cor$NO_2 = NULL
madrid_cor$O_3 = NULL
madrid_cor$PM10 = NULL
madrid_cor$SO_2 = NULL

head(madrid_cor)
# Calculamos el maximo entre todos los gases del dia siguiente y la polucion total del dia siguiente
madrid_cor$maxNextDay <- apply(madrid_cor[,2:6], 1, FUN=max)
head(madrid_cor)
madrid_cor$totalPolNextDay <- rowSums(madrid_cor[,2:6])
head(madrid_cor)<jupyter_output><empty_output><jupyter_text>Tabla de correlaciones<jupyter_code>#Tabla de correlaciones
res <- cor(madrid_cor, use = "complete.obs")
round(res, 2)
#Mostramos la siguiente gráfica en fucnión de la correlación entre variables.
corrplot(res)
title("Correlacion entre la lluvia y la contaminacion de cada gas")<jupyter_output><empty_output><jupyter_text>Las correlaciones más altas que hemos obtenido con la lluvia del dia anterior han sido con la polución total y con el gas más contaminante del día siguiente, –0,12 y – 0,10 respectivamente.  Preparacion para el modeloSe van a crear tres dataframes, para la aplicación posterior del modelo a cada uno de ellos, estos dataframes serán:- version01: Dataframe que contiene los datos discretizados siguiendo índices obtenidos en fuentes que se indicarán a posteriori.- version02: Dataframe que contiene los datos discretizados por cuartiles, mediciones mensuales- version03: Dataframe que contiene los datos discretizados por cuartiles en tres variables categoricas, mediciones diarias    - Una situada hasta el primer cuartil (_data_ <= Q1)    - Otro siatuada entre el primer y tercer cuartil (Q1< _data_ <= Q3)      - La última para datos superiores al tercer cuartil (_data_> Q3)  Prepación para el modelo version01 Tras los pasos anteriores se procederá a realizar la discretización de las variables dependientes que formarán parte del modelo para preparar el dataframe para el modelo.A lo largo de esta fase se utilizarán los siguientes términos:- AQ (Air Index): Representa mediciones numéricas- AIQ (Air Index Quality): Representa mediciones ya discretizadas  Discretizamos la calidad del airePara ello nos guiraremos de la siguiente tabla, de cada fila se escogerá el mayor de todos los contaminantes y se comparara con la tabla adjunta: ![airIndex.png](attachment:airIndex.png)Fuente: https://www.epa.vic.gov.au/your-environment/air/air-pollution/air-quality-index/calculating-a-station-air-quality-index Para ello cogemos el maximo valor de cada fila entre los gases y lo añadimos como una columna nueva al data frame<jupyter_code>#Obtenemos en una columna el máximo de los gases contaminantes
auxData <- madrid
auxData$date=NULL
auxData$station = NULL
auxData[, "AQ"] <- apply(auxData ,1, max)
auxData$CO=NULL
auxData$NO_2 = NULL
auxData$O_3 = NULL
auxData$PM10 = NULL
auxData$SO_2 = NULL
madrid <- cbind(madrid,auxData)
head(madrid)<jupyter_output><empty_output><jupyter_text>A continuación discretizamos la columna anterior, creando una nueva columna con los siguientes valores:    Very good(VG)      Good (G)      Fair(F)      Poor (P)      Very Poor (VP)<jupyter_code>#Preparamos para discretizar
madrid[,"AIQ"] <- madrid$AQ
head(madrid)
#Discretizamos
madrid <- madrid %>%
        mutate(AIQ = case_when (
        
        AIQ>=0 & AIQ <=33       ~"VERY GOOD",
        AIQ >=34 & AIQ <=66     ~"GOOD",
        AIQ >=37 & AIQ <=99     ~"FAIR",
        AIQ >=100 & AIQ <=149   ~"POOR",
        AIQ >=150               ~"VERY POOR",
        TRUE                    ~"NA",
        ))

    
head(madrid)
#Preparamos para discretizar el siguiente dia
madrid[,"AQNextDay"] <- NA
madrid[,"AIQNextDay"] <- NA
head(madrid)
#Discretizamos el día siguiente
madrid$AQNextDay <- c(madrid$AQ[-1],NA)
madrid$AIQNextDay <- c(madrid$AIQ[-1],NA)

madrid02<- madrid #lo guardamos para la siguiente version
head(madrid)<jupyter_output><empty_output><jupyter_text>Discretizamos el volumen de precipitaciones Para ello nos guiraremos de la siguiente tabla, de cada fila se escogerá la cantidad y se comparara con la tabla adjunta: ![rainfall.png](attachment:rainfall.png)<jupyter_code>#Preparamos para discretizar
rainfall[,"rainfallVolume"] <- rainfall$quantity
head(rainfall)
#Discretizamos lluvia
rainfall <- rainfall %>%
        mutate(rainfallVolume = case_when (
             
        rainfall$rainfallVolume <= 2                                      ~ "WEAK",
        rainfall$rainfallVolume > 2 & rainfall$rainfallVolume <= 15       ~ "MODERATE",
        rainfall$rainfallVolume > 15 & rainfall$rainfallVolume <= 30      ~ "STRONG",
        rainfall$rainfallVolume > 30 & rainfall$rainfallVolume <=60       ~ "VERY STRONG",
        rainfall$rainfallVolume > 60                                      ~ "TORRENTIAL RAIN",
        TRUE                                                              ~"NA",
        ))
head(rainfall)
tail(rainfall)<jupyter_output><empty_output><jupyter_text>Discretizamos la fuerza del viento Para ello nos guiraremos de la siguiente tabla, de cada fila se escogerá el viento y se comparara con la tabla adjunta: ![wind.png](attachment:wind.png) A continuación discretizamos la columna anterior, creando una nueva columna con los siguientes valores:    LIGTH WIND (LW)     ADVISORY FORCE (AF)    GALE FORCE (GF)    STORM FORCE (SF)    HURRICANE FORCE (HF)<jupyter_code>#Discretizamos el viento
wind[,"WindIndexQuality"] <- wind$Speed
wind <- wind %>%
        mutate(WindIndexQuality = case_when (
        
        WindIndexQuality >=0.0 & WindIndexQuality  <=10.7       ~"LIGHT WIND",
        WindIndexQuality >=10.8 & WindIndexQuality <=17.1       ~"ADVISORY FORCE",
        WindIndexQuality >=17.2 & WindIndexQuality <=24.4       ~"GALE FORCE",
        WindIndexQuality >=24.5 & WindIndexQuality <=32.6       ~"STORM FORCE",
        WindIndexQuality >=32.7                                 ~"HURRICANE FORCE",
        TRUE                                                    ~"NA",
        ))
head(wind)
tail(wind)<jupyter_output><empty_output><jupyter_text>Cruzamos los datos Finalmente como último paso previo a la aplicación del modelo se unieron todos los datos en un único dataframe<jupyter_code># Cruzamos los datos de las precipitaciones y la contaminacion por fecha 
version01 <- merge ( x = madrid, y = rainfall, by = "date")
version01 <- version01[order(version01$station),]
head(version01)
# Cruzamos los datos con el viento
version01 <- merge ( x = version01, y = wind, by = "date",)
version01 <- version01[order(version01$station),]
head(version01)<jupyter_output><empty_output><jupyter_text>Preparación para el modelo version02 Para esta parte del estudio, se han obtenido los cuartiles correspondientes a la precipitación y contaminación anual para discretizar el mismo<jupyter_code>#Preparamos la versión 02
head(pollutionMonth)<jupyter_output><empty_output><jupyter_text>Discretizacion con cuartiles<jupyter_code>#Obtenemos cuartiles
summary(pollutionMonth$totalPollution)<jupyter_output><empty_output><jupyter_text>Discretizamos cuartiles<jupyter_code>madrid_cuart <- pollutionMonth
#Discretizamos cuartiles
madrid_cuart$calidadAire <- madrid_cuart$totalPollution

madrid_cuart <- madrid_cuart %>%
        mutate(calidadAire = case_when (
        calidadAire <  81.68                               ~"ALTA",
        calidadAire >= 81.68  & calidadAire < 110.60       ~"MEDIA",
        calidadAire >= 110.60 & calidadAire < 136.86       ~"BAJA",
        calidadAire >= 136.86                            ~"MUY BAJA",
        TRUE                                               ~"NA",
        ))
madrid_cuart$date <- as.yearmon(paste(madrid_cuart$year, madrid_cuart$month), "%Y %m")
madrid_cuart$month=NULL
madrid_cuart$year=NULL
head(madrid_cuart)
#Obtenemos cuartiles
summary(rainfallMonth$quantity)
#Discretizamos lluvia por cuartiles
rainfall_cuart <- rainfallMonth
rainfall_cuart$lluvia <- rainfall_cuart$quantity

rainfall_cuart <- rainfall_cuart %>%
        mutate(lluvia = case_when (
        lluvia <= 0.2034                       ~"DEBIL",
        lluvia > 0.2034 & lluvia <=1.5935      ~"MODERADA",
        lluvia > 1.5935                        ~"FUERTE",
        TRUE                                   ~"NA",
        ))

rainfall_cuart$date <- as.yearmon(paste(rainfall_cuart$year, rainfall_cuart$month), "%Y %m")
rainfall_cuart$month=NULL
rainfall_cuart$year=NULL
head(rainfall_cuart, n = 12)
# Cruzamos los datos de las precipitaciones y la contaminacion por fecha 
version02 <- merge ( x = madrid_cuart, y = rainfall_cuart,by = "date")
version02<- version02[order(version02$date),]
head(version02)
cor(x = version02$totalPollution, y = version02$quantity, use="complete.obs")<jupyter_output><empty_output><jupyter_text>Preparación para el modelo version03<jupyter_code>#Preparamos la versión 03
madrid02$AIQ = NULL
madrid02$AIQNextDay = NULL
head(madrid02)<jupyter_output><empty_output><jupyter_text>Obtenemos información sobre la media y los cuartiles de las mediciones sobre la contaminación del aire<jupyter_code>#Obtenemos cuartiles
summary(madrid02$AQ)<jupyter_output><empty_output><jupyter_text>Discretizamos la calidad del aire<jupyter_code>#Discretizamos por cuartiles
madrid02[,"AIQ"] <- madrid02$AQ
madrid02[,"AIQNextDay"] <- madrid02$AQNextDay

madrid02 <- madrid02 %>%
        mutate(AIQ = case_when (
        
        AQ <= 46.03                      ~"GOOD",
        AQ > 46.03 &  AQ <= 72.92      ~"MEDIUM",
        AQ > 72.92                       ~"BAD",
        TRUE                    ~"NA",
        ))

madrid02 <- madrid02 %>%
        mutate(AIQNextDay = case_when (
        
        AQNextDay <= 46.03                      ~"GOOD",
        AQNextDay > 46.03 &  AQNextDay <= 72.92      ~"MEDIUM",
        AQNextDay > 72.92                       ~"BAD",
        TRUE                    ~"NA",
        ))
head(madrid02)<jupyter_output><empty_output><jupyter_text>Obtenemos información sobre la media y los cuartiles de las mediciones de las precipitaciones<jupyter_code>#Obtenemos cuartiles
summary(rainfall$quantity)<jupyter_output><empty_output><jupyter_text>Discretizamos el volumen de precipitaciones<jupyter_code>#Discretizamos precipitaciones
rainfall <- rainfall %>%
        mutate(rainfallVolume = case_when (
             
        rainfall$quantity <= 0                                      ~ "LOW",
        rainfall$quantity > 0       ~ "HIGH",
        TRUE                                                              ~"NA",
        ))
head(rainfall)<jupyter_output><empty_output><jupyter_text>Obtenemos información sobre la media y los cuartiles de las mediciones de la velocidadd el viento<jupyter_code>#Obtenemos por cuartiles
summary(wind$Speed)<jupyter_output><empty_output><jupyter_text>Discretizamos la velocidad del viento<jupyter_code>#Discretizamos velocidad del viento
wind <- wind %>%
        mutate(WindIndexQuality = case_when (
        
        Speed < 0.3                     ~"LOW",
        Speed >= 0.3 & Speed <= 2.2      ~"MEDIUM",
        Speed >2.2                      ~"HIGH",
        TRUE                        ~"NA",
        ))
head(wind)
# Cruzamos los datos de las precipitaciones y la contaminacion por fecha 
version03 <- merge ( x = madrid02, y = rainfall, by = "date")
version03 <- version03[order(version03$station),]
head(version03)
# Cruzamos los datos con el viento
version03 <- merge ( x = version03, y = wind, by = "date",)
version03 <- version03[order(version03$station),]
head(version03)
cor(x = version03$AQNextDay, y = version03$quantity, use="complete.obs")<jupyter_output><empty_output><jupyter_text>4º FASE: MODEL Se entrenan y aplican distintos modelos para el procesamiento de los datos, seleccionando aquel cuya tecnología y/o parámetros produce unos resultados aparentemente más adecuados.En este estudio primero se aplicó un modelo de regresión con el obetivo de obtener la contaminación del día siguiente, el resultado no fue el esperado, por lo que se aplicó durante el resto del estudio el modelo árbol de decisiones  Árbol de decisiones version01 Inicializamos el gestor de memoria<jupyter_code>#Gestor de memoria
gc()
# Establece una semilla fija para garantizar
# la reproducibilidad del ejemplo
head(version01)
tail(version01)
str(version01)
set.seed(1234)
#Preparamos validación simple
separate <- as.integer(nrow(version01)*0.7)
trainData <- version01[1:separate,]
testData <- version01[separate:nrow(version01),]
head(trainData)
head(testData)
myFormula <- AIQNextDay ~ rainfallVolume + WindIndexQuality + AIQ
madrid_rpart <- rpart(myFormula,data=trainData, control = rpart.control(minsplit = 10))<jupyter_output><empty_output><jupyter_text>Comprobación de los datos de entrenamiento<jupyter_code>attributes(madrid_rpart)
summary(madrid_rpart)
fancyRpartPlot(madrid_rpart, caption = NULL)<jupyter_output><empty_output><jupyter_text>Las variables independientes son la calidad del aire(AIQ), las precipitaciones y la velocidad del viento del dia anterior. La unica variable independiente que sale representada en el arbol de decisiones es la calidad del aire del dia anterior(AIQ) al ser la que adquiere mayor importancia,aunque se tienen en cuenta el resto de variables independientes.La variable dependiente es la calidad del aire del dia siguiente.Si la calidad del aire del dia anterior(AIQ) es asumible(fair), pobre o muy pobre teniendo en cuenta el valor del resto de variables independientes, la calidad del aire al dia siguiente será asumible en un 39% de los casos, en caso contrario la calidad del aire al dia siguiente será buena en un 61% de los casos.  Predicción con los datos de test Se añadirá una nueva columna, en esta columna se encontrará el valor predicho por el árbol de decisión<jupyter_code>testData$AIQExpected <- predict(madrid_rpart, newdata = testData, type = "class")
head(testData)
predictTable<- table(testData$AIQExpected,testData$AIQNextDay)
predictTable
accuracy <- sum(diag(predictTable))/sum(predictTable)
accuracy<jupyter_output><empty_output><jupyter_text>Observamos que obtenemos un nivel de confianza del 61 %, este es muy bajo por lo que se probará a utilizar otro modelo predictivo<jupyter_code>rm(madrid_rpart)
rm(testData)
rm(trainData)
rm(predictTable)<jupyter_output><empty_output><jupyter_text>Árbol de decisiones version02 Inicializamos el gestor de memoria<jupyter_code>#Gestor de memoria
gc()
# Establece una semilla fija para garantizar
# la reproducibilidad del ejemplo
head(version02)
tail(version02)
str(version02)
set.seed(1234)
separate <- as.integer(nrow(version02)*0.7)
trainData <- version02[1:separate,]
testData <- version02[separate:nrow(version02),]
head(trainData)
head(testData)
myFormula <- myFormula <- calidadAire ~ lluvia
madrid_rpart <- rpart(myFormula,data=trainData, control = rpart.control(minsplit = 10))<jupyter_output><empty_output><jupyter_text>Comprobación de los datos de entrenamiento<jupyter_code>attributes(madrid_rpart)
summary(madrid_rpart)
fancyRpartPlot(madrid_rpart, caption = NULL)<jupyter_output><empty_output><jupyter_text>La variable independiente son las precipitaciones del dia anterior.La variable dependiente es la calidad del aire al dia siguiente.Si la lluvia del dia anterior es fuerte, la calidad del aire al dia siguiente sera baja en el 26% de los casos, en caso contrario la calidad del aire sera muy baja en el 74% de los casos.  Predicción con los datos de test Se añadirá una nueva columna, en esta columna se encontrará el valor predicho por el árbol de decisión<jupyter_code>testData$AIQExpected <- predict(madrid_rpart, newdata = testData, type = "class")
head(testData)
predictTable<- table(testData$calidadAire,testData$AIQExpected)
predictTable
accuracy <- sum(diag(predictTable))/sum(predictTable)
accuracy
rm(madrid_rpart)
rm(testData)
rm(trainData)
rm(predictTable)<jupyter_output><empty_output><jupyter_text>Árbol de decisiones version03 Inicializamos el gestor de memoria<jupyter_code>#Gestor de memoria
gc()
# Establece una semilla fija para garantizar
# la reproducibilidad del ejemplo

head(version03)
tail(version03)
str(version03)
set.seed(1234)
#Preparamos validación simple
separate <- as.integer(nrow(version03)*0.7)
trainData <- version03[1:separate,]
testData <- version03[separate:nrow(version03),]
head(trainData)
head(testData)
myFormula <- AIQNextDay ~ AIQ + rainfallVolume + WindIndexQuality
madrid_rpart <- rpart(myFormula,data=trainData, control = rpart.control(minsplit = 10))<jupyter_output><empty_output><jupyter_text>Comprobación de los datos de entrenamiento<jupyter_code>attributes(madrid_rpart)
summary(madrid_rpart)
fancyRpartPlot(madrid_rpart, caption = NULL)<jupyter_output><empty_output><jupyter_text>Cada caja de nodo muestra la clasificación, la probabilidad de cada clase en ese nodo (es decir, la probabilidad de la clase condicionada en el nodo) y el porcentaje de observaciones utilizadas en ese nodo.  Las variables independientes son la calidad del aire(AIQ), las precipitaciones y la velocidad del viento del dia anterior. La variable dependiente es la calidad del aire del dia siguiente.Si la calidad del aire del dia anterior(AIQ) es mala teniendo en cuenta el valor del resto de variables independientes, la calidad del aire al dia siguiente será mala en un 28% de los casos, en caso contrario la calidad del aire al dia siguiente será media en un 72% de los casos. Si ademas de no ser mala la calidad del aire del dia anterior es buena, la calidad del aire del dia siguiente sera buena en un 20% de los casos y media en un 52% de los casos.  Predicción con los datos de test Se añadirá una nueva columna, en esta columna se encontrará el valor predicho por el árbol de decisión<jupyter_code>testData$AIQExpected <- predict(madrid_rpart, newdata = testData, type = "class")
head(testData)
predictTable<- table(testData$AIQNextDay,testData$AIQExpected)
predictTable
accuracy <- sum(diag(predictTable))/sum(predictTable)
accuracy<jupyter_output><empty_output><jupyter_text>Liberamos memoria<jupyter_code>rm (list=ls())<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1002635472535602, 'fraction_numerical': 0.03141760680599223, 'mean_word_length': 4.929677861549006, 'pattern_counts': {'":': 0, '<': 406, '<?xml version=': 0, '>': 249, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '47', 'dir': 'jupyter-structured-clean-dedup', 'id': '262051', 'license': 'MIT', 'n_tokens_mistral': 16068, 'n_tokens_neox': 14629, 'n_words': 4040, 'path': 'Estudio.ipynb', 'repo_name': 'manumonforte/DataMining-BigDataParadigm', 'size': '43256'}