Wprowadzenie

Dane dotyczące warunków meteorologicznych bardzo często analizowane są z wykorzystaniem informacji o jakości powietrza, jako że obie te grupy danych są ze sobą wysoko skorelowane. Dane pomiarowe obu tych grup mogą służyć do tworzenia modeli przewidujących wartości stężeń zanieczyszczeń bądź innych powiązanych parametrów. Jednocześnie coraz powszechniejszą praktyką staje się wykorzystywanie do tego metod uczenia maszynowego, które pozwalają na uzyskanie wysokiej jakości wyników modelu statystycznego (w porównaniu np. do prostej regresji).

Stacje meteorologiczne i stacje jakości powietrza zarządzane są w Polsce przez dwie różne instytucje - odpowiednio przez IMGW oraz GIOŚ. Powoduje to, że punkty pomiarowe znajdują się często w odmiennych miejscach. Idealnie byłoby, gdyby zarówno parametry meteorologiczne jak i te dotyczące jakości powietrza były realizowane w jednym i tym samym punkcie w przestrzeni.

Najpowszechniejszą praktyką w tego typu analizach jest zatem wybieranie jak najbardziej reprezentatywnych danych dotyczących jakości powietrza – najprostszym kryterium jest odległość stacji meteorologicznej od stacji jakości powietrza. Można przyjąć, że im ta odległość jest mniejsza, tym łatwiej unikniemy błędu (co niestety nie jest regułą).

1. Dane dotyczące jakości powietrza w Polsce

Podmiotem utrzymującym i zarządzającym pomiarami jakości powietrza w Polsce (zgodnie z wytycznymi Unii Europejskiej) jest Główny Inspektorat Ochrony Środowiska oraz podległe mu Wojewódzkie Inspektoraty. Od kilku lat na stronie GIOŚ - Portal Jakości Powietrza archiwizowane i udostępniane są dane dotyczące m. in. stężeń zanieczyszczeń oraz wyniki oceny jakości powietrza (rocznej, pięcioletniej). To właśnie z tych danych będziemy dzisiaj korzystać, można je pobrać z banku danych pomiarowych GIOŚ.

Dane te niestety nie są przechowywane w zunifikowanym formacie, plik z każdego roku różni się nieco od poprzedniego, co utrudnia wczytywanie tych danych np. do R.

Skupimy się na danych pomiarowych ze stacji:

  • automatycznych,
  • kontenerowych stacjonarnych,
  • tła miejskiego/regionalnego,
  • aktywnych, dla których istnieją pomiary od początku 2015 roku.

Sprawdźmy jak wygląda mapa aktywnych stacji jakości powietrza w Polsce. Aby ułatwić wybór pojedynczej stacji, można zawęzić zbiór o wymienione wyżej warunki (funkcja filter).

# Metadane stacji GIOŚ ####
library(openxlsx)
# Pobranie i wczytanie metadanych dotyczących stacji
setwd("E:/Uczelnia/PNoZII/Projekt_2")
download.file("http://powietrze.gios.gov.pl/pjp/archives/downloadFile/305", destfile="metadane.zip",
              mode = "wb", method = "wininet")
gios_inv <- read.xlsx(xlsxFile = unzip("metadane.zip"))

# Usunięcie polskich znaków z nazw kolumn, których będziemy używać w kodzie
colnames(gios_inv)[c(7,13)] <- c("Data.zamkniecia", "Miejscowosc")
## Konwersja daty Excela
gios_inv$Data.uruchomienia <- convertToDate(gios_inv$Data.uruchomienia, format="%Y-%m-%d")
gios_inv$Data.zamkniecia <- convertToDate(gios_inv$Data.zamkniecia, format="%Y-%m-%d")
## Zmiana nazw kolumn, aby nie zawierały znaków specjalnych
colnames(gios_inv)[15:16] <- c("Szerokosc", "Dlugosc")

