1 R Pakiet

library(tidyverse)
library(plotly)
library(gdata)
library(openair)
library(hexbin)
library(kableExtra)
library(knitr)
library(reshape2)
library(boot)
library(parallel)

2 Wczytywanie, czyszczenie, grupowanie

2.1 Wczytywanie

read.table(file = "T:/2_PROJEKTY_WYNIKI/OSPM_KRAS_PM/Wyniki_ospm_v1.csv", 
           sep = ";", dec = ".", header = T) -> ospm
save(ospm, file = "T:/2_PROJEKTY_WYNIKI/OSPM_KRAS_PM/ospm.RData")

2.2 Czyszczenie i konwersja

Teraz czyscimy i konwertujemy dane. Tworzymy nowe zmienne pomocnicze Ratio (obs/mod) i Ratio (obs/tlo)

ospm[ospm$obs > ospm$back,] -> ospm # usuń niepoprawne tlo
ospm <- na.omit(ospm) # usun dane brakujace

ospm$Date <- as.character(ospm$Date)
ospm$Date <- as.POSIXct(strftime(ospm$Date, format = "%Y-%m-%d", tz="GMT"))

ospm$u_roof <- as.character(ospm$u_roof)
ospm$typ_sL <- as.character(ospm$typ_sL)
ospm$dust <- as.character(ospm$dust)

ospm$mod/ospm$back -> ospm$mq
ospm$obs/ospm$back -> ospm$oq
ospm$ratio <- ospm$obs/ospm$mod
ospm$ratio1 <- ospm$obs/ospm$back

ospm %>% filter(month %in% c(5,6,7,8)) %>% filter(u_roof == "u_roof = 1.0") -> PM
ospm %>% filter(month %in% c(5,6,7,8)) %>% filter(typ_sL %in% c("Qs")) -> PMr

levels(as.factor(PM$typ_sL))

## [1] "Qs"      "Qs15"    "Qs15+Qc" "Qs30"    "Qs30+Qc"

PM$typ_sL  <- factor(PM$typ_sL,  levels = c("Qs","Qs15","Qs30","Qs15+Qc","Qs30+Qc"))
save(ospm,PM,PMr, file = "ospm.RData")

2.3 Grupowanie

2.3.1 Kierunek wiatru co 45 stopi, 8 sektorów

deg <- 45
dir.breaks <- seq(0-(deg/2), 360+(deg/2), deg)
PM$gwd <- cut(PM$wd,breaks = dir.breaks, ordered_result = TRUE)
dir.labels <- as.character(c(seq(0, 360-deg, by = deg), 0))
levels(PM$gwd) <- dir.labels
remove(deg, dir.breaks, dir.labels)
summary(PM %>% filter(dust == "PM10",  typ_sL == "Qs") %>% select(gwd))

##       gwd     
##  270    :812  
##  225    :407  
##  90     :318  
##  180    :150  
##  135    :146  
##  45     :144  
##  (Other):101

levels(PM$gwd)

## [1] "0"   "45"  "90"  "135" "180" "225" "270" "315"

2.3.2 Prędkość wiatru

PM$gws <- cut(PM$ws, breaks = c(0,0.5,1.5,2,2.5,7))
summary(PM %>% filter(dust == "PM10",  typ_sL == "Qs") %>% select(gws))

##         gws     
##  (0,0.5]  :111  
##  (0.5,1.5]:907  
##  (1.5,2]  :356  
##  (2,2.5]  :236  
##  (2.5,7]  :468

2.3.3 Czas

krok co 3 oodziny, time end

PM$godz <- cut(PM$hour, breaks = c(0,3,6,9,12,15,18,21,24), 
               labels = c("1-3","4-6","7-9","10-12","13-15","16-18","19-21","22-24"))
summary(PM %>% filter(dust == "PM10",  typ_sL == "Qs") %>% select(godz))

##       godz    
##  7-9    :280  
##  4-6    :272  
##  19-21  :264  
##  10-12  :258  
##  1-3    :254  
##  22-24  :254  
##  (Other):496

