Przygotowanie danych

d1 <- dane_1g_mazowieckie
d2 <- dane_24g_mazowieckie

d1 <- d1 %>% filter(!(kod == "MzBialaKmiciMOB" | kod == "MzMinMazKaziMOB" | kod == "MzWarBajkowa" | kod == "MzWarChrosci"))
d2 <- d2 %>% filter(!(kod == "MzMinMazKaziMOB" | kod == "MzWarBajkowa" | kod == "MzWarChrosci"))

d1 <- d1 %>%
  spread(key=sub, value=obs)

d1s <- d1 %>% openair::timeAverage(avg.time="year",type="kod")

d2 <- d2 %>%
  spread(key=sub, value=obs)

dane_mazowieckie_ogolnie <- openair::aqStats(mydata = d1, pollutant = c("NO2","SO2","NOx","O3","PM10","PM2.5"), type = "kod")

dane_mazowieckie_ogolnie_war_plo <- dane_mazowieckie_ogolnie %>% filter(grepl("MzWar", kod, fixed = TRUE) | grepl("MzPlo", kod, fixed = TRUE) | grepl("MzRad", kod, fixed = TRUE))

dane_mazowieckie_ogolnie_reszta <- dane_mazowieckie_ogolnie %>% filter(!(grepl("MzWar", kod, fixed = TRUE) | grepl("MzPlo", kod, fixed = TRUE) | grepl("MzRad", kod, fixed = TRUE)))

dane_mazowieckie_ogolnie_war_plo <- dane_mazowieckie_ogolnie_war_plo %>% mutate(kompletnosc_ponad_90 = if_else(dat.cap >= 90, TRUE, FALSE))

dane_mazowieckie_ogolnie_reszta <- dane_mazowieckie_ogolnie_reszta %>% mutate(kompletnosc_ponad_90 = if_else(dat.cap >= 90, TRUE, FALSE))

d1_NO2 <- d1 %>% mutate(przekroczenie_NO2 = if_else(NO2 > 200,1,0))
d1_NO2 <- d1_NO2 %>% filter(przekroczenie_NO2 == 1)
d1_NO2 <- d1_NO2 %>% openair::timeAverage(avg.time="year",type="kod", statistic = "frequency")

d1_SO2 <- d1 %>% mutate(przekroczenie_SO2 = if_else(SO2 > 350,1,0))
d1_SO2 <- d1_SO2 %>% filter(przekroczenie_SO2 == 1)
d1_SO2 <- d1_SO2 %>% openair::timeAverage(avg.time="year",type="kod", statistic = "frequency")

d2_PM10 <- d2 %>% mutate(przekroczenie_PM10 = if_else(PM10 > 50,1,0))
d2_PM10 <- d2_PM10 %>% filter(przekroczenie_PM10 == 1)
d2_PM10 <- d2_PM10 %>% openair::timeAverage(avg.time="year",type="kod", statistic = "frequency")

kwantyle_O3 <- as.data.frame(d1$O3 %>% quantile(., na.rm=T))
kwantyle_NOx <- as.data.frame(d1$NOx %>% quantile(., na.rm=T))
kwantyle_NO <- as.data.frame(d1$NO %>% quantile(., na.rm=T))
kwantyle_NO2 <- as.data.frame(d1$NO2 %>% quantile(., na.rm=T))
kwantyle_SO2 <- as.data.frame(d1$SO2 %>% quantile(., na.rm=T))
kwantyle_PM10 <- as.data.frame(d1$PM10 %>% quantile(., na.rm=T))
kwantyle_PM2_5 <- as.data.frame(d1$PM2.5 %>% quantile(., na.rm=T))

kwantyle <- cbind(kwantyle_O3,kwantyle_NOx,kwantyle_NO,kwantyle_NO2,kwantyle_SO2,kwantyle_PM10,kwantyle_PM2_5)
colnames(kwantyle) <- c("kwantyle O3","kwantyle NOx","kwantyle_NO","kwantyle NO2","kwantyle SO2","kwantyle PM10","kwantyle PM2.5")

kwantyle <- kwantyle[-c(1,5),]

kwantyle_PM10_24h <- as.data.frame(d2$PM10 %>% quantile(., na.rm=T))
kwantyle_PM2_5_24h <- as.data.frame(d2$PM2.5 %>% quantile(., na.rm=T))

