PILDORA R-1. CALIDAD DEL AIRE CON R.

Hace tiempo que quería realizar una serie artículos relacionados con la minería de datos de calidad del aire con R, en formato de ejercicio práctico, que sirviese de referencia a nuevos usuarios y mentes inquietas para iniciarse de forma practica en este lenguaje y en los temas de calidad del aire. Aquí empiezan esas Píldoras de información que complementan al ya antiguo y algo obsoleto Manual R y Openair aplicados a la Calidad del Aire que os podrán resultar de interés. Espero que os gusten.


IMPORTANDO ARCHIVOS FINN-A DE CALIDAD DEL AIRE.

Es más que habitual que en materia de Calidad del Aire en España los datos se pongan a disposición del "público en general" a través del mismo formato normalizado de intercambio de datos que utilizan los organismos competentes para compartirlos entre ellos, el conocido como formato FINN-A que no es más que un archivo de texto plano organizado en horizontal en base a los contaminantes y las fechas.

Esta organización horizontal con varios datos por fila no es compatible con el estudio o manejo de datos por dataframe de R o con Openair, por lo que a la hora de tener estos datos tendremos que pensar en alguna forma de convertirlos o pasarlos a un formato R compatible. Aquí os propongo una de las posibles formas de hacerlo.

Los datos para el ejercicio se pueden encontrar y descargar muy fácilmente desde la plataforma siguiente de Medio Ambiente de la Generalitat, en el apartado de datos abiertos, buscando siempre en "veure les dades" y filtrando por los que necesitemos, que si no el archivo será gigante.



Una vez tenemos el archivo de datos, que no es más que un archivo de texto plano, con los datos de la estación, organizado por días y parámetros, poniendo en cada fila los datos horarios (es decir 24 datos con sus 24 flags en cada fila) procedemos a su importación a nuestro proyecto R/ RStudio para empezar a manejarlo.

Tan sencillo como ir a "Import Dataset" desplegar el menu y decirle a RStudio que nos vamos a traer un archivo de texto plano "From Text (base)..."




Ya tenemos nuestro archivo en R, podemos echarle un ojo para ver qué es lo que nos interesa, y ¿qué campos tienen el qué?.... también sería interesante que comprobásemos que la estructura del archivo es la que nos interesa de cara a que cada campo sea como nosotros queremos. Hay veces que si alguna columna trae algún dato extraño podemos encontrarnos con que R la reconozca como "character", "numeric", "factor", etc. en cuyo caso tendríamos que forzar su conversión.

   #Compruebo la estructura de los datos del archivo de calidad del aire
   str(Sabadell)

El archivo hay que filtrarlo para quedarnos con lo que nos interesa, en este punto, yo me voy a quedar con un archivo con menos columnas seleccionando la columna del contaminante "CONTAMINANT", que es la 14, la columna de la fecha "DATA", que es la columna 20, y las columnas correspondientes a los datos, que son las columnas 21 a 68.

   #Seleccciono las columnas que me interesa traerme al archivo nuevo.
   Sabadell2 <- Sabadell [ , c(14,20:68) ]

Si quisiera también podría meter condiciones de filtrado para las filas, poniendo la fila o filas que quisiese, o la condición que vea más interesante. Le puedo decir, por ejemplo que me de solo las filas correspondientes al año 2.019.

   #Selecciono las filas correspondientes al año 2.019 y las columnas que me interesan
   Sabadell2 <- Sabadell [ Sabadell$ANY == 2019 , c(14,20:68) ]

Los datos en calidad del aire siempre vienen de dos en dos. Primero el dato numérico, y luego el "flag" o letra que lo califica. Dicho flag establece si el dato es válido "V", inválido "N", si está corregido "O", si pertenece a una calibración "C", etc. La situación puede ser muy variopinta.

En nuestro caso, lo que nos interesa es coger el dato validado "V" que es el que sirve para establecer la calidad del aire en la zona, así que para dejar sólo estos datos en nuestro conjunto utilizaremos una función muy interesante de la librería "dplyr" de R que se llama "mutate".

