Autorem jest studentka pierwszego roku na Akademii Górniczo-Hutniczej im. Stanisława Staszica w Krakowie

Wydział: Geodezji Górniczej i Inżynierii Środowiska

Kierunek: Geoinformacja

E-mail:

Wybranie stacji jakości powietrza na terenie Polski do analizy

Poniżej zostały wypisane wszystkie pakiety potrzebne do wykonania tego projektu:

Stworzono mapkę interaktywną aktywnych stacji w pakiecie Leaflet w oparciu o wymienione poniżej warunki. Stacje mają być/mieć:

  • automatyczne,
  • kontenerowe stacjonarne,
  • tło miejskie/regionalne,
  • aktywne (pomiary istnieją od początku 2015 roku).

Poniższy kod przedstawia sposób pobrania i wczytania metadanych stacji. Jako, że RStudio nie obsługuje polskich znaków, usunięto je z nazw kolumn. Wyselekcjonowano stacje kontenerowe i odrzucono te o niepoprawnej lokalizacji geograficznej. Wybrano kolory, jakie będzie miała legenda mapy.

Wybrano stację znajdującą się nad Morzem Bałtyckim, dokładniej mówiąc – nad Zatoką Gdańską, w mieście Gdańsk. Stacja ma kod PmGdaPoWie01. Za pomocą poniższego kodu sprawdzono jej metadane:

Nr Kod.stacji Kod.międzynarodowy Nazwa.stacji Stary.Kod.stacji Data.uruchomienia Data.zamkniecia Status Typ.stacji Typ.obszaru Rodzaj.stacji Województwo Miejscowosc Ulica Szerokosc Dlugosc
744 PmGdaPoWie01 PL0045A AM1 Gdańsk Śródmieście Pm.a01a 1996-10-01 NA aktywny tło miejski kontenerowa stacjonarna POMORSKIE Gdańsk ul. Powstańców Warszawskich 54.35334 18.63528

Pobranie danych dotyczących jakości powietrza

Do pobrania całej bazy danych (lata 2000-2018) i jej przygotowania wykorzystano poniższe funkcje. Najpierw pobrano linki ze strony GIOŚ, w których znajdują się pliki z poszczególnymi latami danych:

Pierwsza funkcja pobiera archiwum danych z portalu GIOŚ, rozpakowuje je do wskazanego folderu (argument rok) oraz zwraca wektor zawierający wykaz plików znajdujących się w danym archiwum.

Funkcja została przetestowana przez pobranie danych dla roku 2018

Przetwarzanie danych o jakości powietrza

Kolejna funkcja zakłada, że pozyskane wcześniej dane znajdują się w folderze roboczym. Pliki dla każdego roku znajdują się w osobnym folderze, którego nazwą jest kolejny rok.

Najpierw wczytano pliki do R, a następnie użyto drugiej funkcji, aby dane mogły być gotowe do przetwarzania. Wynikiem jest obiekt typu tibble.

Wykonano test dla 1-godzinnych danych, które będą potrzebne przy dalszej analizie. Wyszukano plików z pomiarami PM10 co godzinę, a następnie wczytano ich zawartość do jednej zmiennej.

Użyto funkcji filter w celu znalezienia obserwacji pasujących do wybranej stacji. Jako że miała ona kiedyś inną nazwę trzeba uwzględnić też tą już nieobowiązującą:

Łączenie danych o jakości powietrza z danymi meteorologicznymi

Za pomocą pakietu worldmet znalezino stację meterologiczną znajdującą się najbliżej wybranej stacji jakości powietrza. Następnie pobrano dane meterologiczne z bazy ISD NOAA dla lat z przedziału 2005-2018. W tym przypadku najbliższa stacja znajdowała się w pobliżu Portu Lotniczego Gdańsk im. Lecha Wałęsy, oddalona od stacji jakości powietrza o 11.3 km.

## # A tibble: 10 x 15
##    USAF  WBAN  STATION CTRY  ST    CALL  latitude longitude `ELEV(M)` BEGIN     
##    <chr> <chr> <chr>   <chr> <chr> <chr>    <dbl>     <dbl>     <dbl> <date>    
##  1 1215~ 99999 LECHA ~ PL    <NA>  EPGD      54.4      18.5     149   1931-01-01
##  2 1214~ 99999 PRUSZC~ PL    <NA>  EPPR      54.2      18.7       6.4 2002-06-01
##  3 1215~ 99999 GDANSK~ PL    <NA>  <NA>      54.3      18.9       9   1982-09-01
##  4 1214~ 99999 OKSYWIE PL    <NA>  <NA>      54.6      18.5      70   2002-06-01
##  5 1213~ 99999 HEL     PL    <NA>  <NA>      54.6      18.8       3   1955-07-01
##  6 1237~ 99999 MALBORK PL    <NA>  EPMB      54.0      19.1       4.9 2005-01-01
##  7 1213~ 99999 LEBUNIA PL    <NA>  EPCE      54.4      17.8     151   2002-06-01
##  8 1216~ 99999 ELBLAG~ PL    <NA>  EPEL      54.2      19.5      43   1933-01-02
##  9 1212~ 99999 LEBORK  PL    <NA>  <NA>      54.6      17.8      41   1953-04-01
## 10 1212~ 99999 LEBA    PL    <NA>  <NA>      54.8      17.5       6   1973-01-01
## # ... with 5 more variables: END <date>, code <chr>, longR <dbl>, latR <dbl>,
## #   dist <dbl>