kwantyle_24h <- cbind(kwantyle_PM10_24h,kwantyle_PM2_5_24h)
colnames(kwantyle_24h) <- c("kwantyle PM10","kwantyle PM2.5")

dair_mazowieckie <- dair_mazowieckie %>%
  mutate(rok = lubridate::year(date),
         month = lubridate::month(date)) %>%
  group_by(poll, EoICode, rok) %>%
  mutate(rok = as.character(rok),
         month = as.factor(month))

d1s_dt <- d1s %>% mutate_if(is.numeric,funs(round(.,2)))

Ocena jakości powietrza

Podstawowe statystyki danych 1-godzinnych

datatable(d1s_dt)

Poglądowe wykresy zanieczyszczeń dla konkretnych stacji

d1s %>%
  timePlot(pollutant = c("NO","NO2","NOx","O3","PM10","PM2.5", "SO2"), type = "kod", date.format = "%Y")

Zróżnicowanie stężeń zanieczyszczeń w zależności od godziny, dnia, dnia tygodnia i miesiąca dla wybranych substancji

Dane były zbierane co 1 godzinę.

Kwantyle 25%, 50% i 75% dla danych substancji z pomiarów 1h

Zróżnicowanie stężeń zanieczyszczeń miesięcznie i tygodniowo dla PM10 i PM2.5

Dane były zbierane co 24 godziny.

d2 %>% timeVariation(pollutant = "PM10", group="kod")

d2 %>% timeVariation(pollutant = "PM2.5", group="kod")

Kwantyle 0%, 25%, 50%, 75%, 100% dla danych substancji z pomiarów 1h

Kompletność danych na stacjach dla lat 2015-2019

Dane zostały podzielone na dwie grupy dla większej czytelności wykresów. Grupa pierwsza to stacje w Warszawie, Radomiu oraz Płocku. Grupa druga to pozostałe stacje w województwie mazowieckim.

ggplot(data = dane_mazowieckie_ogolnie_war_plo, aes(x = year, y = dat.cap, fill=pollutant)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge') +
  facet_wrap(~pollutant) +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 90, color = 'red', size = 1.5) +
  labs(
    title = "Kompletność danych na stacjach z grupy 1 dla lat 2015-2019",
    x = "Data",
    y = "Kompletność danych [%]",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_war_plo, aes(x = year, y = dat.cap)) +
  geom_bar(stat='identity', aes(position=kod, fill=kompletnosc_ponad_90), position = 'dodge', color = 'black') +
  facet_wrap(~pollutant) +
  scale_fill_manual(values = c("red","green"), labels = c("Nie","Tak")) +
  geom_hline(yintercept = 90, color = 'red', size = 1.5) +
  labs(
    title = "Wykres kompletności danych dla progu 90% (grupa 1)",
    x = "Data",
    y = "Kompletność danych [%]",
    fill = "Kompletność >90%?"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta, aes(x = year, y = dat.cap, fill=pollutant)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge') +
  facet_wrap(~pollutant) +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 90, color = 'red', size = 1.5) +
  labs(
    title = "Kompletność danych na stacjach z grupy 2 dla lat 2015-2019",
    x = "Data",
    y = "Kompletność danych [%]",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta, aes(x = year, y = dat.cap)) +
  geom_bar(stat='identity', aes(position=kod, fill=kompletnosc_ponad_90), position = 'dodge', color = 'black') +
  facet_wrap(~pollutant) +
  scale_fill_manual(values = c("red","green"), labels = c("Nie","Tak")) +
  geom_hline(yintercept = 90, color = 'red', size = 1.5) +
  labs(
    title = "Wykres kompletności danych dla progu 90% (grupa 2)",
    x = "Data",
    y = "Kompletność danych [%]",
    fill = "Kompletność >90%?"
  )

Kompletność danych na stacjach jest bardzo dobra. Tylko w kilku przypadkach zebrane dane nie osiągnęły progu 90% (MzWarKondrat - O3, PM10 i SO2, MzGranicaKPN - SO2 i MzGutyDurCzer - NO2, NOx), nie biorąc pod uwagę stacji, na których dane zanieczyszczenia nie były zbierane.

Do kolejnych analiz usuwamy stacje w których kompletność danych jest mniejsza niż 90% dla danego zanieczyszczenia:

