Łączenie danych meteorologicznych z danymi jakości powietrza z bazy GIOŚ

Projekt 2.

---

Wybranie stacji jakości powietrza

Dane dotyczące stacji zlokalizowanych na terenie Polski pobrano (z banku danych pomiarowych GIOŚ) i wczytano. Następnie uporządkowano je i dokonano odpowiedniej selekcji.

Spis stacji znajduje się w tabeli poniżej.

Przy pomocy pakietu Leaflet stworzono mapę stacji z podziałem na ich typy.

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)

Wybrana do analizy stacja znajduje się w Zdzieszowicach w województwie opolskim.

Poniżej znajdują się szczegółowe informacje na jej temat.

---

Przetwarzanie danych

Pobranie danych o jakości powietrza

Ze strony GIOŚ pobrano pliki z kolejnymi latami danych.

Funkcja gios_download 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.

gios_download <- function(url, rok = "2018", ...) {
  
  if (!is.character(url) | !is.character(rok)) {
    stopifnot("Wartości argumentów: `url` i `rok` muszą być typu character") }
  
  if (length(url) != 1 | length(rok) != 1) {
    stopifnot("Argumenty: `url` i `rok` przyjmują tylko pojedyncze wartości") }

  zip <- paste0(rok,".zip")
  
  download.file(url, zip, quiet = T, mode = "wb", method = "wininet")

  unzip(zipfile = zip, exdir = rok)
  file.remove(zip)

  pliki <- dir(paste0(rok, "/"))
  
  pliki <- pliki[!str_detect(pliki, "Jony")]
  pliki <- pliki[!str_detect(pliki, "epozycj")]
  
  return(pliki) 
}

Wykonano test funkcji dla roku 2018.

pliki_2018 <- gios_download(url = "http://powietrze.gios.gov.pl/pjp/archives/downloadFile/303", rok = "2018")
knitr::kable(pliki_2018)

x
2018_As(PM10)_24g.xlsx
2018_BaA(PM10)_24g.xlsx
2018_BaP(PM10)_24g.xlsx
2018_BbF(PM10)_24g.xlsx
2018_BjF(PM10)_24g.xlsx
2018_BkF(PM10)_24g.xlsx
2018_C6H6_1g.xlsx
2018_C6H6_24g.xlsx
2018_Cd(PM10)_24g.xlsx
2018_CO_1g.xlsx
2018_DBahA(PM10)_24g.xlsx
2018_formaldehyd_24g.xlsx
2018_Hg(TGM)_1g.xlsx
2018_Hg(TGM)_24g.xlsx
2018_IP(PM10)_24g.xlsx
2018_Ni(PM10)_24g.xlsx
2018_NO2_1g.xlsx
2018_NO2_24g.xlsx
2018_NOx_1g.xlsx
2018_O3_1g.xlsx
2018_Pb(PM10)_24g.xlsx
2018_PM10_1g.xlsx
2018_PM10_24g.xlsx
2018_PM25_1g.xlsx
2018_PM25_24g.xlsx
2018_SO2_1g.xlsx
2018_SO2_24g.xlsx

Następnie pobrano całą bazę i zapisano ją jako pliki_all.

---

Skorzystano z funkcji, której zadaniem jest stworzenie obiektu typu tibble, zawierającego dane gotowe do przetwarzania. Użyto pozyskane wcześniej dane znajdujące się w folderze roboczym. Zapewniając, że pliki dla każdego roku zostały zlokalizowane w osobnym folderze, którego nazwą jest kolejny rok.

Wszystkie daty są w UTC+01, co odpowiada czasowi zimowemu w Polsce.

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), "/"))
  
  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 
  }
  
  dane <- openxlsx::read.xlsx(nazwa,          
                              startRow = startRow, 
                              colNames = T)  
  colnames(dane)[1] <- "date"                  
  sub =  dane[1,2]                             
  
  dane <- map_df(.x = dane[-c(1:end_row),],    
                 .f = str_replace,             
                 pattern = ",",                
                 replacement = ".") 
  
  
  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")) 
    
    
  } else if(czas_mu == "24g") {
    
    dane$date <- as.Date(as.numeric(dane$date), 
                         origin = "1899-12-30", 
                         tz = "UTC")
  }
  
  dane <- dane %>% 
    mutate_if(is.character,             
              as.double)  %>% 
    gather(key = "kod",                    
           value = "obs", -date) %>% 
    mutate(sub = sub, 
           obs = round(obs, 6)) %>%     
    .[,c("kod", "sub", "date", "obs")]
  
  setwd(lok)
  
  return(dane)
  
}

Przygotowanie pomiarów wybranej wcześniej stacji.

PM10_1h %>%
  filter(kod=="OpZdziePiast") -> stacja_PM10

Łączenie danych o jakości powietrza z meteorologicznymi

Pobrano metadane stacji meteorologicznych znajdujących się w pobliżu Zdziszowic.

library(worldmet)
noaa_isd <- getMeta(end.year="current", lon=18.12074, lat=50.42353, returnMap=F)

Do dalszej analizy wybrano stację oddaloną o 24.1 km, która znajduje się w Opolu.

opole_met <- importNOAA(code="125300-99999", year=seq(2005,2018, by=1))

Połączono pomiary meteorologiczne z danymi o stężeniach PM10.

opole_met$date <- opole_met$date + 3600
dane <- inner_join(opole_met, stacja_PM10)