library(tidyverse)
# Selekcja stacji kontenerowych, odrzucenie stacji o niepoprawnej lokalizacji geograficznej
gios_inv %>% 
  filter(Szerokosc != -999 & Dlugosc != -999 & Rodzaj.stacji == "kontenerowa stacjonarna") -> gios_inv

# Tworzenie mapki stacji w pakiecie Leaflet ####
library(leaflet)
## Deklarowanie zawartości popupu
desc <- paste(paste(
  gios_inv$Kod.stacji,
  paste("Miejscowość:", gios_inv$Miejscowosc),
  paste("Data uruchomienia:", gios_inv$Data.uruchomienia),
  paste("Data zamknięcia:", gios_inv$Data.zamkniecia),
  paste("Typ stacji:", gios_inv$Typ.stacji),
  paste("Typ obszaru:", gios_inv$Typ.obszaru),
  paste("Współrzędne:", gios_inv$Dlugosc,"E, ", gios_inv$Szerokosc,"N"),
  sep = "<br/>"
))

## Deklarowanie palety legendy dla typów stacji
pal <- colorFactor(c("blue", "red", "darkgreen"), domain = c("tło", "komunikacyjna", "przemysłowa"))

## Tworzenie mapki
leaflet() %>%
  addTiles() %>%
  addCircleMarkers(data=gios_inv,
                   lng= ~ Dlugosc,
                   lat= ~ Szerokosc,
                   popup = desc,
                   fillOpacity = 1,
                   radius = 3,
                   color = ~ pal(Typ.stacji)) %>%
  addLegend('bottomright', colors = c("green", "red", "blue"), labels = c("tło", "przemysłowa", "komunikacyjna"),
            opacity = 1)

Do dalszych analiz wybieram stację o kodzie SlDabro1000L. Sprawdźmy metadane:

Pobieranie danych o jakości powietrza

Ze względu na niejednorodność konstrukcji plików przechowujących informacje o stężeniach w serwisie GIOŚ, wczytywanie ich do R nie jest takie łatwe. W poniższych krokach, do pobrania całej bazy danych (lata 2000-2018) i jej przygotowania wykorzystamy funkcje przygotowane i udostępnione dzięki uprzejmości Pana dr. Mateusza Rzeszutka. Pozwoli nam to na zaoszczędzenie czasu i skupienie się na właściwym celu projektu.

Skompresowane pliki z kolejnymi latami danych są dostępne na stronie GIOŚ pod poniższymi linkami:

Poniższa 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.

Wykonujemy test funkcji w następujący sposób, pobierając dane dla roku 2018:

Powyższą funkcję możemy wykorzystać do pobrania całej bazy naraz. Uwaga, cała baza zajmuje 630MB, więc należy upewnić się, że takim miejscem dysponujemy w folderze roboczym.

Zanim przejdziemy dalej musimy zdawać sobie sprawę z kilku problemów związanych z plikową bazą danych GIOŚ.

  • pliki *.xlsx przechowujące dane o pomiarach stężeń nie są jednorodne dla każdego roku.
  • pliki *.xlsx podzielone są na stężenia 24-godzinne i 1-godzinne, co wymaga stosowania różnego formatowania daty.
  • we wskazanym okresie zmienił się system nazwy stacji. W związku z tym w różnych latach będą stosowane różne oznaczenia stacji monitoringu jakości powietrza. Oznaczenia stacji zmieniły się począwszy od 2015 r. Ten problem można rozwiązać tylko w połączeniu z importem metadanych. Ten problem dotyczy tylko sytuacji, w której chcemy otrzymać szereg czasowy danych z wielolecia. Została do tego celu opracowana dodatkowa funkcja.
  • Każdy folder zawiera różną liczbę plików.
  • Każdy plik przechowuje informacje o stężeniach wybranego zanieczyszczenia powietrza dla różnej liczby stacji.

Przetwarzanie danych o jakości powietrza

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