dane_mazowieckie_ogolnie_war_plo <- dane_mazowieckie_ogolnie_war_plo %>% filter(dat.cap >= 90)
dane_mazowieckie_ogolnie_reszta <- dane_mazowieckie_ogolnie_reszta %>% filter(dat.cap >= 90)

dane_mazowieckie_ogolnie_reszta_NOx <- dane_mazowieckie_ogolnie_reszta %>% filter(pollutant == "NOx")
dane_mazowieckie_ogolnie_reszta_PM10 <- dane_mazowieckie_ogolnie_reszta %>% filter(pollutant == "PM10")
dane_mazowieckie_ogolnie_reszta_PM2_5 <- dane_mazowieckie_ogolnie_reszta %>% filter(pollutant == "PM2.5")
dane_mazowieckie_ogolnie_reszta_SO2 <- dane_mazowieckie_ogolnie_reszta %>% filter(pollutant == "SO2")
dane_mazowieckie_ogolnie_reszta_NO2 <- dane_mazowieckie_ogolnie_reszta %>% filter(pollutant == "NO2")

dane_mazowieckie_ogolnie_war_plo_NOx <- dane_mazowieckie_ogolnie_war_plo %>% filter(pollutant == "NOx")
dane_mazowieckie_ogolnie_war_plo_PM10 <- dane_mazowieckie_ogolnie_war_plo %>% filter(pollutant == "PM10")
dane_mazowieckie_ogolnie_war_plo_PM2_5 <- dane_mazowieckie_ogolnie_war_plo %>% filter(pollutant == "PM2.5")
dane_mazowieckie_ogolnie_war_plo_SO2 <- dane_mazowieckie_ogolnie_war_plo %>% filter(pollutant == "SO2")
dane_mazowieckie_ogolnie_war_plo_NO2 <- dane_mazowieckie_ogolnie_war_plo %>% filter(pollutant == "NO2")

Wykresy średniej rocznej ilości zanieczyszczeń na stacjach dla lat 2015-2019

Czerwona linia pokazuje dopuszczelną wartość stężenia danej substancji w powietrzu z rozporządzenia ministra środowiska w sprawie poziomów niektórych zanieczyszczeń w środowisku (link).

ggplot(data = dane_mazowieckie_ogolnie_war_plo_NOx, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 30, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie NOx w powietrzu w zależności od roku dla stacji z grupy 1",
    x = "Data",
    y = "Średnia roczna NOx",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta_NOx, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 30, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie NOx w powietrzu w zależności od roku dla stacji z grupy 2",
    x = "Data",
    y = "Średnia roczna NOx",
    fill = "Nazwa stacji"
  )

Z rozporządzenia próg dopuszczalnego rocznego stężenia NOx to 30 µg/m3. Roczne stężenie NOx w Warszawie bardzo przekracza normę (stacje MzWarAlNiepo i MzWarMarszal), na innych stacjach norma nie jest lub delikatnie przekroczona (MzWarKondrat, MzRadTochter, MzWarWokalna i MzPiasPulask). Wysokie stężenie może być spowodowane dużym ruchem samochodów w okolicy stacji, które głównie odpowiadają za występowanie tego zanieczyszczenia.

ggplot(data = dane_mazowieckie_ogolnie_war_plo_SO2, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 20, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie SO2 w powietrzu w zależności od roku dla stacji z grupy 1",
    x = "Data",
    y = "Średnia roczna SO2",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta_SO2, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 20, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie SO2 w powietrzu w zależności od roku dla stacji z grupy 2",
    x = "Data",
    y = "Średnia roczna SO2",
    fill = "Nazwa stacji"
  )

Z rozporządzenia próg dopuszczalnego rocznego stężenia SOx to 20 µg/m3. Roczne stężenie SO2 nie było przekroczone na żadnej stacji.

ggplot(data = dane_mazowieckie_ogolnie_war_plo_NO2, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 40, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie NO2 w powietrzu w zależności od roku dla stacji z grupy 1",
    x = "Data",
    y = "Średnia roczna NO2",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta_NO2, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 40, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie NO2 w powietrzu w zależności od roku dla stacji z grupy 2",
    x = "Data",
    y = "Średnia roczna NO2",
    fill = "Nazwa stacji"
  )

Z rozporządzenia próg dopuszczalnego rocznego stężenia NO2 to 40 µg/m3. Stężenie NO2 wygląda podobnie do stężenia NOx, ale została przekroczona norma tylko na dwóch stacjach - MzWarAlNiepo i MzWarMarszal.

