Ocena prognoz jakości powietrza
Zarządzanie emisjami i jakością powietrza
Adriana Szulecka, szulecka@agh.edu.pl
05 listopad, 2019
Wprowadzenie
Na podstawie powszechnie dostępnych w Internecie krótkoterminowych prognoz jakości powietrza należy wykonać krótką ocenę statystyczną uzyskanych danych. Obliczenia polegają na porównaniu wartości prognostycznych z rzeczywistymi pomiarami realizowanymi na stacjach monitoringu ciągłego jakości powietrza w regionie przy użyciu wybranych metod statystycznych.
Prognozy jakości powietrza dla Małopolski
Małopolska w zdrowej atmosferze
Dane prognoz krótkoterminowych dla obszaru Małopolski i poszczególnych miast w jego zasięgu zamieszczane są na stronie Małopolska w zdrowej atmosferze od roku 2016. Wcześniej prognozy te były publikowane w ramach serwisu Wrota Małopolski od roku 2011. Obecnie do ich przygotowania wykorzystywany jest model GEM-AQ przez Instytut Ochrony Środowiska - Państwowy Instytut Badawczy. Dla obszaru województwa wykonywane są mapy maksymalnej wartości indeksu jakości powietrza w ciągu dnia. Dla poszczególnych miast w Małopolsce na stronie udostępniane są przewidywane 24-godzinne stężenia substancji (PM10, PM2,5 i O3) osobno na trzy nadchodzące doby (dzisiaj, jutro i pojutrze). Wyniki prognozowania aktualizowane są codziennie.
EKO-prognoza
W serwisie EKO-prognoza możemy znaleźć mapy indeksu jakości powietrza dla całego kraju, a także Europy, również sporządzane w oparciu o model fotochemii GEM-AQ. Wytwarzaniem i nadzorowaniem prognoz krótkoterminowych zajmuje się fundacja EKO-prognoza od roku 2009.
System FAPPS
System prognozowania rozprzestrzeniania zanieczyszczeń powietrza - FAPPS to system prognoz krótkoterminowych dla obszaru Małopolski oraz Krakowa dla czterech substancji zanieczyszczających: PM10, PM2,5, NO2 oraz SO2. W serwisie przedstawiane są średnie 24-godzinne oraz 1-godzinne stężenia zanieczyszczeń w postaci map. Wytwarzaniem oraz utrzymywaniem prognoz zajmuje się IMGW-PIB od roku 2013, a cały system powstał w oparciu o własny projekt rozwojowy finansowany z Narodowego Centrum Badań i Rozwoju. Prognozowanie odbywa się w systemie kilku modeli meteorologicznych (ALADIN oraz MM5), preprocesora meteorologicznego CALMET oraz właściwego dyspersyjnego modelu obłoku CALPUFF.
Ocena statystyczna prognoz jakości powietrza
Projekt wykonamy w oparciu o typowe miary błędów prognoz ex post (wygasłych), które mierzą w jakim stopniu wartości prognozowane odchylone są od rzeczywistych danych pomiarowych. Są to:
- FAC2 - ilość wartości prognozowanych spełniających warunek:
\[\begin{equation*} 0,5 \le \frac{P_{i}}{O_{i}} \le 2,0 \end{equation*}\] - MB - Mean Bias - Błąd Średni,
\[\begin{equation*} MB = \frac{1}{n} \sum_{i=1}^n P_{i}-O_{i} \end{equation*}\] - MGE - Mean Gross Error - Błąd Średni Bezwzględny,
\[\begin{equation*} MGE = \frac{1}{n} \sum_{i=1}^n |P_{i}-O_{i}| \end{equation*}\] - NMB - Normalised Mean Bias - Znormalizowany Błąd Średni,
\[\begin{equation*} NMB = \frac{\sum_{i=1}^n P_{i}-O_{i}}{\sum_{i=1}^n O_{i}} \end{equation*}\] - NMGE - Normalised Mean Gross Error - Znormalizowany Błąd Średni Bezwzględny,
- RMSE - Root Mean Squared Error - Pierwiastek Błędu Średniokwadratowego, \[\begin{equation*} RMSE = \left(\frac{\sum_{i=1}^n ({P_{i}-O_{i}})^2}{n}\right)^\frac{1}{2} \end{equation*}\]
gdzie \(P_{i}\) oznacza kolejną wartość prognozowaną, \(O_{i}\) oznacza kolejną wartość obserwowaną, a \(n\) jest liczebnością par prognoza-obserwacja. Ponadto policzymy:
współczynnik korelacji Pearsona (r),
współczynnik korelacji Spearmana (rho) - dla porównania.
Oceny siły korelacji będziemy dokonywać w oparciu o poniższą skalę:
- r < 0,1 – brak korelacji,
- 0,1 < r < 0,3 – korelacja słaba,
- 0,3 < r < 0,5 – korelacja umiarkowana,
- 0,5 < r < 0,7 – korelacja wysoka,
- 0,7 < r < 0,9 – korelacja bardzo wysoka,
- 0,9 < r < 1 – korelacja pełna.
W celu umożliwienia wizualnej interpretacji wyników i sporządzania dodatkowych wniosków wspomożemy się:
wykresami trendu serii czasowych porównujących wartości w czasie,
wykresami rozrzutu, gdzie na osi X znajdować będą się obserwacje a na osi Y - wartości prognozowane,
diagramem Taylora, który pozwala na wybranie najlepszej prognozy (dzisiaj, jutro lub pojutrze),
wykresem kwantylowym, który pokazuje dla jakich wartości stężeń obserwowane są największe odchylenia od modelu idealnego.
Pobieranie i przygotowanie danych dotyczących prognoz i pomiarów
Pobieranie prognoz krótkoterminowych dla wybranych miast w Małopolsce
Ze względu na to, że jedynie w serwisie Małoposka w zdrowej atmosferze udostępniane są dane liczbowe dotyczące prognozowanych stężeń zanieczyszczeń, to dzisiaj skupimy się na tych właśnie danych (analiza danych w postaci w mapy jest bardzo problematyczna).
Dostęp do prognoz jakości powietrza na stronie Małopolska w zdrowej atmosferze można uzyskać po kliknięciu w odpowiednie miasto na liście, a następnie przejściu do pozycji “Prognozy średniodobowej jakości powietrza” na trzy kolejne dni. Oprócz ogólnej jakości powietrza, w rozwijanej liście “Pokaż szczegóły” możemy dostać informacje dotyczące prognozowanych średniodobowych stężeń substancji.
Ponieważ prognozy są codziennie aktualizowane, to ważna jest systematyczność przy ich kolekcjonowaniu. Aby ułatwić pracę skorzystamy z pliku przygotowanego przez Prowadzącego, o nazwie prognozy_2018.xlsx dla kilku wybranych tygodni. W pliku w poszczególnych arkuszach kalkulacyjnych znajdziemy prognozy stężenia średniodobowego dla trzech substancji: PM10, PM2,5 oraz O3 dla dziewięciu miast w Małopolsce, dla których równolegle prowadzone są pomiary zanieczyszczeń powietrza (lub w ich pobliżu):
- Kraków
- Tarnów
- Nowy Sącz
- Chrzanów
- Gorlice
- Nowy Targ
- Sucha Beskidzka
- Wieliczka
- Zakopane
Pobieranie pomiarów z odpowiednich stacji pomiarowych WIOŚ
Ponieważ obserwacji nie jest dużo, to przy pobieraniu danych pomiarowych z odpowiednich stacji skorzystamy ze strony WIOŚ w Krakowie. Wykorzystamy pomiary automatyczne, które są raportowane na bieżąco na stronie internetowej. Wybieramy dane miesięczne dla odpowiedniej stacji, a następnie zapisujemy plik z pomiarami do pliku .csv.
Można wykonać wstępną edycję pliku .csv w celu ułatwienia dalszej pracy (zmiana nazw kolumn, poprawienie daty).
Przygotowanie danych
Przygotowanie danych opiera się w głównej mierze na wczytaniu ich i uporządkowaniu w programie RStudio. Ze względu na specyficzny wygląd pliku z prognozami jest to konieczne - musimy połączyć je z pomiarami na podstawie odpowiedniej daty. Robienie tego ręcznie np. w programie Excel byłoby bardzo żmudne i niepotrzebne.
Dalsza praca będzie odbywać się w oparciu o jeden z zestawów prognoz - dla miasta Nowy sącz i susbtancji PM10.
# Ustawianie folderu roboczego
setwd("E:/01_Ocena_prognoz/2018_2019")
#Instalowanie i wczytywanie bibliotek - również z myślą o analizach statystycznych
library(readxl)
#install.packages("openair")
library(openair)
library(tidyr)
library(ggplot2)
#install.packages("ggpmisc")
library(ggpmisc)
#Wczytywanie arkusza z Excela i zmiana formatu daty
nowy_sacz_p <- read_excel("prognozy_2018.xlsx", sheet = "Nowy Sącz")
nowy_sacz_p$Data <- as.Date(nowy_sacz_p$Data)
#Wczytywanie z pliku *.csv i zmiana formatu daty
nowy_sacz_o <- read.csv("nowy_sacz_o.csv", sep=";", as.is=T)
nowy_sacz_o$Data <- as.Date(nowy_sacz_o$Data, "%d.%m.%Y")
#Sprawdzenie typów danych
str(nowy_sacz_o)## 'data.frame': 28 obs. of 2 variables:
## $ Data : Date, format: "2018-10-24" "2018-10-25" ...
## $ PM10_o: int 13 12 50 50 31 19 34 39 69 65 ...str(nowy_sacz_p)## Classes 'tbl_df', 'tbl' and 'data.frame': 90 obs. of 5 variables:
## $ Data : Date, format: "2018-10-24" "2018-10-25" ...
## $ Prognoza: chr "Dzisiaj" "Jutro" "Pojutrze" "Dzisiaj" ...
## $ PM10_p : num 5 5 7 6 6 26 8 23 11 8 ...
## $ PM2.5_p : num 4 6 7 6 7 29 8 25 12 8 ...
## $ O3_p : num 46 47 49 43 54 37 52 45 26 52 ...#Połączenie obu zestawów danych po dacie
nowy_sacz <- merge(x=nowy_sacz_p, y=nowy_sacz_o, by="Data", all.x=T)
write.csv(nowy_sacz, "nowy_sacz_op.csv")W tym miejscu można kontynuować pracę w R (zalecany i szybszy sposób) bądź przejść do pracy w Excelu (wymagane ręczne wpisywanie wzorów).
Analiza statystyczna danych w R
Policzenie miar błędów i współczynników korelacji
Policzenie wszystkich wspomnianych wcześniej miar błędów można zrealizować w R za pomocą jednej funkcji z pakietu openair:
modStats(nowy_sacz, obs="PM10_o", mod="PM10_p", type="Prognoza", statistic = c("n", "FAC2", "MB", "MGE", "NMB", "NMGE", "RMSE", "r"))Dzięki rozdzieleniu słowem kluczowym type typu prognozy, dostajemy od razu odpowiedź jak wyglądają poszczególne współczynniki w zależności od rodzaju prognoz. Funkcja modStats zwraca również wartości współczynnika korelacji Pearsona, ale jednak bez określenia jego istotności statystycznej.
Do policzenia współczynników korelacji polecam funkcję cor.test z pakietu stats:
cor.test(x=nowy_sacz$PM10_o[nowy_sacz$Prognoza == "Dzisiaj"], y=nowy_sacz$PM10_p[nowy_sacz$Prognoza == "Dzisiaj"],
method = 'pearson')##
## Pearson's product-moment correlation
##
## data: nowy_sacz$PM10_o[nowy_sacz$Prognoza == "Dzisiaj"] and nowy_sacz$PM10_p[nowy_sacz$Prognoza == "Dzisiaj"]
## t = 4.4103, df = 26, p-value = 0.0001594
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3718590 0.8257299
## sample estimates:
## cor
## 0.6541802cor.test(x=nowy_sacz$PM10_o[nowy_sacz$Prognoza == "Dzisiaj"], y=nowy_sacz$PM10_p[nowy_sacz$Prognoza == "Dzisiaj"],
method = 'spearman')##
## Spearman's rank correlation rho
##
## data: nowy_sacz$PM10_o[nowy_sacz$Prognoza == "Dzisiaj"] and nowy_sacz$PM10_p[nowy_sacz$Prognoza == "Dzisiaj"]
## S = 1318.5, p-value = 0.000251
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho
## 0.6391557Za pomocą tej funkcji możemy policzyć zarówno współczynnik korelacji Pearsona jak i Spearmana. W odpowiedzi dostaniemy również wartość tzw. p-value, które informuje nas o istotności policzonego parametru. Nie wchodząc w szczegóły, jeżeli p-value jest poniżej 5% (0,05) to uznajemy, że współczynnik korelacji jest istotny statystycznie (przy przyjętym poziomie istotności 5%).
Niestety obliczenia współczynników musimy wykonać dla każdego typu prognozy osobno i spisać ręcznie wyniki.
Sporządzenie wykresów
Do sporządzenia wykresów trendu oraz rozrzutu wykorzystamy pakiet ggplot2.
Wykres rozrzutu
Ponieważ mamy różne substancje oraz różne typy prognoz to należy zastanowić się nad najlepszą metodą wizualizacji danych tak, aby widoczne były różnice ułatwiające sporządzenie wniosków dotyczących sprawdzalność prognozowanych wartości stężeń.
Możemy zamieścić wykresy obok siebie:
ggplot(data = nowy_sacz, aes(x = PM10_o, y = PM10_p)) +
geom_point(shape=1, size=5) +
geom_smooth(method=lm, se=FALSE, formula=y~x-1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0.25,
formula = y~x-1, parse = TRUE, size = 6) +
facet_grid(rows=vars(Prognoza)) +
scale_y_continuous(name = 'Prognoza PM10 [μg/m3]',
limits = c(0,40),
expand = c(0, 0)) +
scale_x_continuous(name = 'Pomiar PM10 [μg/m3]',
limits = c(0, 100),
expand = c(0, 0)) +
theme_minimal()
Możemy je też nałożyć na siebie i oznaczyć innym kolorem:
ggplot(data = nowy_sacz, aes(x = PM10_o, y = PM10_p, color=Prognoza)) +
geom_point(shape=1, size=5) +
geom_smooth(method=lm, se=FALSE, formula=y~x-1) +
stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
label.x.npc = "right", label.y.npc = 0.25,
formula = y~x-1, parse = TRUE, size = 8) +
scale_y_continuous(name = 'Prognoza PM10 [μg/m3]',
limits = c(0,40),
expand = c(0, 0)) +
scale_x_continuous(name = 'Pomiar PM10 [μg/m3]',
limits = c(0, 100),
expand = c(0, 0)) +
theme_minimal()
Wykres trendu serii czasowej
Dla serii czasowej zrobimy trzy wykresy dla prognoz na dzisiaj, jutro i pojutrze, z zaznaczeniem poziomu dopuszczalnego PM10 oraz linii pokazującej wartości stosunku prognozy do pomiaru (czyli bezpośrednio warunku zawartego we wzorze na parametr FAC2).
ggplot(data = subset(nowy_sacz, Prognoza == "Dzisiaj"), aes(x = Data)) +
geom_line(aes(y=PM10_o, color="Pomiar"), size=1.1) +
geom_line(aes(y=PM10_p, color="Prognoza na 'Dzisiaj'"), size=1.1) +
geom_line(aes(y=(PM10_p/PM10_o)*20, color="Stosunek prognoza/pomiar"), size=0.5) +
geom_hline(yintercept = 50, color="red", size=1) +
scale_x_date(date_labels = "%m.%d", breaks="1 day") +
scale_y_continuous(name = 'Stężenie PM10 [μg/m3]',
limits = c(0, 100),
expand = c(0, 0),
breaks = seq(0, 100, by = 20),
sec.axis=sec_axis(~./20,
name='Prognoza/pomiar',
breaks = seq(0, 5, by = 0.5))) +
theme_minimal() +
theme(legend.position = "bottom", legend.direction = "horizontal") +
scale_colour_manual("",
breaks = c("Prognoza na 'Dzisiaj'", "Pomiar","Stosunek prognoza/pomiar"),
values = c("darkgreen", "orange", "blue")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Tutaj zielona linia oznacza pomiary, pomarańczowa - prognozy na “dzisiaj”, przerywana niebieska - stosunek prognozy do pomiaru, a czerwona - 24-godzinny poziom dopuszczalny PM10.
Diagram Taylora
Diagram Taylora pozwala na pokazanie trzech uzupełniających się statystyk jednocześnie: współczynnika korelacji Pearsona R, odchylenia standardowego σ oraz scentrowanego błędu średniokwadratowego RMSE. Takie graficzne przedstawienie kilku modeli naraz pozwala na ich względne porównanie i stwierdzenie, których z nich jest najlepszy. Aby wykreślić w łatwy sposób diagram Taylora w R, skorzystamy z funkcji pakietu openair.
TaylorDiagram(nowy_sacz, obs="PM10_o", mod="PM10_p", group="Prognoza")
Jak widać na wykresie powyżej, najlepszy model powinien znajdować się w miejscu fioletowej kropki z napisem “observed” - tak by było, gdyby wartości prognozowane były równe obserwacjom. Na tej podstawie możemy stwierdzić, że prognozy na poszczególne dni nie różnią się znacząco między sobą i charakteryzują się podobnymi wartościami odchylenia standardowego. Jednakże, nieco lepiej pod względem wartości współczynnika korelacji R oraz błędu RMSE sprawdza się prognoza na “dzisiaj”, następnie na “jutro”, a najgorzej na “pojutrze”, co potwierdza intuicyjne wnioski.
Warunkowy wykres kwantylowy
Wykres, którym się posłużymy nie jest stricte wykresem kwantylowym, ponieważ znajduje się na nim o wiele więcej informacji. Pokazuje on zarówno rozkłady obserwowanych i prognozowanych stężeń w postaci histogramów, jak i linię modelu idealnego i tego, który analizujemy. Sporządźmy warunkowy wykres kwantylowy dla wszystkich trzech rodzajów prognoz jakimi dysponujemy.
conditionalQuantile(nowy_sacz, obs="PM10_o", mod="PM10_p", type="Prognoza", min.bin=2, bins=10, xlab="Prognozowana wartość", ylab="Obserwowana wartość")
Jak widać, linie wyznaczone przez prognozowane wartości wraz z obszarami kolejnych percentyli nie pokrywają się z linią modelu idealnego. Znajdują się ponad nią, co oznacza, że model jest mocno niedoszacowany. Jednocześnie wartości, jakie przyjmuje prognoza, są o wiele niższe niż obserwacje, co potwierdza wnioski zauważone wcześniej. Różnice między poszczególnymi prognozami nie są znaczące, jednakże można powiedzieć, że linia prognozy na “dzisiaj” wydaje się być najbliżej modelu idealnego w porównaniu do dwóch pozostałych.
Zakres projektu
Kompletny projekt powinien zawierać:
- Plik .doc/.docx/.pdf zawierający:
- krótki opis systemu prognoz realizowanych na stronie Małopolska w zdrowej atmosferze;
- krótki opis stacji pomiarowej (lokalizacja, mierzone substancje), na której rejestrowano pomiary wykorzystane w obliczeniach (punkty a i b razem na maksymalnie 1 stronę);
- wyszczególnienie, jakich substancji dotyczą obliczenia;
- wyniki wartości współczynników korelacji Pearsona (r) oraz Spearmana (rho) wraz z podaniem, czy wyniki są istotne statystycznie;
- tabelę zbiorczą z wynikami miar błędów dla wszystkich ocenianych substancji;
- dla każdej substancji wykresy zmienności stężeń prognozowanych i zmierzonych w analizowanym okresie (na osi x – data, na osi y – stężenie [jednostka]); dla pyłu PM10 i SO2 zaznaczyć dobowy poziom dopuszczalny na wykresie;
- do każdego wykresu zmienności stężeń dołączyć wykres rozrzutu wartości wraz z linią trendu (na osi x – stężenie zmierzone na stacji, na osi y – prognozowane stężenie, początek linii trendu w punkcie 0,0);
- wykreślone diagramy Taylora i wykresy kwantylowe z krótkimi komentarzami dotyczącymi wyboru najlepszej prognozy wraz z uzasadanieniem;
- krótkie podsumowanie uzyskanych wyników z wyszczególnieniem (jeśli wystąpiły) dni, w których nastąpiło ponad dwukrotne niedoszacowanie/przeszacowanie wartości stężeń danej substancji, określeniem siły korelacji danych oraz ogólną oceną działania modelu (również wizualną na podstawie wykresów).
- Plik .R/.xls/.xlsx bądź inny, który zawiera:
- serię danych prognozowanych i pomiarowych wraz z przypisaną datą dla każdej z analizowanych substancji;
- obliczenia cząstkowe i wyniki końcowe;
- wykresy zmienności stężeń w czasie oraz wykresy rozrzutu;
- funkcje diagramu Taylora i wykresu kwantylowego (w przypadku pliku .R).
W przypadku skryptu .R wymagane są również działające obliczenia wraz z komentarzami, a także dostarczenie plików wczytywanych do środowiska R.
Gotowy projekt należy przesyłać w wersji elektronicznej na adres szulecka@agh.edu.pl w terminie do dwóch tygodni od daty zajęć (tj. 20.11.2019, 23:59).
Literatura uzupełniająca
[1] Chang J. C., Hanna S. R.; Technical Descriptions and User’s Guide for the BOOT Statistical Model Evaluation Software Package, Version 2.0, 2005.
[2] Derwent D., Fraser A., Abbott J., Jenkin M., Willis P., Murrells T.; Evaluating the Performance of Air Quality Models, 2010.
[3] Kamiński J. W., Neary L., Strużewska J., McConnell J. C., Lupu A., Jarosz J., Toyota K., Gong S. L., Côté J., Liu X., Chance K., Richter A.; GEM-AQ, an on-line global multiscale chemical weather modelling system: model description and evaluation of gas phase chemistry processes [W:] Atmos. Chem. Phys., 8, s. 3255-3281, 2008.
[4] Strużewska-Krajewska J., Kamiński J. W., Durka P., Jefimow M.; Sprawdzalność operacyjnej prognozy stanu jakości powietrza dla obszaru Województwa Małopolskiego w okresie październik 2012 – wrzesień 2013, Warszawa 2013.
[5] Szczygłowski, P.; Ocena przydatności wybranych modeli gaussowskich w obliczeniach stanu zanieczyszczenia powietrza, Dysertacja doktorska, AGH, WGGiIŚ, Kraków 2007.
[6] Szulecka, A., Mazur, M.; Application of the Statistical Error and Quantitative Performance Measures in the Evaluation Process of Short-Term Air Quality Forecasts for Krakow (Poland) [W:] Geomatics and Environmental Engineering, Vol. 10, no. 3, s. 87-99, 2016.
[7] Carslaw D. C., K. Ropkins; openair — an R package for air quality data analysis [W:] Environmental odelling & Software, vol. 27-28, s.52–61, 2012.
[8] Carslaw D.; The openair manual — open-source tools for analysing air pollution data. Manual for version 1.1-4, King’s College London, 2015.