PILDORA R-2. CALIDAD DEL AIRE CON R.

CAPTURANDO DATOS METEOROLÓGICOS PARA NUESTRO CONJUNTO.

Muchas estaciones de calidad del aire están dirigidas en exclusiva a registrar los datos de calidad del aire para dar cumplimiento a la legalidad vigente y comprobar si se superan los valores límite y objetivos de calidad del aire, disponiendo de datos de los distintos parámetros de contaminación, sin parámetros adicionales de meteorología.

Tener datos de calidad del aire sin disponer a la vez de datos de meteorología es como tener una ventana en  una casa y no disponer de manija para abrirla. Puedes ver si fuera  hace buen o mal tiempo, pero es imposible abrir la puerta y conocer bien qué tiempo hace y qué te puedes poner para salir.

 

Esto, no obstante, se puede solventar muchas veces bajando e integrando datos de meterología de estaciones próximas, en ocasiones de estaciones que están justo al lado de la propia estación, pero que reportan datos al organismo meteorologico correspondiente por otro lado, sin integrarse a la misma.

Para ilustrar este ejercicio seguiremos con nuestra estación de Sabadell de la Píldora 1, donde aprendimos a extraer datos en R de un archivo estandar de reporte de datos en texto plano. Ahora nos iremos a buscar un archivo de meteorología de esta misma localidad, pero al repositorio de datos abiertos de meteorología, muy similar de funcionamiento al que vimos para la Calidad del Aire.

Una vez en la página nos iremos al apartado de "veure les dades", activaremos la opción "Filtre" y nos lanzaremos a filtrar la estación de "Sabadell", que es el "CODI_ESTACIO" con las letras "XF". Esta información la podríamos consultar en la página de metadatos de la propia web de datos meteorológicos abiertos.

Una vez filtrado el archivo, y exportado a un csv a nuestro ordenador, nos podemos poner manos a la obra con R. Echarle un repaso a la Píldora 1 para recordar cómo se importa el archivo de datos csv a nuestro sistema,  a parti de ahí el archivo con el que empezaríamos a trabajar en R sería algo de este aspecto general de partida:

Veremos que en el caso de los archivos de meteorología los datos suelen estar organizados de forma distinta a la calidad del aire. En este caso nos viene un dato por variable y fecha en cada línea, y es posible que nos encontremos con que la agregación de los datos sea incluso diferente (datos diezminutales, semihorarios, etc). Más diversión para nuestro trabajo de integración.

Lo primero que vamos a hacer es asegurarnos que los datos que tenemos son validados, convirtiendo en "NA" (datos nulos) todos aquellos que no se correspondan con el flag de dato válido "V", al igual que ya hicimos en el caso de la Calidad del Aire.  En este caso, y puesto que se trata sólo de una columna, pasaremos directamente a utilizar la función "ifelse", aunque bien podríamos utilizar como antes la función "mutate".

   #Reviso los datos en función del flag y muto los que no son válidos a NA.
   Sabadellmeteo$VALOR_LECTURA<-ifelse(Sabadellmeteo$CODI_ESTAT=="V",
   Sabadellmeteo$VALOR_LECTURA,NA)

Lo segundo que vamos a hacer es seleccionar una a una las variables meteorológicas que nos interesan. Podremos ver los metadatos de las estaciones meteorológicas, como os decía, en la propia web de datos abierto, para nuestro ejercicio cogeremos no obstante los siguientes datos, a los que os adjunto el código al lado: velocidad del viento (30), dirección del viento (31), temperatura (32), humedad relativa (33), presión(34), precipitación (35), y la radiación solar total (36).

Recordaréis que para seleccionar y dividir un conjunto de datos utilizábamos la función "subset", que nos permitía exponer un criterio de selección, que en nuestro caso será el parámetro que queremos extraer, como por ejemplo "CODI_VARIABLE==30" para extraer la velocidad, y que además nos permitía seleccionar los campos que queríamos extraer, que en nuestro caso bastará con los campos de fecha y valor de lectura.

Aprovecharemos que extraemos el archivo para cambiar los nombres de los campos, de forma que luego podamos agruparlos bien en un sólo archivo, con el comando "names". Y ya que estamos, aprovecharemos para convertir la fecha, que nos hemos descargado como un factor, en una campo de fecha indexado que luego nos pueda venir bien para trabajar. Utilizaremos para ello el paquete "lubridate" que como sabéis es un paquete muy útil de R, y dada la fecha del archivo, seleccionaremos la función "dmy_hms" que ya reconoce automáticamente si es AM / PM y y nos da la fecha final en formato adecuado para el manejo de R.

