Model Performance Statistics for meteorogical data (MPSmet)

Title of the article: …

Posted in …

aAGH University of Science and Technology, Faculty of Mining Surveying and Environmental Engineering, Department of Environmental Management and Protection, Av. A. Mickiewicz 30, Krakow, Poland, contact: rzeszut@agh.edu.pl.

2. Przygotowanie danych

Wczytujemy zestawy danych z pliku binarnego .RData.

Podgląd danych. Ponizsze dane zawierają informację o wynikach pomiarów prowadzonych w ramach eksperymentów MC, LV i TR oceny modeli dyspersji. W ramach tych eksperymentów prowadzone były pomiary wybranych parametrów meteorologicznych (Prędkości i kierunku wiatru oraz temperatury) w pobliżu źródła emisji. W ramach nieniejszego opracowania wykorzystamy te dane do oceny modelu WRF. Symulacje modelem WRF były prowadzone w dwóch rozdzielczosciach 1 km WRF1 i 4 km WRF4. Ponadto pozyskano również dane z stacji synoptycznych, które znajdowały się wewnątrz domeny obliczeniowej. Done te pozyskano z bazy danych [ISD] (https://www.ncdc.noaa.gov/isd).

W obiektach typu tibble tj. synop i synop_g znjadują się zestawione dane pomiarowe i wyniki symulacji. Dane zostały przygotowane od razu w postaci dwóch reprezentacj tj. wąskiej i szerokiej. Wiecej na temat reprezentacji danych i sposoób przetwarzania baz danych w R za pomocą pakiety tidy możesz poczytać w książce R for Data Science.

Zmienne:

  • date - data wykonania pomiaru
  • MED - oznaczenie eksperymentu
  • MOD - oznaczenie rozdzielczości siatki obliczeniowej dla prognoz z modelu WRF
  • ws, wd, tt - odpowiednio prędkość, kierunke wiatru i temperatura. Indeksy o - obserwacje, indeks m prognozy.
  • station - Skrócona nazwa stacji, jeśli dane pozyskane z oksperymentu (onsite-data) to oznaczenie ONS.
  • station_typ - przyjmyje wartosci ISD dla stacji synoptycznych, a ONS dla danych pozyskanych w ramach eksperymentu.
  • id_var - połącznie zmiennych grupujących MED, station, MOD, w wartosci unikalne.

Podgląd danych:

# A tibble: 6 x 12
  date                MED   MOD    ws_m  wd_m   tt_m  ws_o  wd_o  tt_o
  <dttm>              <chr> <chr> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>
1 1988-01-01 00:00:00 LV    WRF1   4.49  242. -0.650   2.6   270   3.3
2 1988-01-01 01:00:00 LV    WRF1   5.06  242. -0.25    5.1   240   3.3
3 1988-01-01 02:00:00 LV    WRF1   5.21  249   0.95    2.6   250   3.3
4 1988-01-01 03:00:00 LV    WRF1   4.8   254   1.35    3.6   240   3.3
5 1988-01-01 04:00:00 LV    WRF1   3.94  253.  1.25   NA      NA  NA  
6 1988-01-01 05:00:00 LV    WRF1   3.51  248.  1.25   NA      NA  NA  
# ... with 3 more variables: station <chr>, station_typ <chr>,
#   id_var <chr>

Wąska reprezentacja zawiera analogiczne inf. Dane różnią się tylko strukturą. Dane są przechowywane w kolumnach obs i mod. Czynniki meteorologiczne można zidentyfikować po zmiennej param, którea przechowuje ich oznaczenia tj. ws, wd i tt.

Podgląd danych

# A tibble: 6 x 9
  date                MED   MOD   station station_typ id_var param   mod
  <dttm>              <chr> <chr> <chr>   <chr>       <chr>  <chr> <dbl>
1 1984-08-06 00:00:00 TR    WRF1  RTIA    ISD         TR-RT~ tt    17.2 
2 1984-08-06 00:00:00 TR    WRF1  RTIA    ISD         TR-RT~ wd    75.8 
3 1984-08-06 00:00:00 TR    WRF1  RTIA    ISD         TR-RT~ ws     2.79
4 1984-08-06 00:00:00 TR    WRF4  RTIA    ISD         TR-RT~ tt    16.0 
5 1984-08-06 00:00:00 TR    WRF4  RTIA    ISD         TR-RT~ wd    68.2 
6 1984-08-06 00:00:00 TR    WRF4  RTIA    ISD         TR-RT~ ws     3.34
# ... with 1 more variable: obs <dbl>

3. Funkcje oceny modeli (stats param)

Funkcję wczytujemy z pliku źródłowego zgodnie z poniższym przykładem. Upewnij się, że ten plik znajduje się w folderze roboczym. np. poleceniami getwd i dir.

Powyższe polecenie spowodowało wykonanie kodu zawartego w pliku. Teraz w oknie Environment oprócz przykłądowych danych znjadują się również

3.1. eval_all - Model Performance Statistics for met. data

Funkcja eval_all oblicza szereg statystycznych parametrów oceny modeli prognostycznych. Niniejsza funkcja była pisanna pod ocenę danych meteorologicznych, ale z można ją stosowac do oceny innych modeli. Funkja ta oblicza następujace parametry:

  • n = Count - NA
  • mean_o = Mean Observed
  • mean_m = Mean Predicted
  • SD_o = Standard Deviation Observed value
  • SD_m = Standard Deviation Predicted value
  • BIAS = Mean Bias (BIAS)
  • MAE = Gross Error (MAE - Mean Absolut Error)
  • RMSE = Root Mean Square Error
  • CRMSE = Centralised root-mean-square error
  • MFB = Mean fractional bias
  • MNB = Mean normalised bias
  • MNE = Mean normalised error
  • NMB = Normalised mean bias
  • NME = Normalised mean error
  • FAC2 = Fraction of predictions within a factor of two
  • IOA = The Index of Agreement based on Willmott et al. (2011)
  • R = The Pearson correlation coefficient

Opis powyższych parametrów znajdziesz w poniższych publikacjach, raportach, czy instrukcjach:

  • Cambridge Environmental Research Consultants Ltd., 2016. Model Evaluation Toolkit User Guide 1–63. link
  • Emery, C., Liu, Z., Russell, A.G., Odman, M.T., Yarwood, G., Kumar, N., 2017. Recommendations on statistics and benchmarks to assess photochemical model performance. J. Air Waste Manag. Assoc. 67, 582–598. link
  • Emery, C., Tai, E., Yarwood, G., 2001. Enhanced Meteorological Modeling and Performance Evaluation for Two Texas Ozone Episodes. Env. Int. Corp. 235.link
  • Lingard, J., Labrador, L., Brookes, D., Fraser, A., 2013. Statistical evaluation of the input meteorological data used for the UK air quality forecast (UK-AQF). link

Uwaga:: Parametry MFB, MNB oraz MNE wymagają ustalenia wartości progowej. np. dla prędkosci wiatru można założyć ws > 0.01. W przeciwnym wypadku funkcja zwróci wartości INF, -INF. Są one wrażliwe na wartości bliskie zeru.

Skłądnia funkcji eval_all jest następująca:

eval_all(data = synop_g, o = obs,m = mod, date = date, ...)
  • data - zestaw danych w pstaci df lub tbl,
  • o - nazwa kolumny w której przechowywane są obserwacje (pomiary),
  • m - nazwa kolumny w której przechowywane są wyniki prognoz,
  • date - nazwa kolumny przechowującej informacje na temat daty w formacie POSIXct.
  • ... - wprowadzamy nazwy kolumn które przechowują zmienne kategoryczne, dzielące nasze obserwacje na podzbiory, np. stacje, bazy danych, źródła danych, czy reprezentujace wyniki z róznych modeli.

Najpier przyjrzymy się naszym danym synop_g.

      date                         MED                MOD           
 Min.   :1984-08-06 00:00:00   Length:90880       Length:90880      
 1st Qu.:1988-05-28 15:45:00   Class :character   Class :character  
 Median :1988-11-02 10:00:00   Mode  :character   Mode  :character  
 Mean   :1990-03-18 22:02:42                                        
 3rd Qu.:1992-09-25 07:15:00                                        
 Max.   :1993-05-19 23:00:00                                        
                                                                    
   station          station_typ           id_var         
 Length:90880       Length:90880       Length:90880      
 Class :character   Class :character   Class :character  
 Mode  :character   Mode  :character   Mode  :character  
                                                         
                                                         
                                                         
                                                         
    param                mod              obs        
 Length:90880       Min.   : 0.000   Min.   : 0.000  
 Class :character   1st Qu.: 2.040   1st Qu.: 1.000  
 Mode  :character   Median : 3.050   Median : 2.400  
                    Mean   : 3.258   Mean   : 2.697  
                    3rd Qu.: 4.170   3rd Qu.: 4.000  
                    Max.   :21.290   Max.   :24.200  
                                     NA's   :4394

Zestaw danych skłąda się z ncol(synop_g) kolumn, które przechowują informacje o dacie (date), bazie danych (MED), stacji, typie stacji, czy o mierzonych parametrach (param) itd..

Wykonamy pierwszy test z grupowaniem po jednej zmiennej

Wykonamy drugi test z grupowaniem po wielu zmiennych

3.2. eval_wd - Model Performance Statistics for wind direction

Funkcja eval_wd oblicza szereg statystycznych parametrów oceny modeli prognostycznych, ale tylko dla kierunku wiatru. Funkja ta oblicza następujace parametry:

  • n - Count - NA
  • mean_o - Mean Observed
  • mean_m - Mean Predicted
  • BIAS - Mean Bias
  • MAE - Gross Error (MAE - Mean Absolut Error)
  • RMSE - Root Mean Square Error
  • SD - Standard Deviation for x = (Predicted value - Observed value)

Funkcje opracowano na podstawie programu metstat-9dec13_1, który jest napisany w fortranie i udostepniany przez RAMBOLL ENVIRON link.

W poniższym raporcie znajdziesz metodologię obliczania wybranych parametrów.

  • Emery, C., Tai, E., Yarwood, G., 2001. Enhanced Meteorological Modeling and Performance Evaluation for Two Texas Ozone Episodes. Env. Int. Corp. 

Skłądnia funkcji eval_wd jest następująca:

eval_wd(data = synop, wdm = wd_m, wdo = wd_o,
                      wsm = ws_m, wso = ws_o,
                      date = date, only_diff = T, ...)
  • data - zestaw danych w pstaci df lub tbl,
  • wdo - nazwa zmiennej w której przechowywane są obserwacje kierunku wiatru (pomiary),
  • wdm - nazwa zmiennej w której przechowywane są wyniki prognoz kierunku wiatru,
  • wso - nazwa zmiennej w której przechowywane są obserwacje prędkości wiatru (pomiary),
  • wsm - nazwa zmiennej w której przechowywane są wyniki prognoz prędkości wiatru,
  • date - nazwa kolumny przechowującej informacje na temat daty w formacie POSIXct.
  • only_diff - Przyjmuje argumnety T lub F, jeśli T - utworzy w zbiorze danych dodatkową zmienną przechowującą informację o różnicy wdm - wdo, w przedziale od -180 do 180.
  • ... - wprowadzamy nazwę zmiennej, która przechowują zmienne kategoryczne, dzielące nasze obserwacje na podzbiory, np. stacje, bazy danych, źródła danych, czy reprezentujace wyniki z róznych modeli.

Najpier przyjrzymy się naszym danym synop.

      date                         MED                MOD           
 Min.   :1984-08-06 00:00:00   Length:90880       Length:90880      
 1st Qu.:1988-05-28 15:45:00   Class :character   Class :character  
 Median :1988-11-02 10:00:00   Mode  :character   Mode  :character  
 Mean   :1990-03-18 22:02:42                                        
 3rd Qu.:1992-09-25 07:15:00                                        
 Max.   :1993-05-19 23:00:00                                        
      ws_m             wd_m            tt_m              ws_o       
 Min.   : 0.000   Min.   :  0.1   Min.   :-27.950   Min.   : 0.000  
 1st Qu.: 2.040   1st Qu.:104.9   1st Qu.:  0.950   1st Qu.: 1.000  
 Median : 3.050   Median :215.8   Median :  9.450   Median : 2.400  
 Mean   : 3.258   Mean   :196.3   Mean   :  9.472   Mean   : 2.697  
 3rd Qu.: 4.170   3rd Qu.:285.5   3rd Qu.: 18.350   3rd Qu.: 4.000  
 Max.   :21.290   Max.   :360.0   Max.   : 38.060   Max.   :24.200  
      wd_o            tt_o          station          station_typ       
 Min.   :  0.5   Min.   :-26.60   Length:90880       Length:90880      
 1st Qu.:137.3   1st Qu.:  2.80   Class :character   Class :character  
 Median :240.0   Median : 10.60   Mode  :character   Mode  :character  
 Mean   :214.9   Mean   : 10.61                                        
 3rd Qu.:295.0   3rd Qu.: 18.49                                        
 Max.   :360.0   Max.   : 36.10                                        
    id_var         
 Length:90880      
 Class :character  
 Mode  :character  
                   
                   
                   
 [ osiągnięto getOption("max.print") -- pominięto 1 wiersz]

Zestaw danych zawiera wymagane zmienne, możeymy wykonać testy. Zauważ, że w tym przypadku potrzebujemy szerkoiej reprezentacji danych. Gdzie zarówno ws i wd są przechowywane w osobnych zmiennych.

Wykonamy pierwszy test z grupowaniem po jednej zmiennej

Wykonamy drugi test z grupowaniem po wielu zmiennych

4. Funkcje oceny modeli (stats plot)

Druga grupa funkcji tworzy graficzne raporty przedstawiajace dokłdność prognoz meteorologicznych. Ananlogicznie jak w poprzednim przypadku zastosowano rozdział na wrzystkie parmetry z wykluczeniem kierunku wiatru eval_met_plot, jak i przygotowano funkcję dedykowaną dla kierunku wiatru eval_wd_plot.

4.1. eval_met_plot - Plot Model Performance Statistics for met. data

Funkcja eval_met_plot tworzy gtaficzny raport oceny modelu dla jednej grupy danych. W funkcji tej nie można stosować grupowania po zmiennej kategorycznej. Można również stosować ją tylko dla jednego parametru, np. prędkości wiatru lub temperatury i innnych.

Funkcja tworzy raport skłądający się z 4 wykresów.

  1. P1 - Observed vs Predicted Scatter Plot - Wykres pierwszy oznaczony jako P1 przedstawia wykres rozrzutu wartości obserwowanych do wartości prognozowanych. W przypadku modeli prgonoz meteorologicznych z reguły dyspinujemy długimi szeregami czasowymi danych związku z tym zamiast punktów prezentuje on sześcioboki do których przypisana jest kolorowa skala opisująca liczność obserwacji znajdujacych się w danym obszarze wykresu.
  2. P2 - Model Performance statistics - jest to graficzne przedstawienie podstawowych statystyk oceny modelu w postaci tabeli. 
  • Name = Unikalna nazwa np. stacji meteorologicznej, lub zestawu danych
  • n = Count - NA
  • d_satrt = Początek okresu pomiarowego
  • d_end = Koniec okresu pomiarowego
  • mean_o = Mean Observed
  • mean_m = Mean Predicted
  • SD_o = Standard Deviation Observed value
  • SD_m = Standard Deviation Predicted value
  • BIAS = Mean Bias (BIAS)
  • MAE = Gross Error (MAE - Mean Absolut Error)
  • RMSE = Root Mean Square Error
  • NMB = Normalised mean bias
  • NME = Normalised mean error
  • FAC2 = Fraction of predictions within a factor of two
  • IOA = The Index of Agreement based on Willmott et al. (2011)
  • R = The Pearson correlation coefficient
  1. P3 - statistical metrics vs. observation range - wykres przedstawijący zmienność niepewności modelu (BIAS, MAE, RMSE) w przedziałach wartości obserowwanej, pozwala określić w których przedziałach wystepują istotne odchylenia względem wartości otrzymanej dla całego zbioru.
  2. P4 - Error histogram przedstawia rozkład reszt diff = (mod - obs) dla każdej wartości.

Składnia funkcji eval_met_plot:

eval_met_plot(data = synop, m = ws_m, o = ws_o, date = date,
          scater_bins = 25, scater_pos_legend = c(0.1, 0.8),
          p3_probs = seq(0,1,0.1), p3_vjust = -3,
          p3_size = 3, p3_pos_legend = c(0.3, 0.9),
          p4_bindwidth = 0.25, p4_pos_sd = 0.15,
          station_name = "MC-LVIA-WRF1 (ws)")
  • data - zestaw danych w postaci df lub tbl,
  • o - zmienna w której przechowywane są obserwacje (pomiary),
  • m - mienna w której przechowywane są wyniki prognoz,
  • date - nazwa kolumny przechowującej informacje na temat daty w formacie POSIXct,
  • scater_bins - numeric vector giving number of bins in both vertical and horizontal directions. Set to 25 by default.
  • scater_pos_legend - dwuelementowy numeryczny wektor c(x,y) w zakresie wartosci od 0 do 1, orientujacy połozenie legendy w obrębie wykresy rozrzutu,
  • p3_probs = sekwencja prawdopodobieństwa dla wyznaczenia kwantyli, działa jak argument probs w funkcji quantile,
  • p3_vjust - liczba rzeczywista, zmiana położenia etykiety względem osi poziomej,
  • p3_size - liczba rzeczywista, rozmiar tekstu etykiet,
  • p3_pos_legend - dwuelementowy numeryczny wektor c(x,y) w zakresie wartosci od 0 do 1, orientujacy połozenie legendy w obrębie wykresy metryk statystycznych
  • p4_bindwidth - liczba całkowiata, szerokość wg której mają zostac wyznaczone przedziały na histogramie,
  • p4_pos_sd - liczba rzeczywista, pozycja względem osi y położenia lini reprezentujacej odchylenie standardowe od średniej,
  • station_name - vektor character przechowujący nazwę stacji, lub inne słowne oznaczenie dla przeprowadzanej analizy

Wykonamy pierwszy test funkcji dla prędkości wiatru mierzonej na wysokości (10 m):

Zapis grafiki do pliku w odpowiednich proporcjach.

Wykonamy drugi test funkcji dla tmperatury przy powierzchni terenu (2m),

Zapis grafiki do pliku w odpowiednich proporcjach.

4.2. eval_wd_plot - Plot Model Performance Statistics for wind direction

Funkcja eval_wd_plot tworzy gtaficzny raport oceny modelu dla jednej grupy danych. W funkcji tej nie można stosować grupowania po zmiennej kategorycznej. Można również stosować ją tylko dla jednego parametru, tj.: kierunku wiatru.

Funkcja tworzy raport skłądający się z 4 wykresów.

  1. P1 - Wind direction error - Wykres pierwszy oznaczony jako P1 przedstawia wykres rozrzutu wartości obserwowanych do wartości prognozowanych (w zakresie od -180, do 180) w układzie biegunowym. Odległość od punktu centralnego oznacza bład predkości czyli (mod_ws - obs_ws). Takie podejście pozwala jednocześnie okreslić, czy występują duże ochylenia prędkości wiatru.

  2. P2 - Model Performance statistics - jest to graficzne przedstawienie podstawowych statystyk oceny modelu w postaci tabeli. 
  • Name = Unikalna nazwa np. stacji meteorologicznej, lub zestawu danych
  • n = Count - NA
  • d_satrt = Początek okresu pomiarowego
  • d_end = Koniec okresu pomiarowego
  • mean_o = Mean Observed
  • mean_m = Mean Predicted
  • BIAS = Mean Bias (BIAS)
  • MAE = Gross Error (MAE - Mean Absolut Error)
  • RMSE = Root Mean Square Error
  1. P3 - statistical metrics vs. observation range - wykres przedstawijący zmienność niepewności modelu (BIAS, MAE, RMSE) w przedziałach wartości obserowwanej, pozwala określić w których przedziałach wystepują istotne odchylenia względem wartości otrzymanej dla całego zbioru.
  2. P4 - Error histogram przedstawia rozkład reszt diff = (mod - obs) dla każdej wartości.

Składnia funkcji eval_wd_plot:

eval_wd_plot(data = synop, wdm = wd_m, wdo = wd_o,
             wsm = ws_m, wso = ws_o, date = date,
             station_name = "LVIA",
             p1_max_wind = 7, p1_size_point = 0.15, p1_degree_pos = 90,
             p3_pos_legend = c(0.1, 0.8), p3_vjust = -3, p3_size = 3,
             p4_bindwidth = 10, p4_pos_sd = 0.0035)
  • data - zestaw danych w postaci df lub tbl,
  • wdo - zmienna w której przechowywane są obserwacje (pomiary) kierunku wiatru,
  • wdm - mienna w której przechowywane są wyniki prognoz kierunku wiatru,
  • wso - zmienna w której przechowywane są obserwacje (pomiary) prędkości wiatru,
  • wsm - mienna w której przechowywane są wyniki prognoz prędkości wiatru,
  • date - nazwa kolumny przechowującej informacje na temat daty w formacie POSIXct,
  • p1_max_wind - maksymalna wartosć etykiety na osi x (różnicy prędkości wiatru),
  • p1_size_point - rozmiar punktów na wykresie p1,
  • p1_degree_pos - położenie etykiet predkości wiatru na osi kierunku wiatru w stopniach w zakresie (-180, 180),
  • p3_vjust - liczba rzeczywista, zmiana położenia etykiety względem osi poziomej,
  • p3_size - liczba rzeczywista, rozmiar tekstu etykiet,
  • p3_pos_legend - dwuelementowy numeryczny wektor c(x,y) w zakresie wartosci od 0 do 1, orientujacy połozenie legendy w obrębie wykresy metryk statystycznych
  • p4_bindwidth - liczba całkowiata, szerokość wg której mają zostac wyznaczone przedziały na histogramie,
  • p4_pos_sd - liczba rzeczywista, pozycja względem osi y położenia lini reprezentujacej odchylenie standardowe od średniej,
  • station_name - vektor character przechowujący nazwę stacji, lub inne słowne oznaczenie dla przeprowadzanej analizy

Wykonamy pierwszy test funkcji:

Zapis grafiki do pliku w odpowiednich proporcjach.

Wykonamy drugi test funkcji:

Zapis grafiki do pliku w odpowiednich proporcjach.

4.3. eval_target_plot - Target Plot

Funkcja eval_target_plot tworzy target plot (wykres celu).

Składnia funkcji eval_wd_plot:

eval_target_plot(data = synop_g, o = obs, m = mod, date = date,
                 title = "Target",  color = MED, shape = station, label = station_typ,
                 size = 4, ...)
  • data - zestaw danych w postaci df lub tbl,
  • o - zmienna w której przechowywane są obserwacje (pomiary),
  • m - mienna w której przechowywane są wyniki prognoz,
  • date - nazwa kolumny przechowującej informacje na temat daty w formacie POSIXct,
  • title - “Target”,
  • color - jak w geom_point(),
  • shape - jak w geom_point(),
  • label - Nazwy zawarte w etykietach,
  • size - jak w geom_point,
  • ... - zmienne kategoryczne dzielace zestaw danych na podzbiory

Wykonamy pierwszy test funkcji:

Zapis grafiki do pliku w odpowiednich proporcjach.

Wykonamy drugi test funkcji:

Zapis grafiki do pliku w odpowiednich proporcjach.

Wykonamy trzeci test funkcji:

Zapis grafiki do pliku w odpowiednich proporcjach.

Mateusz Rzeszuteka

2019-08-05