A la función mutate sólo  hay que decirle el nombre del data frame y el cambio que queremos hacerle a la columna correspondiente, para lo cual es habitual usar un condicional "ifelse". En nuestro caso queremos decirle que mute el dataframe Sabadell2, cogiendo las columnas de las horas "H01","H02".... y conservando el número si "if" la columna de validados es "V" o "else" descartando el número y poniendo "NA" si la columna es otro "flag". El código sería el siguiente:

   #Llamo a la librería que tiene el comando que voy a utilizar.
   library(dplyr)

   #Propongo la mutación de filas que quiero hacer.
   Sabadell2<-mutate(Sabadell2, H01=ifelse(V01=="V",H01,NA), H02=ifelse(V02=="V",H02,NA),
   H03=ifelse(V03=="V",H03,NA), H04=ifelse(V04=="V",H04,NA), H05=ifelse(V05=="V",H05,NA), 
   H06=ifelse(V06=="V",H06,NA), H07=ifelse(V07=="V",H07,NA), H08=ifelse(V08=="V",H08,NA), 
   H09=ifelse(V09=="V",H09,NA), H10=ifelse(V10=="V",H10,NA), H11=ifelse(V11=="V",H11,NA), 
   H12=ifelse(V12=="V",H12,NA), H13=ifelse(V13=="V",H13,NA), H14=ifelse(V14=="V",H14,NA), 
   H15=ifelse(V15=="V",H15,NA), H16=ifelse(V16=="V",H16,NA), H17=ifelse(V17=="V",H17,NA), 
   H18=ifelse(V18=="V",H18,NA), H19=ifelse(V19=="V",H19,NA), H20=ifelse(V20=="V",H20,NA), 
   H21=ifelse(V21=="V",H21,NA), H22=ifelse(V22=="V",H22,NA), H23=ifelse(V23=="V",H23,NA), 
   H24=ifelse(V24=="V",H24,NA))

Si lo hemos hecho bien, nuestro dataframe tiene que haber sustituido todos los valores con flag "N" o similar en valores "NA", tal y como podemos observar. Esto nos permite contar con los datos sólo de validados.



Ahora el problema es conseguir transformar un archivo de valores en horizontal en un archivo de valores hábil para utilizar en R con los datos indexados por la fecha y los valores dispuestos por columnas en función del contaminante. Aquí es donde está nuestro principal reto de programación en R. Yo os voy a dejar una de las formas en la que lo podemos hacer, quizás la más intuitiva.

En primer lugar usamos el comando "subset" para seleccionar las columnas de fecha y datos de nuestro archivo en base al contaminante que queramos extraer. El comando "subset" de R nos permitirá llevar a cabo este trabajo de una forma sencilla, y sólo nos requiere que le digamos el data frame, la condición y la columna a extraer. Haremos tantos archivos como datos horarios haya que extraer para cada contaminante. El código sería el siguiente:

   #Voy a extraer los datos por contaminantes y hora para luego juntarlos.
   NOX01<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H01"))
   NOX02<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H02"))
   NOX03<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H03"))
   NOX04<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H04"))
   NOX05<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H05"))
   NOX06<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H06"))
   NOX07<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H07"))
   NOX08<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H08"))
   NOX09<-subset(Sabadell2, CONTAMINANT=="NOX", select=c("DATA","H09"))
   (...)


El siguiente punto es manejar la fecha de cada uno de los archivos que hemos creado para conseguir que se interprete como una fecha (ahora mismo es sólo un factor), y que además de la fecha se corresponda con la hora. Para ello usaremos el paquete de R "lubridate", un paquete muy útil que nos permite manejar rapidamente las horas.

El primer paso, tras llamar al paquete, es pedirle que reconozca la hora en función del formato que nos viene dado. En nuestro caso, el archivo viene como "día/mes/año", por lo que la función de lubridate a utilizar sera "dmy" y se aplicaría sobre la columna "DATA" de cada archivo.

De esta forma ya tendíramos la columna "DATA" reconocida por R como fecha. Para poder sumarle una hora, lo que tenemos que hacer es convertirla en un número con el comando "as.numeric", indicándole la columna "DATA".

Al convertirla en número R lo interpreta como el número de días desde el 1 de enero de 1970, por lo que habrá que pasar el número a segundos, para conseguir que R lo interprete como fecha y hora en el siguiente comando. Multiplicaremos por tanto el número por 86.400 y le sumaremos las horas que correspondan a cada archivo convertidas en segundos: 3600, 7200, ....

Una vez hecho esto, usaremos el comando "as.datetime" para convertir el número de nuevo en una fecha y hora, y unifico los nombres de las columnas con el comando "names" para que luego no me den problemas en la unificación. El código quedaría como sigue, para cada uno de los archivos de NOX que tenemos.

   #Llamo al paquete para el manejo de fechas
   library(lubridate)
   #convierto el factor DATA en fecha
   NOX01$DATA<-dmy(NOX01$DATA)
   #transformo la fecha en un numérico
   NOX01$DATA<-as.numeric(NOX01$DATA)
   #Convierto el fichero a segundos y le sumo una hora en segundos
   NOX01$DATA<-(NOX01$DATA*86400)+3600
   #Convierto la columna numérica en una fecha con hora.
   NOX01$DATA<-as_datetime(NOX01$DATA)

   #Estandarizo los nombres del archivo para unificar los datos más adelante.
   names(NOX01)<-c("date","nox")