Kolejnym krokiem było połączenie pomiarów meteorologicznych z danymi o stężeniach PM10 dla wybranych lat. Aby czas w obu danych się zgadzał trzeba dodać godzinę w danych meteorologicznych.

Prognoza stężenia PM10 na 2018 rok

Wykorzystano las losowy (random forest) w celu zaprogramowania wartości stężenia pyłu PM10 w zależności od warunków meteorologicznych na danym obszarze.

Braki danych i selekcja

Pierwszym krokiem było sprawdzenie, czy nie ma braków danych, gdyż w tej metodzie wymagane jest posiadanie kompletnych wierszy obserwacji oraz zmiennych, które sa istotne dla tworzonego modelu. Jako, że kolumny visibility oraz atmos_pres posiadały znaczne braki danych, zdecydowano się je usunąć. Następnie za pomocą funkcji complete.cases() wybrano tylko te wiersze, które mają kompletne dane.

Wykres obrazujący ilość rekordów w danym roku

Utworzono wykres, który pokazuje, ile pomiarów było w poszczególnych latach od 2015 do 2018.

Jak widać, ilość rekordów pomiarowych wacha się od około 3700 do 5200. Podzielono zbiór na dwa podzbiory:

  • treningowy - obejmuje dane, na których nauczono model przewidywać stężenia PM10 na podstawie danych meteorologicznych,
  • testowy - służy do weryfikowania modelu.

Korelacje między zmiennymi

Macierz korelacji liniowej Pearsona służy do pokazania zależności między parametrami. Współczynnik korelacji może przybierać wartości w zakresie od -1 do 1. W celu zwiększenia czytelności zmieniono nazwy kolumn na polskie odpowiedniki.

Można zauważyć, że wszystkie korelacje między stężeniem pyłu PM10 a innymi zmiennymi są statystycznie istotne. Jednak żaden związek się nie wyróżnia, wszystkie są na umiarkowanym poziomie. Najwyższa korelacja pułu PM10 jest z podstawą chmur.

Wykonanie predykcji

Do uruchomienia lasu losowego wykorzystano pakiet caret. Ustawiono metodę tzw. “kroswalidacji”, czyli samo-testowania się modelu w trakcie jego budowy. W modelu wykorzystano wszystkie zmienne oprócz tej objaśnianej - PM10. Uruchomiono także klaster obliczeniowy z wszystkich dostępnych rdzeni.

## Random Forest 
## 
## 9804 samples
##    6 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 7843, 7843, 7844, 7843, 7843 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   1     12.67305  0.5689482  8.312695
##   2     12.32689  0.5857032  7.978540
##   3     12.32706  0.5853855  7.977057
##   4     12.42131  0.5793072  8.005592
##   5     12.43681  0.5787013  8.031056
##   6     12.52416  0.5730743  8.100871
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 2.

Testowanie modelu

Za pomocą poniższego kodu sprawdzono, jak zachowuje się model na danych, które nie były wykorzystane do jego uczenia. Wykorzystano go do oszacowania stężenia pyłu PM10 na zbiorze testowym.

Wykres rozrzutu

Utworzono dwa wykresy obrazujące zależność wartości stężenia PM10 zaobserwowanego na stacji jakości powietrza od zamodelowanej jego wartości.

Najpierw sporządzono wykres rozrzutu, który pokazuje, jak mocno odchylone są wyniki prognoz od pomiarów stężeń PM10.

Przerywana linia obrazuje model idealny, gdzie prognoza jest równa pomiarowi, a ciągła pokazuje linię trendu uzyskaną z badanych danych, o równaniu widocznym na wykresie.

Można łatwo zauważyć, że prognoza jest zaniżona w stosunku do pomiarów stężenia PM10 na stacji w Gdańsku, co jest szczególnie widoczne w przypadku najwyższych stężeń.

Wykres trendu

Ostatnim wykresem jest wykres trendu, który pozwala na ogólną ocenę zgodności prognoz z obserwacjami stężeń PM10.

Prognozy na wykresie są zakreślone linią czerwoną, natomiast obserwacje stężeń PM10 - linią szarą. Można potwierdzić wniosek wysunięty na podstawie poprzedniego wykresu - najwyższe wartości stężeń są zaniżane, Można też zauważyć, że najniższe wartości są nieco zawyżane. Jednak można stwierdzić, że ogólny trend jest zachowany.