Skip navigation

Medi ambient

Objectius

1. Crear un pòster en anglès analitzant les dades de contaminació atmosfèrica i metereològiques des del 1991 fins al 2021 hora a hora, dia a dia.

2. Treure conclusions científiques amb estadística i gràfics avançats amb Openair R i RStudio

1. Obtenció de les dades

Les dades es poden trobar googlejant "dades obertes atmòsfera" i "dades obertes meteocat":

Dades de contaminació atmosfèrica

Dades meteorològiques

Hem de tenir en compte que les dades són horàries per contaminació i les dades meteorològiques són semihoràries (cada 30 minuts) que hem de combinar en un únic dataframe per fer gràfics que relaciones ambdós tipus de dades. Les dades de meteo seràn de l'estació més propera a la nostra de contaminació.

De les dades meteorològiques els codis de les variables que necessitem són 30 per wind speed que posarem ws, i 31 per wind direction que posarem wd. Les dates seran en format ISO8601  2021-03-15 16:00:00 amb format POSIXct amb el nom date. Els noms date, ws i wd són requisits de Openair

2. Programari necessari

Es necessari instal.lar el programari RStudio i els paquets tidyverse i openair.

Tidyverse permet ordenar les dades. Recordeu la famosa dita científica que diu que el 80% estem ordenant les dades i el 20% de temps restant queixant-nos de la seva ordenació.

Openair permet fer estudis de l'aire, hora a hora, dia a dia, any a any de forma específica i avançada.

Per instal.lar aquest 2 paquets de R has d'escriure:

install.packages (c("tidyverse","openair"))

3. Ordenant les dades atmosfèriques

Hem de llegir les dades del nostre ordinador:

> city<-read.csv("C://Users/YOURCOMPUTERNAME/Documents/city.csv")

> View(city)

Canviem les hores de les columnes per les files

>city1<-pivot_longer(city,cols=c(h01,h02,h03,h04,h05,h06,h07,h08,h09,h10,h11,h12,h13,h14, h15,h16,h17,h18,h19,h20,h21,h22,h23,h24), names_to="hour", values_to = "value")

>city2<-city1[-c(1,2,4,6:16)]

>write.csv(city2,"C:\\Users\\YOURCOMPUTERNAME\\Documents\\city2.csv")

4. Modifica les dades atmosfèriques amb LibreOffice Calc

Esborrem T00.00.00.000 i substituim les hores h01, etc  per 01:00:00 fins obtenir les dates en format ISO:

martorell3

Combinem data i hora junts a city4.csv emprant LibreOffice Calc o RStudio

> city4 <- city3 %>% mutate(name=paste0(data, " ", hour))

Creem city5 juntant dia i hora sota el nom de columna date

city5

5. Convertim el format de les dades

library(openair)

city5NO2 <- subset(city5, pollutant=="NO2")

city5NO2$date<-as.POSIXct(city5NO2$date,"%Y-%m-%d %H:%M:%S", tz="Europe/Madrid")

class(city5NO2$date)

[1] "POSIXct" "POSIXt"

View(city5NO2)

S'ha de vigilar que la data no sigui un conjunt de caràcters (character) sino un POSIX, es a dir, una data
> class(city5NO2$date)
[1] "character"
> city5NO2$date<-as.POSIXct(city5NO2$date,"%Y-%m-%d %H:%M:%S", tz="Europe/Madrid")
> class(city5NO2$date)
[1] "POSIXct" "POSIXt"
> class(city5NO2$pollutant)
[1] "character"
city5h2s$pollutant<-as.factor(city5NO2$pollutant)
> class(city5NO2$pollutant)
[1] "factor"

6. Creant gràfics amb Openair R

> timeVariation(city5NO2, pollutant="value")
timevariationh2s
> trendLevel(city5NO2, pollutant = "value", main="Hydrogen sulfide evolution in MYCITYNAME")
trendLevel Martorell
> daily<-timeAverage(city5NO2,avg.time = "day")
> View(daily)
> calendarPlot(city%NO2, pollutant="value", year="2020")
NO2 Martorell 2020
> yearly<-timeAverage(city5no2,avg.time = "year")
> View(yearly)
date                 annual mean NO2
1991-01-01      41.99301
1992-01-01      31.39408
1993-01-01      42.22537
1994-01-01      49.77212
1995-01-01      74.34786 
1996-01-01      60.54973
1997-01-01       49.65882
1998-01-01      47.21167
1999-01-01      40.91106
2000-01-01     40.07582
2001-01-01      30.32435
2002-01-01      32.65402
2003-01-01      34.38977
2004-01-01      36.17926
2005-01-01      39.91667
2006-01-01      NA
2007-01-01       40.96774
2008-01-01      17.58333
2009-01-01      40.40244
2010-01-01       59.78571
2011-01-01       43.76190
2012-01-01       39.73099
2013-01-01       40.99248
2014-01-01       36.41213
2015-01-01       41.03577
2016-01-01       34.98688
2017-01-01        36.85954
2018-01-01       31.52100
2019-01-01       29.44134
2020-01-01       21.85310
2021-01-01        32.40141 (fins 11 març 2021)