2.3.4 Opad atmosferyczny

Jesli wieksze 0.24 to nie liczono wtornej emisji pylu zgodnie z metodyka

PM$gopad <- cut(PM$opad, breaks = c(-1,0.24,max(PM$opad)), labels = c("Nie","Tak"))
summary(PM %>% filter(dust == "PM10",  typ_sL == "Qs") %>% select(gopad))

##  gopad     
##  Nie:1614  
##  Tak: 464

2.3.5 Wysokość mieszania

Na wszelki wypadek, ale ogolnie nie tedy droga

quantile(PM$H_mix, probs = seq(0,1,0.2)) -> x ; x[1]<-0 ;x

##      0%     20%     40%     60%     80%    100% 
##    0.00   46.68  165.08  737.03 1466.21 1997.98

PM$gh_mix <-cut(PM$H_mix, breaks = x) ; remove(x)
summary(PM %>% filter(dust == "PM10",  typ_sL == "Qs0") %>% select(gh_mix))

##               gh_mix 
##  (0,46.7]        :0  
##  (46.7,165]      :0  
##  (165,737]       :0  
##  (737,1.47e+03]  :0  
##  (1.47e+03,2e+03]:0

2.3.6 Wilgotność powietrza

Na wszelki wypadek, ale ogolnie nie tedy droga

quantile(PM$R_hum, probs = seq(0,1,0.2)) -> x ; x[1]<-0 ;x

##     0%    20%    40%    60%    80%   100% 
## 0.0000 0.4617 0.6070 0.7261 0.8191 0.9187

PM$ghum <-cut(PM$R_hum, breaks = x) ; remove(x)
summary(PM %>% filter(dust == "PM10",  typ_sL == "Qs") %>% select(ghum))

##             ghum    
##  (0,0.462]    :478  
##  (0.462,0.607]:419  
##  (0.607,0.726]:408  
##  (0.726,0.819]:384  
##  (0.819,0.919]:389

3 Analiza podstawowa

3.1 Parametry statystyczne (funkcje)

NAD <- function(obs = obs, mod = mod) {
  NAD = round((mean(Mod(obs-mod))) / (mean(obs) + mean(mod)),3)}

FB <- function(obs = obs, mod = mod) {
  FBFP = round((0.5*(sum(Mod(obs - mod)  +  (mod - obs)))) / 
                 (0.5 * sum(obs+mod)),3)
  FBFN = round((0.5*(sum(Mod(obs - mod) + (obs - mod)))) / 
                 (0.5 * sum(obs+mod)),3)
  FB = round(FBFN-FBFP,3)}

FBFN <- function(obs = obs, mod = mod) {
  FBFN = round((0.5*(sum(Mod(obs - mod) + (obs - mod)))) / 
                 (0.5 * sum(obs+mod)),3)}

FBFP <- function(obs = obs, mod = mod) {
  FBFP = round((0.5*(sum(Mod(obs - mod)  +  (mod - obs)))) / 
                 (0.5 * sum(obs+mod)),3)}

NMSE <- function(obs = obs, mod = mod) {
  NMSE = round(mean((obs - mod)^2) / (mean(obs) * mean(mod)),3)}

VG <- function(obs = obs, mod = mod) {
  VG   = round(exp(mean( (log(obs)-log(mod))^2 ) ),3)}

MG <- function(obs = obs, mod = mod) {
  MG   = round(exp(mean(log(obs)) - mean(log(mod))),3)}

R <- function(obs = obs, mod = mod) {
  R    = round(cor(x = mod, y = obs, method = c("pearson")),3)}

IA <- function(mod=mod, obs=obs) {
  IA = sum((mod-obs)^2) / 
    sum((abs(mod - mean(obs)) + 
           abs(obs - mean(obs)))^2)
  IA = round(1-IA,3)
}

3.2 Funkcja liczenia parametrów staystytcznych.