ggplot(data = dane_mazowieckie_ogolnie_war_plo_PM10, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 40, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie PM10 w powietrzu w zależności od roku dla stacji z grupy 1",
    x = "Data",
    y = "Średnia roczna PM10",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta_PM10, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 40, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie PM10 w powietrzu w zależności od roku dla stacji z grupy 2",
    x = "Data",
    y = "Średnia roczna PM10",
    fill = "Nazwa stacji"
  )

Z rozporządzenia próg dopuszczalnego rocznego stężenia PM10 to 40 µg/m3. Stężenie PM10 było przekroczone nieznacznie tylko w jednej stacji w Warszawie (MzWarAlNiepo). Duża część stacji miała poziom zanieczyszczenia przy granicy dopuszczalnego stężenia.

ggplot(data = dane_mazowieckie_ogolnie_war_plo_PM2_5, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 20, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie PM2.5 w powietrzu w zależności od roku dla stacji z grupy 1",
    x = "Data",
    y = "Średnia roczna PM2_5",
    fill = "Nazwa stacji"
  )

ggplot(data = dane_mazowieckie_ogolnie_reszta_PM2_5, aes(x = year, y = mean)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 20, color = 'red', size = 1.5) +
  labs(
    title = "Stężenie PM2.5 w powietrzu w zależności od roku dla stacji z grupy 2",
    x = "Data",
    y = "Średnia roczna PM2_5",
    fill = "Nazwa stacji"
  )

Z rozporządzenia próg dopuszczalnego rocznego stężenia PM2.5 to 20 µg/m3. Roczne stężenie PM2.5 przekraczało normę na większości stacji we wszystkich latach oprócz roku 2019, w którym jest widoczny spadek stężeń w porównaniu do poprzednich lat.

Na wszystkich wykresach jest widoczny spadek stężeń zanieczyszczeń na przestrzeni kolejnych lat.

Wykresy liczby przekroczeń danych zanieczyszczeń godzinowo lub dziennie w latach 2015-2019 na danych stacjach

Brak stacji na wykresach wskazuje na brak przekroczeń godzinowych lub dziennych.

ggplot(data = d1_NO2, aes(x = date, y = przekroczenie_NO2)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 18, color = 'red', size = 1.5) +
  labs(
    title = "Dzienne stężenie NO2 w powietrzu w zależności od roku",
    x = "Data",
    y = "Dziennie przekroczony dopuszczalny poziom NO2 [n]",
    fill = "Nazwa stacji"
  )