El código quedaría conforme muestro a continuación, teniendo en cuenta que habría que repetirlo para  cada parámetro:

   #seleccionamos los datos uno a uno en función del parámetro que nos interesa.
   #Llamo a la librería para el manejo de fechas
   library(lubridate)
   #saco la velocidad del viento con la fecha del dataframe general
   ws<-subset(Sabadellmeteo, CODI_VARIABLE==30,select=c("DATA_LECTURA",
   "VALOR_LECTURA"))
   #le cambio los nombres a las columnas
   names(ws)<-c("date","ws")
   #uso lubridate para conseguir que la columna tenga formato de fecha
   ws$date<-dmy_hms(ws$date)

   #Hago lo mismo para la dirección del viento.
   wd<-subset(Sabadellmeteo, CODI_VARIABLE==31, select=c("DATA_LECTURA",
   "VALOR_LECTURA"))
   names(wd)<-c("date","wd")
   wd$date<-dmy_hms(wd$date)

   (...)

Tenemos que juntar ahora los archivos por la fecha. Ya sabemos que para ello tenemos la función "merge" que ya utilizamos en la Píldora 1 para juntar los datos de calidad del aire en un sólo archivo de dos en dos. En esta ocasión tenemos un montón de data frame, por lo que puede ser un dolor de cabeza andar haciendo combinaciones a pares. Esta es la oportunidad perfecta para hablar del comando base "Reduce".

El comando "Reduce" te permite usar una función binaria como "merge" para combinar elementos de una lista o vector, que en nuestro caso serán nuestras data frame. Para ello podemos utilizar la siguiente combinación:

   #Genero un data frame de meteorológicas de Sabadell extraídas por fecha.
   meteo<-Reduce(merge, list(wd,ws,t,p,r,rh,srt))
   #Borro los data frame intermedios que he creado. 
   rm(wd,ws,t,p,r,rh,srt,Sabadellmeteo)

Si nos fijamos bien en nuestro archivo de datos de meteorología ya unificado, podremos ver que los datos meteorológicos los tenemos en periodo "semihorario", cuando los datos de calidad del aire los tenemos agregados en horas. Si los unificásemos a las bravas, estaríamos comparando los datos semihorarios frente a los horarios, y estaríamos incluyendo un sesgo importante.

Es imperativo, por tanto, que ambos conjuntos de datos se encuentren en el mismo periodo de integración, y en este caso, lo suyo es convertir el archivo de meteorología en datos horarios, para lo cual utilizaremos la función "timeAverage" del paquete "Openair" que nos permite hacerlo de una forma sencilla y rápida.

Esta función, además, incorpora una opción específica que nos permite reconocer automáticamente la dirección y velocidad del viento y hacer el promedio vectorial. Para ello, eso sí, tendremos que tener "en inglés" nuestros campos, definiendo la velocidad del viento como "ws" (wind speed) y la dirección del viento como "wd" (wind direction).

   #Llamo a la librería Openair.
   library(openair)
   #Saco un dataframe calculando la media horaria de los parámetros meteorológicos
   meteoh<-timeAverage(meteo, avg.time="hour", vector.ws = TRUE)
   rm(meteo) #Borro el antiguo data frame de datos semihorarios.

Si lo hemos hecho bien, tendríamos que tener nuestro dataframe de "Sabadellaire" de la primera píldora, y nuestro dataframe de datos de meteorología de Sabadell una vez filtrado y calculado a datos horarios "meteoh". Con una última función "merge" mezclaremos ambos data frame, solo que en esta ocasión, dado que ambos conjuntos de datos pueden tener periodos temporales distintos, lo que nos interesa es decirle al comando que nos conserve los datos de partida de calidad del aire, con un "all.x=TRUE".

   #Unifico Calidad del Aire y meteorológicas.
   Sabadellfinal<-merge(Sabadellaire,meteoh, by="date", all.x=TRUE)
   rm(meteoh) #borro el de meteorológica de transición.

Tendremos finalmente un data frame de datos de calidad del aire y de meteorología para Sabadell que podremos empezar a trabajar tranquilamente en R.

Con los datos de contaminación y los datos de meteorología en el mismo archivo, ya podremos empezar a estudiar la calidad del aire y a trabajar los datos para conocer cómo evoluciona la estación y qué peculiaridades presenta. Podremos, por ejemplo, comparar parámetros meteorológicos, como la radiación solar, con parámetros vinculados, como el ozono, consiguiendo ver gráficas de evolución temporal como la que se expone.

Pero esto ya queda para próximas píldoras de R....

 

Esta segunda píldora va dedicada a Ricardo Ruiz. Felicidades por ese sobresaliente. Que nadie te quite nunca esa ilusión por aprender cosas nuevas y hacer siempre el mejor trabajo.