---

Prognoza stężenia PM10 na 2018 rok

Poniższa analiza będzie polegać na zaprognozowaniu wartości stężenia pyłu PM10 w zależności od warunków meteorologicznych na danym obszarze. W tym celu zostanie wykorzystany las losowy (random forest).

dane_rf <- dane[,c(1, 9:16,26)]
dane_rf <- dane_rf[complete.cases(dane_rf),]

Podzielono zbiór na dwa podzbiory. Pierwszy obejmujący dane, które zostaną przeznaczone do nauki modelu przewidywania stężenia PM10 na podstawie danych meteorologicznych (lata 2016-2017). Drugi - testowy, jego celem jest weryfikacja modelu (rok 2018).

dane_rf %>%
  filter(years>2015 & years<2018) %>%
  select(-date, -years) -> dane_rf_train

dane_rf %>%
  filter(years==2018) %>%
  select(-years) -> dane_rf_test

Korelacje między zmiennymi

Następnym etapem jest sprawdzenie prostych korelacji liniowych między zmiennymi.

Korelacje między stężeniem pyłu PM10 a innymi zmiennymi są na umiarkowanym poziomie.

Wykonanie predykcji

library(caret)
library(parallel)
library(doParallel)

cross.walid <- trainControl(method = "cv",
                            number = 5,
                            allowParallel = T,
                            returnResamp = 'all')

tunegrid <- expand.grid(.mtry = c(1:8))

cluster <- makeCluster(detectCores() - 1)
registerDoParallel(cluster)

model_rf <- train(obs ~ ., 
                  data = dane_rf_train, 
                  method = "rf", 
                  replace = T,
                  ntree = 100,
                  metric = 'RMSE',
                  tuneGrid = tunegrid, 
                  trControl = cross.walid)

stopCluster(cluster)
registerDoSEQ() 

model_rf

## Random Forest 
## 
## 3100 samples
##    8 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 2480, 2480, 2480, 2480, 2480 
## Resampling results across tuning parameters:
## 
##   mtry  RMSE      Rsquared   MAE     
##   1     20.65952  0.7158069  12.62585
##   2     19.50344  0.7393187  11.47035
##   3     19.35469  0.7414440  11.25426
##   4     19.36096  0.7408737  11.14210
##   5     19.27020  0.7432377  11.09962
##   6     19.41861  0.7387918  11.15198
##   7     19.49188  0.7365411  11.23137
##   8     19.84891  0.7264744  11.37127
## 
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 5.

Testowanie modelu

dane_rf_test$mod <- predict(model_rf, dane_rf_test)

Sporządzono wykres rozrzutu obrazujący zależność wartości stężenia PM10 zaobserwowanego na stacji jakości powietrza od zamodelowanej jego wartości. Na jego podstawie możliwe jest sprawdzenie, jak dobrze model odwzorowuje rzeczywistość.

library(lubridate)
library(ggpmisc)

ggplot(data=dane_rf_test, aes(x=obs, y=mod))+
  geom_point(alpha=0.8, color="#CE9DDA") +
  geom_smooth(method="lm", formula = y~x-1, colour = "#1FA2C3") +
  geom_abline(intercept=0, col="grey", linetype="dashed", size=1) +
  stat_poly_eq(aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")), 
               label.x.npc = "right", label.y.npc = 0,
               formula = y~x-1, parse = TRUE, size = 5) +
  xlab("Stężenie PM10 [µg/m3] (Obserwacja)") +
  ylab("Stężenie PM10 [µg/m3] (Prognoza)") +
  ggtitle("Wykres rozrzutu prognozy od obserwacji dla stężeń PM10\nna stacji w Zdzieszowicach w roku 2018") +
  theme_minimal()

Szara przerywana linia obrazuje model idealny (tzn. gdy prognoza równa jest pomiarowi), a niebieska – linię trendu uzyskaną z wybranych danych (jej równanie znajduje się na wykresie). Otrzymany współczynnik determinacji (R2) wynosi 0.78.

Na wykresie można zaobserować, że prognoza jest zaniżona w stosunku do pomiarów stężenia.

W celu potwierdzenia analizy, stworzono wykres trendu prognoz (linia niebieska) od obserwacji stężeń PM10 (linia szara), który pozwala na ogólną ocenę zgodności tych dwóch serii danych.

ggplot(data = dane_rf_test) +
  geom_path(aes(x = date, y = obs, color = "Obserwacje"), size = 1.5) + 
  geom_path(aes(x = date, y = mod, color = "Model"), size = 1.5, alpha=0.6) +
  scale_x_datetime(limits = ymd_h(c("2018-01-01 00", "2018-01-31 23"))) +
  scale_y_continuous(limits=c(0,250)) +
  scale_color_manual(values=c("Obserwacje"="#CACACA","Model"="#1FA2C3"))+
  xlab("Data") +
  ylab("Stężenie PM10 [µg/m3]") +
  ggtitle("Wykres trendu prognozy od obserwacji dla stężeń PM10\nna stacji w Zdzieszowicach w roku 2018") +
  theme_minimal() +
  theme(legend.position="top",legend.title=element_blank())

Z wykresu można wywnioskować, że najwyższe wartości stężeń są zaniżane w stworzonym modelu. Niewątpliwie został zachowany ogólny trend, niemniej dla niższych stężeń także zdarzają się odchylenia.