ggplot(data = d1_SO2, aes(x = date, y = przekroczenie_SO2)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff")) +
  geom_hline(yintercept = 24, color = 'red', size = 1.5) +
  labs(
    title = "Dzienne stężenie SO2 w powietrzu w zależności od roku",
    x = "Data",
    y = "Dziennie przekroczony dopuszczalny poziom SO2 [n]",
    fill = "Nazwa stacji"
  )

Dzienne stężenie NO2 i SO2 było przekroczone kilka razy, ale liczba przekroczeń nie przewyższa dopuszczalnej normy. NO2 byo przekroczone tylko w Warszawie na stacji MzWarAlNiepo, a SO2 w Płocku na dwóch stacjach.

ggplot(data = d2_PM10, aes(x = date, y = przekroczenie_PM10)) +
  geom_bar(stat='identity', aes(fill=kod), position = 'dodge', color = 'black') +
  scale_fill_manual(values = c("#000000", "#ff33ff", "#0066ff", "#33ffff", "#666600", "#33ff33", "#ffff33", "#8C1717", "#6600ff", "#990033", "#99ccff", "#cccc66", "#ffcc33", "#008080")) +
  geom_hline(yintercept = 35, color = 'red', size = 1.5) +
  labs(
    title = "Miesięczne stężenie SO2 w powietrzu w zależności od roku",
    x = "Data",
    y = "Miesięcznie przekroczony dopuszczalny poziom PM10 [n]",
    fill = "Nazwa stacji"
  )

W 2015 i 2018 stężenie dopuszczalne miesięczne było przekroczone na większości stacji, nastomiast w 2016 i 2017 stężenie było trochę mniejsze. W 2019 widać wyraźny spadek stężenia zanieczyszczeń w porównaniu do poprzednich lat.

Patrząc na wszystkie wykresy można zauważyć, że najbardziej zanieczyszczone są okolice stacji MzWarAlNiepo w Warszawie.

Porównanie danych automatycznych i referencyjnych

W województwie mazowieckim byly wykonywane jednoczenie pomiary metod automatyczna jak i referencyjna na stacjach:

dair_mazowieckie$EoICode %>% unique()
##  [1] "PL0673A" "PL0129A" "PL0753A" "PL0539A" "PL0134A" "PL0538A" "PL0140A"
##  [8] "PL0717A" "PL0739A" "PL0143A" "PL0141A"

Do oceny czy metoda automatyczna jest wystarczajca w stosunku do metody referencyjnej porównamy ze soba dwa zanieczyszczenia powietrza:

dair_mazowieckie$poll %>% unique()
## [1] "PM10"  "PM2.5"

Średnie roczne stężenia PM10 oraz PM2.5 w obu metodach oraz ich minimum i maksimum dla poszczególnych stacji oraz lat prezentują się nastepująco:

datatable(dair_mazowieckie %>%
  summarise(sre_a = round(mean(aut, na.rm = T),2),
            sre_r = round(mean(ref, na.rm = T),2),
            min_a = round(min(aut, na.rm = T),2),
            min_r = round(min(ref, na.rm = T),2),
            max_a = round(max(aut, na.rm = T),2),
            max_r = round(max(ref, na.rm = T),2)))

Tabelka przedstawiająca wspólne statystyki oceny modelu:

dair_mazowieckie_dt <- dair_mazowieckie %>%
  modStats(mod = "aut",
           obs = "ref",
           type = c("rok", "poll", "EoICode"))

dair_mazowieckie_dt <- dair_mazowieckie_dt %>% mutate_if(is.numeric,funs(round(.,2)))
datatable(dair_mazowieckie_dt)

Jak możemy zauważyć kompletność naszych danych jest bardzo duża, w najlepszym przypadku wynosi ona aż 98,9% czyli 361 dni w cigu roku dla stacji w PL0539A w 2018r. oraz dla stacji PL0717A w 2019r. Najmniejsz kompletność danych możemy zauważyć dla stacji PL0673A w 2017r., która wynosi 86%. Na stacjach w których kompletnosc danych pomiarowych nie przekracza progu 90% nie powinna być prowadzona ocena poziomów substancji w powietrzu. Sa to stacje o kodach: PL0538A w roku 2016 dla PM10, PL0673A w roku 2017 i 2019 dla PM10. Tabelka również prezentuje srednie odchylenie naszych danych , które oscyluja w okolicach mniejszych niz 0,5, czyli jest to naprawde bardzo dobry wynik. Wartoć wspóczynnika korelacji Pearsona jest w wiekszosci przypadkow powyzej 0,95 , czyli zaleznosc liniowa miedzy dwiema metodami jest silna korelacja.

Usuwanie niekompletnych danych

dair_mazowieckie_1 <- dair_mazowieckie %>% 
  filter(!(EoICode == "PL0538A" & rok == "2016" & poll == "PM10"))
dair_mazowieckie_1 <- dair_mazowieckie_1 %>% 
  filter(!(EoICode == "PL0673A" & rok == "2017" & poll == "PM10"))
dair_mazowieckie_1 <- dair_mazowieckie_1 %>% 
  filter(!(EoICode == "PL0673A" & rok == "2019" & poll == "PM10"))

dair_mazowieckie_1_dobowe <- dair_mazowieckie_1 %>%
  mutate(doba = lubridate::day(date)) %>%
  group_by(poll, EoICode, rok, month) %>%
  mutate(doba = as.factor(doba))

Tabelka przedstawiajaca wspólne statystyki oceny modelu w poszczególnych porach roku

dair_mazowieckie_1_dt <- dair_mazowieckie_1 %>%
  modStats(., mod = "aut", 
           obs = "ref", type = c("year", "EoICode",  "poll", "season"))

dair_mazowieckie_1_dt <- dair_mazowieckie_1_dt %>% mutate_if(is.numeric,funs(round(.,2)))
datatable(dair_mazowieckie_1_dt)

Wykresy rozrzutu dla poszczególnych lat

dair_mazowieckie_1_dobowe %>%
  na.omit() %>%
  filter(lubridate::year(date) == 2015,
         poll == "PM10") %>%
  ggplot(aes(x = ref, y = aut)) +
  ggtitle("2015 r.")+
  geom_point(pch = 20, size = 2, aes(fill = doba, color = doba)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_dark() +
  facet_wrap(~EoICode) +
  geom_abline(intercept = c(0,0,0),
              slope = c(0.85,1,1.15))

dair_mazowieckie_1_dobowe %>%
  na.omit() %>%
  filter(lubridate::year(date) == 2016,
         poll == "PM10") %>%
  ggplot(aes(x = ref, y = aut)) +
  ggtitle("2016 r.")+
  geom_point(pch = 20, size = 2, aes(fill = doba, color = doba)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_dark() +
  facet_wrap(~EoICode) +
  geom_abline(intercept = c(0,0,0),
              slope = c(0.85,1,1.15))

dair_mazowieckie_1_dobowe %>%
  na.omit() %>%
  filter(lubridate::year(date) == 2017,
         poll == "PM10") %>%
  ggplot(aes(x = ref, y = aut)) +
  ggtitle("2017 r.")+
  geom_point(pch = 20, size = 2, aes(fill = doba, color = doba)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_dark() +
  facet_wrap(~EoICode) +
  geom_abline(intercept = c(0,0,0),
              slope = c(0.85,1,1.15))

dair_mazowieckie_1_dobowe %>%
  na.omit() %>%
  filter(lubridate::year(date) == 2018,
         poll == "PM10") %>%
  ggplot(aes(x = ref, y = aut)) +
  ggtitle("2018 r.")+
  geom_point(pch = 20, size = 2, aes(fill = doba, color = doba)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_dark() +
  facet_wrap(~EoICode) +
  geom_abline(intercept = c(0,0,0),
              slope = c(0.85,1,1.15))

dair_mazowieckie_1_dobowe %>%
  na.omit() %>%
  filter(lubridate::year(date) == 2019,
         poll == "PM10") %>%
  ggplot(aes(x = ref, y = aut)) +
  ggtitle("2019 r.")+
  geom_point(pch = 20, size = 2, aes(fill = doba, color = doba)) +
  scale_x_continuous(expand = c(0,0)) +
  scale_y_continuous(expand = c(0,0)) +
  theme_dark() +
  facet_wrap(~EoICode) +
  geom_abline(intercept = c(0,0,0),
              slope = c(0.85,1,1.15))

Patrząc na wykresy kolejnych lat jako pierwsze możemy zauważyc, że liczba stacji w których wykonywane sa pomiary obiema metodami wzrasta, co jest dobrą wiadomością, ponieważ dzięki temu można przeprowadzić więcej analiz. Na większości wykresów nie widać znaczących różnic, wartości są sobie równe. Natomiast na niektórych stacjach możemy zaobserwować pojedyncze wartości, które znacząco odbiegają od linii. Są to między innymi stacje takie jak PL0538A w latach 2015 oraz 2017r. , stacje, które wykonywały pomiary w 2017 r. poza stacją o kodzie PL0140A. Najmniejsze zależności między dwiema metodami możemy zauważyć na wykresie stacji PL0753A w 2019 r.

Diagram Taylora

dair_mazowieckie %>%
  filter(lubridate::year(date) == 2019, # filtrowanie
         poll == "PM10",
         EoICode == "PL0753A") %>%
  TaylorDiagram(mod = "aut",
                obs = "ref",
                rms.col  =  "darkgoldenrod",
                normalize  =  FALSE,
                cols = c("pink", "green", "blue"))

Diagram prezentuje przykładowe statystyki dla pomiarów w 2019 roku na stacji PL0753A. Współczynnik korelaci wynosi 0.76, natomiast wartość średniego błędu kwadratowego (RMSE) wynosi 9.34, a odchylenie standardowe jest równe 13.5.

Wnioski

Do badania zawartości pyłu zawieszonego PM10 oraz PM2.5 w powietrzu stosuje się dwie metody, metodę grawimetryczną(refernecyjną) oraz metodę automatyczną. Metoda refernecyjna jest powszechnie uważana za metodę najbardziej skuteczną oraz precyzyjną pod względem wyników pomiarów. W naszej analizie wykonaliśmy parę tabelek oraz wykresów, które wykazują, że metoda automatyczna wcale nie odbiega znacząco od precyzji wykonywania pomiarów jak w przypadku metody manualnej. Możemy zatem stwierdzić, że metoda automatyczna jest równoważna z metodą referencyjną.