Wprowadzono dodatkowo argument sciezka. Jeśli pobrana baza danych znajduje się w innej lokalizacji, to można ją podać jako zmienna character. Jeśli plikowa baza danych znajduje się w katalogu roboczym projektu RStudio, to pozostawiamy wartość domyślną.

# Wczytujemy do R
names(pliki_all) <- paste0("R",linki[,2])

gios_read <- function(nazwa, czas_mu = "24g", sciezka= "") {
    
    lok <- getwd()
    setwd(dir = paste0(sciezka,substr(nazwa, 1, 4), "/"))
  
  # Rózne formaty plików xlsx ustawienia parametrów wczytywania 
  
  if(str_sub(nazwa, 1,4) %in% c(2016, 2017, 2018)) {
    
    startRow = 2 ; end_row = 4 
    
  } else if(str_sub(nazwa, 1,4) %in% c(2000:2015)) {
    
    startRow = 1 ; end_row = 2 
  }
  
  # Wczytywanie danych
  
  dane <- openxlsx::read.xlsx(nazwa,           # Wczytujemy dane z exele
                              startRow = startRow, 
                              colNames = T)  
  colnames(dane)[1] <- "date"                  # modyfikacja nazwy pierwszej kolumny
  sub =  dane[1,2]                             # nazwa substancji
  
  dane <- map_df(.x = dane[-c(1:end_row),],    # usunięcie nagłówka oraz
                 .f = str_replace,             # zamiena separatora z ',' na '.'
                 pattern = ",",                # w kazdej kolumnie
                 replacement = ".") 
  
  # konwersja daty z formatu liczbowego exela
  
  if (czas_mu == "1g") {
    dane <- dane %>%
      mutate(date = as.numeric(date),
             date = janitor::excel_numeric_to_date(date, include_time = T,
                                          date_system = "modern",
                                          tz = "UTC") - 3600,
             date = lubridate::round_date(date, "hour")) # zaokrąglanie daty, problem dokładności daty w excel
    
    
  } else if(czas_mu == "24g") {
    
    dane$date <- as.Date(as.numeric(dane$date), 
                         origin = "1899-12-30", 
                         tz = "UTC")
  }
  
  # ostatni szlift
  
  dane <- dane %>% 
    mutate_if(is.character,             # Zamiena typ danych na double 
              as.double)  %>% 
    gather(key = "kod",                 # Wąski układ danych (wygodny do łączenia tych danych)    
           value = "obs", -date) %>% 
    mutate(sub = sub, 
           obs = round(obs, 6)) %>%     # Nazwa substancji
    .[,c("kod", "sub", "date", "obs")]
  
  setwd(lok)
  
  return(dane)
  
}

Wynikiem jest obiekt tibble, który zawiera dane gotowe do przetwarzania. Wszystkie daty są w UTC+01, co odpowiada czasowi zimowemu w Polsce.

Wykonajmy test dla 1-godzinnych danych, które będą nam potrzebne w niniejszym projekcie.

I na koniec, znajdźmy tylko te obserwacje, które pasują do naszej wybranej stacji, czyli SlDabro1000L. Niestety, ten kod stacji obowiązuje dopiero od roku 2015, stary kod stacji to SlDabroDabr_1000L. Żeby wyciągnąć z naszych danych odpowiednie pomiary, musimy ten fakt uwzględnić:

Łączenie danych o jakości powietrza z meteorologicznymi

Sprawdziliśmy nieco wcześniej lokalizację geograficzną stacji jakości powietrza w Dąbrowie Górniczej. Wykorzystajmy pakiet worldmet do pobrania metadanych stacji meteorologicznych znajdujących się w jej pobliżu. Następnie, pobierzemy dane meteorologiczne dla najbliższej dostępnej stacji dla dostępnych lat z przedziału 2005-2018 (tam, gdzie dostępne są pomiary jakości powietrza). W moim przypadku będzie to stacja w Katowicach, oddalona o 17.7km od stacji jakości powietrza w Dąbrowie Górniczej (mogło być lepiej).