stat.b <- function(mdata = mdata, obs = "obs", mod = "mod",
                 typ  = c("dust","typ_sL")) {
  group_by_(mdata, .dots = typ) %>%
    summarise(lobs  = n(),
               sr.o = round(mean(obs),2), 
               sr.m = round(mean(mod),2),
               FBFP = round(FBFP(obs= obs, mod = mod),2), 
               FBFN = round(FBFN(obs = obs, mod = mod),2),
               FB   = round(FB(obs  = obs, mod = mod),2), 
               NMSE = round(NMSE(obs = obs, mod = mod),2),
               NAD  = round(NAD(obs = obs, mod = mod),2), 
               MG   = round(MG(obs   = obs, mod = mod),2), 
               VG   = round(VG(obs  = obs, mod = mod),2),
               IA   = round(IA(obs  = obs, mod = mod),2)) -> P
  
  modStats(mdata, mod = mod, obs = obs, type = typ) -> P1
  
  P  <- arrange_(P, .dots = typ)  
  P1 <- arrange_(P1, .dots = typ)
  X <- merge(P,P1, by = typ)

  round(X$FAC2,2) -> X$FAC2
  round(X$IOA ,2) -> X$IOA
  round(X$COE ,2) -> X$COE
  round(X$r   ,2) -> X$r
  X
}

3.3 Parametry statystyczne

Najpierw liczymy parametry statystyczne. zgodnie z utworzonymi funkcjami.