Si lo hemos hecho bien, el resultado debería ser el siguiente para el fichero de NOX01, y exactamente igual para los sucesivos.



Una vez tenemos el archivo de datos de NOX de las distintas horas, es tan sencillo como unir los distintos dataframes creados por filas, con el comando "rbind" y borrar los archivos intermedios que hemos creado, con el comando "rm", para limpiar el espacio de trabajo de RStudio.

   #Unifico todos los data frames individuales creados para Sabadell
   NOX<-rbind(NOX01,NOX02, NOX03,NOX04,NOX05,NOX06,NOX07,NOX08,NOX09,NOX10,

   NOX11,NOX12,NOX13,NOX14,NOX15,NOX16,NOX17,NOX18,NOX19,NOX20,NOX21,
   NOX22,NOX23,NOX24)
 
   #Borro los archivos creados intermedios
   rm(NOX01,NOX02, NOX03,NOX04,NOX05,NOX06,NOX07,NOX08,NOX09,NOX10,

   NOX11,NOX12,NOX13,NOX14,NOX15,NOX16,NOX17,NOX18,NOX19,NOX20,
   NOX21,NOX22,NOX23,NOX24)

Repitiendo el mismo código para el resto de parámetros que queremos extraer (tan sencillo como copiar el código sustituyento el parámetro "NOX" por el que corresponda "PM10", "O3"...) tendremos un archivo por cada contaminante de la estación de Sabadell. Aquí ya sólo nos quedaría unir los distintos data frame para conseguir tener un solo data frame de la estación de Sabadell ya preparado para trabajar bajo R con total fluidez.

Para hacerlo podemos usar el comando "merge", que nos permite unir data frames en función de una columna común de datos, que en este caso sería "date". La función sólo permite unir dos data frames, que se nombran al principio, poniendo luego la columna común.

El último parámetro que le marcamos al comando le indica qué columna conservará todos los datos. En este caso, nos interesa conservar todos los datos de todos los data frame, por lo que R rellenará con "NA" aquellas columnas que no coincidan en fecha por uno u otro lado.

El código será el siguiente:

   #Utilizo el comando merge para juntar los dataframe de contaminantes por fecha

   #Junto los datos de PM10 y O3 por fecha conservando las fechas de ambos dos
   sabadell01<-merge(PM10,O3, by="date", all=TRUE)


   #Junto el archivo anterior con el de NOx conservando todas las fechas
   sabadell02<-merge(sabadell01,NOX, by="date", all=TRUE)


   #Junto el anterior con el de NO2 conservando todas las fechas
   Sabadellaire<-merge(sabadell02,NO2, by="date", all=TRUE)


   #Borro de mi espacio de trabajo todo lo que me sobra
   rm(sabadell01,sabadell02,Sabadell, Sabadell2, NO2,NOX,O3,PM10)


Si lo hemos hecho bien, el resultado final de nuestro código tiene que ser el siguiente:



Ya tenemos un data frame de datos puros de calidad del aire, en este caso para la estación de Sabadell, que nos permitirán hacer estudios posteirores sobre calidad del aire que podremos ver con posterioridad.

Espero que esta píldora os haya resultado de interés y utilidad. 




Comentarios

  1. Gracias por el gran aporte de tu documento. Me ha surgido una duda durante el intento de seguir los pasos con un fichero de datos de calidad del aire. Me refiero a la parte de "#Convierto el fichero a segundos y le sumo una hora en segundos", NOX01$DATA<-(NOX01$DATA*86400)+3600. Supongo que para el resto de fichero correspondientes a las respectivas horas el último termino lo multiplicas por 2x3600 en el caso de 2 hr. NOX02 y así sucesivamente. Y en el caso de NOX24 lo multiplicas por 3540 para que no pase de 24 hr pues sería otro día. Gracias

    ResponderEliminar

Publicar un comentario

Tus comentarios son siempre bienvenidos. Agrega aquí lo que desees en relación al artículo publicado.

Entradas populares de este blog

TRÁFICO Y CALIDAD DEL AIRE EN MADRID