Następnie łączymy pomiary meteorologiczne z danymi o stężeniach PM10. Pamiętajmy, że dane są w różnych strefach czasowych, dlatego przed połączeniem należy “dodać” godzinę do danych meteorologicznych.

2. Prognoza stężenia PM10 na 2018 rok

Wykonamy ćwiczenie polegające na zaprognozowaniu wartości stężenia pyłu PM10 w zależności od warunków meteorologicznych na danym obszarze. Wykorzystamy do tego las losowy (random forest), będący jedną z metod ensemblingu uczenia maszynowego. Jest to rozszerzenie metody drzew decyzyjnych, która polega na budowaniu modelu w oparciu o zadawanie i dzielenie zbioru według warunków nakładanych na dane zmienne. I chociaż lasy losowe nadają się lepiej do prognozowania danych dyskretnych bądź nominalnych (metoda klasyfikacyjna), to również w naszym przypadku powinna przynieść satysfakcjonujące efekty.

W teorii powinniśmy poświęcić więcej czasu na dobór metody prognozowania i sprawdzenia rzeczywistych regresji między skorelowanymi parametrami naszego zbioru, ale na obecnym etapie nauki samo dostrzeżenie potencjału uczenia maszynowego w dziedzinie nauk o środowisku jest bardzo wartościowe.

Wstępna analiza danych

Braki danych i selekcja

W metodach uczenia maszynowego wymagane jest posiadanie kompletnych wierszy obserwacji, a także tylko tych zmiennych, które wniosą informację dodaną do naszego modelu.

Wykreślmy już znany nam wykres obrazujący ilość rekordów w danym roku:

Podzielmy zbiór na dwa podzbiory. Jeden z nich, tzw. “zbiór treningowy”, będzie obejmować dane, na których nauczymy model przewidywać stężenia PM10 na podstawie danych meteorologicznych. Nie wybierajmy zbyt wiele rekordów, bo zwiększy to wykładniczo czas obliczeniowy. Zdecydujmy się na 2-3 lata najbardziej kompletnych danych. Ja do zbioru treningowego wybieram lata 2016-2017. Do zbioru testowego, czyli takiego, na którym będziemy weryfikować nasz model, wybierzmy najnowszy rok, 2018.

Korelacje między zmiennymi

Sprawdźmy proste korelacje liniowe między zmiennymi. Aby mówić w ogóle o wykorzystaniu modelu, musimy mieć solidne podstawy, aby sądzić, że związek między zmiennymi jest silny i jedną z nich (PM10 - zmienna objaśniana) da się przedstawić jako funkcja pozostałych zmiennych (zmienne objaśniające).

Wykorzystamy znany nam już corrplot.

Jak widzimy, korelacje liniowe Pearsona między stężeniem pyłu PM10 a innymi zmiennymi są jedynie na umiarkowanym poziomie, ale wszystkie są statystycznie istotne. To w zupełności wystarczające, aby zbudować na tym model. Nie przejmujmy się, jeżeli jakaś korelacja okaże się nieistotna statystycznie - model prawdopodobnie sam ją odrzuci na etapie uczenia.

Wykonanie predykcji

Do uruchomienia lasu losowego wykorzystamy pakiet caret, który jest bardzo bogatym pakietem do modelowania statystycznego w R. Aby oszczędzić nieco czasu uruchomimy także obliczenia równoległe, wykorzystując dostępne rdzenie procesora.

## Random Forest 
## 
## 11184 samples
##     8 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 8947, 8946, 8948, 8947, 8948 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   1     22.46841  0.7965890  13.06427
##   2     21.35689  0.8095574  12.29682
##   3     21.03905  0.8136887  12.16124
##   4     20.90376  0.8154691  12.07384
##   5     20.83905  0.8162529  12.07695
##   6     21.03553  0.8125722  12.12584
##   7     21.06461  0.8117339  12.14637
##   8     21.20962  0.8091120  12.22436
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 5.