Air quality standards indica diferents nivells límits de contaminació, p.e. per NO2; els nivells han de ser inferior a 40 µg/m3
Veure la informació en el pdf al final d'aquesta pàgina
Comprovar per exemple el número de superacions de les mitjanes octohoràries de l'ozó amb aquest codi:
library(openair)
episode<-selectRunning(city6, pollutant="O3",threshold=120, run.len=8)
nrow(episode)
> nrow(episode)
[1] 95
R ens diu que s'ha superat el nivell de 120 micrograms/metre cúbic en la mitjana de 8h un total de 95 vegades.

7. Ordenant les dades amb pivot_wider

> library(tidyverse)

> library (openair)

> city6<-pivot_wider(city5, names_from= pollutant, values_from =value)

> View (city6)

> write.csv(city6,"C:\\Users\\ELMEUNOM\\Documents\\LAMEVACIUTAT\\city6.csv")

pivotwider

8. Variació temporal contaminants 1991-2021

class(city6$date)
[1] "character"
city6$date<- as.POSIXct(city6$date,format="%Y-%m-%d %H:%M:%S",tz="Europe/Madrid")

timeVariation(city6, pollutant=c("O3","NO2","H2S","NO","HCNM","CO","SO2","HCT", "NOX","PM10"), main="Air pollution in MYCITYNAME (1991-2021)")

timeVariation

9. Ordenant les dades meteorològiques

> wind<-read.csv("C://Users/YOURCOMPUTERNAME/Documents/wind.csv")

> View(wind)

wind1←wind[-c(1,2,5,7,8)]

wind2<-pivot_wider(wind1,names_from = CODI_VARIABLE, values_from = VALOR_LECTURA)

names(wind2)[names(wind2) == "31"] <- "wd"

names(wind2)[names(wind2) == "30"] <- "ws"

names(wind2)[names(wind2) == "DATA_LECTURA"] <- "date"

write.csv(wind2,"C:\\Users\\YOURCOMPUTERNAME\\Documents\\wind2.csv")

10. Unir les dos bases de dades (contaminants i meteo)

Les dades de descarreguem de la base de dades de meteorologia indicades en l'apartat anterior son semihoràries i expressades en AM i PM i necessitem procedir amb Cerca i Reemplaça de Calc perquè tingui el mateix format que les dades de contaminació.

Amb library (openair) i wind6<-timeAverage(wind5, time.avg="hour") creem les dades horàries a partir de semihoràries.

Hem d'estar segurs que la classe de la date de la base de dades del vent es una data tipus POSIXct per poder combinar-la

Es pot unir amb diferents instruccions per exemple amb una d'openair:

> cityall<-merge(city6, wind6, by ="date")

> View (cityall)

I ja podem fer pollutionRose i veure d'on vé la contaminació.

11. Analitzar les dades futures

Podem fer-ho amb diferents biblioteques de R per fer "forecast" o previsió de dades futures a partir de dades presents:

tot2<-read.csv("C://Users/ELTEUNOM/Documents/LATEVACIUTAT/tot2.csv")
View(tot2)
m<-prophet(tot2)
future <- make_future_dataframe(m, periods = 365)
forecast <- predict(m, future)
plot(m,forecast)

library(prophet)
prophet_plot_components(
+     m,
+     forecast,
+     uncertainty = TRUE,
+     plot_cap = TRUE,
+     weekly_start = 0,
+     yearly_start = 0,
+     render_plot = TRUE
+ )

View (forecast) es veuen totes les dades previstes anteriors i futures amb tots els paràmetres i components temporals emprats per l'algorisme prophet. PODEM COMPROVAR SI ELS VALORS PREVISTOS S'HAN COMPLERT I AMB QUIN ERROR RESPECTE ELS NOUS VALORS REALS MESURATS

Per veure les dades finals:

tail(forecast[c('ds', 'yhat',
+ 'yhat_lower', 'yhat_upper')])

TAULA OBTINGUDA: amb valors previstos (yhat) , previst màxim (yhat_lower) i previst mínim (yhat_lower) 

ds yhat yhat_lower yhat_upper
77223 2022-03-04 23:00:00 36.04179 13.277273 58.87805
77224 2022-03-05 23:00:00 25.11000 0.532229 49.23281
77225 2022-03-06 23:00:00 26.99720 3.531979 50.60021
77226 2022-03-07 23:00:00 34.88387 11.911203 58.48056
77227 2022-03-08 23:00:00 33.60286 8.947426 56.70240
77228 2022-03-09 23:00:00 35.81467 12.322405 58.83855