stat.b(mdata = PM, obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL")) ->  b.all

stat.b(mdata = PMr, obs = "obs", mod = "mod",
     typ  = c( "dust","u_roof")) ->  b.f

Drukujemy tabele z wynikami.

kable(b.all[,c(1:2,14,6:13,15,21:23)], caption = "Tabela 3.1.Parametry staystyczne dla u/u = 1.0")

Tabela 3.1.Parametry staystyczne dla u/u = 1.0
dust	typ_sL	n	FBFP	FBFN	FB	NMSE	NAD	MG	VG	IA	FAC2	r	COE	IOA
PM10	Qs	2078	0.07	0.20	0.13	0.17	0.14	1.10	1.15	0.78	0.94	0.70	0.33	0.66
PM10	Qs15	2078	0.10	0.18	0.08	0.15	0.14	1.04	1.15	0.79	0.94	0.68	0.31	0.66
PM10	Qs15+Qc	2078	0.11	0.14	0.03	0.13	0.13	0.99	1.12	0.83	0.96	0.72	0.36	0.68
PM10	Qs30	2078	0.12	0.15	0.03	0.15	0.14	0.99	1.15	0.79	0.94	0.66	0.29	0.64
PM10	Qs30+Qc	2078	0.13	0.13	-0.01	0.13	0.13	0.96	1.13	0.82	0.95	0.71	0.32	0.66
PM2.5	Qs	2136	0.05	0.25	0.20	0.18	0.15	1.18	1.17	0.80	0.92	0.75	0.33	0.66
PM2.5	Qs15	2136	0.07	0.22	0.16	0.16	0.15	1.13	1.17	0.81	0.92	0.73	0.33	0.67
PM2.5	Qs15+Qc	2136	0.07	0.19	0.12	0.12	0.13	1.09	1.12	0.86	0.96	0.80	0.40	0.70
PM2.5	Qs30	2136	0.07	0.22	0.14	0.16	0.15	1.12	1.17	0.81	0.92	0.72	0.33	0.66
PM2.5	Qs30+Qc	2136	0.07	0.18	0.11	0.12	0.13	1.08	1.12	0.85	0.96	0.79	0.40	0.70

kable(b.f[,c(1:2,14,6:13,15,21:23)], caption = "Tabela 3.3.Parametry staystyczne porównania u/u = 1.0 dla Qs")

Tabela 3.3.Parametry staystyczne porównania u/u = 1.0 dla Qs
dust	u_roof	n	FBFP	FBFN	FB	NMSE	NAD	MG	VG	IA	FAC2	r	COE	IOA
PM10	u_roof = 0.6	2078	0.10	0.17	0.07	0.15	0.14	1.02	1.15	0.78	0.93	0.68	0.31	0.66
PM10	u_roof = 1.0	2078	0.07	0.20	0.13	0.17	0.14	1.10	1.15	0.78	0.94	0.70	0.33	0.66
PM2.5	u_roof = 0.6	2136	0.08	0.21	0.14	0.16	0.14	1.10	1.16	0.80	0.93	0.72	0.33	0.67
PM2.5	u_roof = 1.0	2136	0.05	0.25	0.20	0.18	0.15	1.18	1.17	0.80	0.92	0.75	0.33	0.66

3.4 Scatter plot dla Qs przy różnym u_roof

Współczynnik 0.6 jest domyślnym ustawienim OSPM, gdy stosujemy dane z stacji mierzacej ws i wd na wysokosci okolo 10 m. Korzystalismy z stacji AGH, która mierzy predkosc wiatru na wysokosci budynku stąd przyjęto, że u_mast/u_roof = 1.0

scatterPlot(PMr%>% filter(dust == "PM10"), x = "obs", y = "mod", type = c("u_roof"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(2,1), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 1 ScatterPlot PM10 dla Qs

scatterPlot(PMr%>% filter(dust == "PM2.5"), x = "obs", y = "mod", type = c("u_roof"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(2,1), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 2 ScatterPlot PM2.5 dla Qs

3.5 Scatter plot PM10 dla u_roof = 1.0 przy róznym typ_sL

scatterPlot(PM%>% filter(dust == "PM10"), x = "obs", y = "mod", type = c("typ_sL"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,2), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 1 ScatterPlot PM10 all data

scatterPlot(PM%>% filter(dust == "PM10") %>% filter(hour %in% c(1:6)), 
            x = "obs", y = "mod", type = c("typ_sL"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,2), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 2 ScatterPlot PM10 dla hour(1:6)

scatterPlot(PM%>% filter(dust == "PM10") %>% filter(hour %in% c(1:6)), 
            x = "obs", y = "mod", type = c("typ_sL"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,2), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 3 ScatterPlot PM10 dla hour(7:24)

scatterPlot(PM%>% filter(dust == "PM10") %>% filter(hour %in% c(1:24)),
            x = "obs", y = "mod", type = c("typ_sL", "dni"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,6), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys.32 ScatterPlot PM10 dla dla hour(1:24)

3.6 Scatter plot PM2.5 dla u_roof = 1.0 przy róznym typ_sL

scatterPlot(PM %>% filter(dust == "PM2.5"), x = "obs", y = "mod", type = c("typ_sL"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,2), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 1 ScatterPlot PM2.5 dla Qs

scatterPlot(PM%>% filter(dust == "PM2.5") %>% filter(hour %in% c(1:6)), 
            x = "obs", y = "mod", type = c("typ_sL"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,2), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 2 ScatterPlot PM2.5 dla hour(1:6)

scatterPlot(PM%>% filter(dust == "PM2.5") %>% filter(hour %in% c(7:24)),
            x = "obs", y = "mod", type = c("typ_sL"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,2), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 3 ScatterPlot PM2.5 dla dla hour(7:24)

scatterPlot(PM%>% filter(dust == "PM2.5") %>% filter(hour %in% c(7:24)),
            x = "obs", y = "mod", type = c("typ_sL", "dni"), 
            method = "hexbin", cols = "jet", mod.line = T, layout = c(3,6), 
            xlab = "Obesrvations [ug/m3]", ylab = "Modelled [ug/m3]")

Rys. 4 ScatterPlot PM2.5 dla dla hour(1:24)

3.7 Ratio > 2 (czestość)

PM %>% filter(PM$ratio > 2.0) -> two_fac
two_fac %>% group_by(dust, typ_sL, hour) %>% summarise(n = n()) -> two_fac

ggplot(data = two_fac, aes(x = hour, y = n, fill = typ_sL)) + 
  labs(x = "Hours in a day", y = "Frequency") +
  geom_col() +facet_grid(dust~typ_sL, scales = "free_y") + 
  theme_bw() + scale_x_continuous(breaks = seq(1,24,2)) + 
  theme(legend.position = "top")

PM %>% filter(PM$ratio < 0.5) -> two_fac
two_fac %>% group_by(dust, typ_sL, hour) %>% summarise(n = n()) -> two_fac

ggplot(data = two_fac, aes(x = hour, y = n, fill = typ_sL)) + 
  labs(x = "Hours in a day", y = "Frequency") +
  geom_col() +facet_grid(dust~typ_sL, scales = "free_y") + 
  theme_bw() + scale_x_continuous(breaks = seq(1,24,2)) + 
  theme(legend.position = "top")

3.8 Reszty

PM$mod - PM$obs -> PM$dif
ggplot(PM, aes(x = dif, fill = typ_sL)) + geom_histogram(bins = 10) +facet_grid(dust~typ_sL) + xlim(-60,60)

4 Warunki meteorologiczne

4.1 Kierunke wiatru (histogramy)

rbind(data.frame(PM %>% filter(dust == "PM10", typ_sL == "Qs"), czas = "all hours"),
data.frame(PM %>% filter(dust == "PM10", typ_sL == "Qs") %>% filter(hour  %in% c(1:6)), czas = "only 1-6"),
data.frame(PM %>% filter(dust == "PM10", typ_sL == "Qs") %>% filter(hour %in% c(7:24)), czas = "only 7-24")) -> px


ggplot(px, aes(x = wd)) + geom_histogram(binwidth = 10) + labs(x = "wind direction") -> p
p + facet_grid(.~czas) -> p
ggplotly(p)

Rys. 1 Histogram kierunku wiatru

4.2 Prędkość wiatru (histogramy)

ggplot(px, aes(x = ws)) + geom_histogram(binwidth = 0.2) + facet_grid(.~czas) + labs(x = "wind speed")->p
ggplotly(p)

Rys. 1 Histogram prędkości wiatru

4.3 Róża wiatrów

pollutionRose(px, pollutant = "ws", breaks = c(0,0.5,1,1.5,2,4,6.87), 
              type = "czas", layout = c(3,1), angle =  360/16)

Rys. 1 Wind Rose

5 Analiza danychw grupach

Liczymy parametry staystyczne wzgledem różnych grup.

PM -> PMstare
PM %>% filter(typ_sL %in% c("Qs","Qs15+Qc", "Qs30+Qc")) -> PM

names(PM) -> x
x[1] <- "data"
x[5] <- "hh"
names(PM) <- x

PM$hh <- as.factor(PM$hh)
stat.b(mdata = PM, obs = "obs", mod = "mod",
     typ  = c("dust","typ_sL","hh")) -> b.hour

stat.b(mdata = PM, obs = "obs", mod = "mod",
     typ  = c("dust","typ_sL","dni","hh")) -> b.hourz

stat.b(mdata = PM %>% select(-data), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","godz")) -> b.godz

stat.b(mdata = PM %>% select(-data), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gws")) -> b.gws

stat.b(mdata = PM %>% filter(hh %in% c(1:6)), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gws")) -> b.h1.gws

stat.b(mdata = PM %>% filter(hh %in% c(7:24)), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gws")) -> b.h7.gws

stat.b(mdata = PM, obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gwd")) -> b.gwd

stat.b(mdata = PM %>% filter(hh %in% c(1:6)), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gwd")) -> b.h1.gwd

stat.b(mdata = PM %>% filter(hh %in% c(7:24)), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gwd")) -> b.h7.gwd

stat.b(mdata = PM, obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","ghum")) -> b.hum

stat.b(mdata = PM %>% filter(hh %in% c(1:6)), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","ghum")) -> b.h1.hum

stat.b(mdata = PM %>% filter(hh %in% c(7:24)), obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","ghum")) -> b.h7.hum

stat.b(mdata = PM, obs = "obs", mod = "mod",
     typ  = c( "dust","typ_sL","gh_mix")) -> b.mix

5.1 PZREDZIAŁY CI BOOSTRAP 95%

fun.b <- function(x, i,...) {
    mean(x[i])
  }

P1 <- data.frame()
levels(as.factor(PM$dust)) -> s.dust
levels(as.factor(PM$typ_sL)) -> s.typ ; s.typ[c(1,4,5)] -> s.typ
levels(as.factor(PM$hh)) -> s.hh
  
  for (i in 1:length(s.dust)) {
    for (j in 1:length(s.typ)) {
      for(k in 1:length(s.hh)) {
        
    PM %>% filter(dust == s.dust[i]) %>% filter(typ_sL == s.typ[j])  %>% filter(hh == s.hh[k]) -> P

    b.mod <- boot(P$mod, statistic = fun.b, R = 200)
    b.obs <- boot(P$obs, statistic = fun.b, R = 200)
    b <- b.obs$t - b.mod$t

rbind(data.frame(dust = s.dust[i], typ_sL = s.typ[j], hour = s.hh[k], type = "mod",
            Mean   =  b.mod$t0, 
            Lower =  boot.ci(b.mod,conf = 0.95, type = "norm")$normal[2],
            Upper =  boot.ci(b.mod,conf = 0.95, type = "norm")$normal[3]),
     data.frame(dust = s.dust[i], typ_sL = s.typ[j], hour = s.hh[k], type = "obs",
            Mean   =  b.mod$t0, 
            Lower =  boot.ci(b.obs,conf = 0.95, type = "norm")$normal[2],
            Upper =  boot.ci(b.obs,conf = 0.95, type = "norm")$normal[3]),
     data.frame(dust = s.dust[i], typ_sL = s.typ[j], hour = s.hh[k], type = "obs-mod",
            Mean  = b.obs$t0-b.mod$t0, 
            Lower = (b.obs$t0-b.mod$t0) - qnorm(p = 0.95)*sd(b),
            Upper = (b.obs$t0-b.mod$t0) + qnorm(p = 0.95)*sd(b))) -> Z
  P1 <- rbind(P1, Z)
}}} # END for()

remove(Z, P, s.dust, s.typ, s.hh)

P1$hour <- as.numeric(P1$hour)
ggplot(P1, aes(y = Mean, x = hour)) + facet_grid(typ_sL~dust, scales = "free_y") + 
  geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = type), alpha = 0.3) +
  geom_line(aes(y = Mean, x = hour,col = type)) + 
  theme_bw() + 
  scale_x_continuous(breaks = seq(0,24,2)) + theme(legend.position = "top")

CI95 BOOSTRAP

5.2 PRZEDZIAŁY CI NORMALNE 95%

PM %>% group_by(dust, typ_sL, hh) %>% summarise(n = n(),
                                                Mean = mean(mod),
                                                sigma = sd(mod)) -> Z1
PM %>% group_by(dust, typ_sL, hh) %>% summarise(n = n(),
                                                Mean = mean(obs),
                                                sigma = sd(obs)) -> Z2
PM %>% group_by(dust, typ_sL, hh) %>% summarise(n = n(),
                                                Mean = mean(obs-mod),
                                                sigma = sd(obs-mod)) -> Z3

rbind(data.frame(Z1, type = "mod"),
      data.frame(Z2, type = "obs"),
      data.frame(Z3, type = "obs-mod")) -> Z

Z$Lower <- round(Z$Mean+c(-1)*Z$sigma/sqrt(32)*qnorm(.975),2)
Z$Upper <- round(Z$Mean+c(1)*Z$sigma/sqrt(32)*qnorm(.975),2)

Z$hh <- as.numeric(Z$hh)
ggplot(Z, aes(y = Mean, x = hh)) + facet_grid(typ_sL~dust, scales = "free_y") + 
  geom_ribbon(aes(ymin = Lower, ymax = Upper, fill = type), alpha = 0.3) +
  geom_line(aes(y = Mean, x = hh,col = type)) + 
  theme_bw() + 
  scale_x_continuous(breaks = seq(0,24,2)) + theme(legend.position = "top")

CI95 NORMAL

5.3 GODZINY 1:24

b.hour$hh <- as.numeric(b.hour$hh)
arrange(b.hour, dust, typ_sL, hh) -> b.hour

ggplot(b.hour, aes(y = FB, x = hh, col = typ_sL)) +
  geom_point(size = 2) + geom_line(aes(group = typ_sL)) +
  facet_grid(dust~.) + theme_bw() +theme(legend.position = "top") + 
  geom_hline(yintercept = c(-0.3,0,0.3), linetype = 2) -> p
ggplotly(p)

Rys. 1 Zmienność FB względem godzin

ggplot(b.hour, aes(y = NMSE, x = FB, shape = typ_sL,col = typ_sL, frame = hh)) +
  geom_point(aes(size = abs(FB))) + geom_vline(xintercept = c(-0.3,0.3), linetype = 2) + 
  theme_bw() + theme(legend.position = "top") + 
  facet_grid(.~dust) + xlim(-0.6,0.6) ->p
ggplotly(p)

Rys. 1 Zmienność FB względem godzin

5.4 GODZINY 1:24 i DNI

b.hourz$hh <- as.numeric(b.hourz$hh)
arrange(b.hourz, dust, typ_sL, dni, hh) -> b.hourz

ggplot(b.hourz %>% filter(dust == "PM10"), aes(y = FB, x = hh, col = typ_sL)) +
  geom_point(size = 2) + geom_line(aes(group = typ_sL)) + labs(x = "hour of day", y = "FB") +
  facet_grid(dni~.) + theme_bw() +theme(legend.position = "top") + 
  geom_hline(yintercept = c(-0.3,0,0.3), linetype = 2) -> p
ggplotly(p)

Rys. 1 Zmienność FB względem godzin i dni

ggplot(b.hourz, aes(y = NMSE, x = FB, shape = typ_sL,col = typ_sL, frame = hh)) +
  geom_point(aes(size = abs(FB))) + 
  geom_vline(xintercept = c(-0.3,0.3), linetype = 2) + 
  theme_bw() + theme(legend.position = "top") + 
  facet_grid(dni~dust) + xlim(-0.6,0.6) -> p
ggplotly(p)

Rys. 1 Zmienność FB względem godzin i dni

5.5 GODZINY CO

ggplot(b.godz, aes(y = FB, x = godz, col = typ_sL)) +
  geom_point(size = 2) + geom_line(aes(group = typ_sL)) +
  facet_grid(dust~.) + theme_bw() +theme(legend.position = "top") + 
  geom_hline(yintercept = c(-0.3,0,0.3), linetype = 2) -> p

ggplotly(p)

Rys. 1 Zmienność FB względem godzin (co 3)

x <- seq(-2,2, by=0.001) ; y <- 4*x^2/(4-x^2) ; yb = 0.45

ggplot(b.godz, aes(y = NMSE, x = FB, shape = typ_sL,col = godz)) +
  labs(x = "Fractional  Bias (FB)", y = "Normalised  Mean  Square  Error (NMSE)", shape = "", col = "") +
  theme_bw() + theme(legend.position = "top") + 
  annotate("polygon",x = c(x[x<=0.3 & x >=-0.3],0.3,-0.3), y = c(y[x<=0.3 & x >=-0.3],yb,yb), 
             alpha = 0.2, fill = "grey", colour = "black") +
  annotate("line",x = x, y = y, colour = "black", linetype = 2) +
  scale_y_continuous(limits = c(0.0, 0.45)) +
  scale_x_continuous(limits = c(-0.5,0.5)) + 
  geom_point(size = 2.5) +
  facet_grid(.~dust)

Rys. 1 Zmienność FB i NMSE względem godzin (co 3)

#ggplotly(p)