Testowanie modelu

Teraz, kiedy już mamy model, możemy sprawdzić jak zachowuje się on na danych, które nie były wykorzystane do jego uczenia. Wykorzystajmy utworzony model do oszacowania stężenia pyłu PM10 w tym roku na zbiorze testowym, czyli danych dla 2018 roku.

Żeby sprawdzić, jak dobrze nasz model odwzorowuje rzeczywistość musimy mieć jakąś wartość odniesienia. Sporządźmy dwa wykresy obrazujące nam zależność wartości stężenia PM10 zaobserwowanego na stacji jakości powietrza od zamodelowanej jego wartości. Pierwszy z nich to wykres rozrzutu.

Wykres rozrzutu pokazuje jak mocno wyniki prognoz odchylone są od pomiarów stężeń PM10. Szara przerywana linia obrazuje model idealny (w sytuacji, gdy prognoza równa jest pomiarowi, czyli x=y), a niebieska – linię trendu uzyskaną z naszych danych, o równaniu wyświetlonym na wykresie.

Według naszych danych, prognoza jest zaniżona w stosunku do pomiarów stężenia PM10 na stacji w Dąbrowie Górniczej, co jest szczególnie widoczne w przypadku najwyższych stężeń. Może to wynikać z faktu, że warunki meteorologiczne nie są jedynymi czynnikami wpływającymi na stężenie zanieczyszczenia na danym obszarze. Pominęliśmy całkowicie kwestię emisji, która, w przypadku sytuacji wyjątkowych i chwilowych, może drastycznie wpłynąć na mierzone stężenia zanieczyszczeń.

Wykres trendu prognoz (linia czerwona) od obserwacji stężeń PM10 (linia szara) pozwala na ogólną ocenę zgodności tych dwóch serii danych. Możemy zauważyć, że wniosek wyciągnięty na podstawie poprzedniego wykresu jest potwierdzony – najwyższe wartości stężeń są zaniżane w naszym modelu. Ogólny trend jest zachowany, chociaż również dla najniższych stężeń zdarzają się odchylenia.

3. Projekt do wykonania

  1. Wybierz stację jakości powietrza na terenie Polski, zlokalizuj ją na mapie i opisz jej położenie wspomagając się metadanymi z bazy GIOŚ.
  2. Pobierz dane z banku danych pomiarowych GIOŚ i wczytaj je do R za pomocą dostarczonych funkcji.
  3. Sprawdź, za pomocą pakietu worldmet, bliskość stacji meteorologicznych w bazie ISD NOAA, dla których istnieją pomiary w analizowanym okresie. Pobierz dane dla odpowiedniego okresu i najbliższej stacji meteorologicznej.
  4. Sprawdź czy istnieją istotne zależności liniowe między parametrami (za pomocą współczynnika korelacji Pearsona).
  5. Sporządź według instrukcji w skrypcie prognozę stężenia wybranego zanieczyszczenia powietrza z wykorzystaniem lasu losowego. Podaj wartość parametru mtry oraz błędu RMSE dla finalnego modelu oraz przygotuj dwa wykresy:
  1. rozrzutu wartości prognozowanych od obserwowanych wraz z równaniem regresji i linią trendu,
  2. trendu zmian prognozy i obserwacji stężeń dla wybranego okresu (tak, aby linie i wartości były dobrze widoczne).
  1. Skomentuj uzyskane wyniki i sporządź odpowiednie wnioski.

Gotowy projekt w formacie *.pdf wraz ze skryptem *.R lub w formacie *.html\*.pdf utworzonym z wykorzystaniem R Markdown (z widocznym kodem R) należy przesłać w w ciągu dwóch tygodni od przeprowadzonych zajęć na adres mailowy prowadzącego.