-
Notifications
You must be signed in to change notification settings - Fork 0
/
00_NO2_Grdz_Confound.R
244 lines (234 loc) · 9.37 KB
/
00_NO2_Grdz_Confound.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
# Beispiel Heizg_daten
library(xts)
library(tidyverse)
library(lubridate)
library(broom)
library(gridExtra)
options(digits = 3)
load("~/Documents/Luftqualitaet/Analysen/NO2/Heizg_data.RData")
load("~/Documents/Luftqualitaet/Analysen/BW_station.RData")
BW_stations_NO2_tbl %>% summary()
BW_stations_NO2_tbl$name %>% levels() # 15 Stationen
Heizungsdaten_15_19 %>% head() # Daten Bad Cannstatt, Bernhausen, Reutlingen
Heizungsdaten_15_19 %>% summary()
NO2_data <- BW_stations_NO2_tbl %>% filter(name %in% c("Can","Brn","Rt_"))
stat_urban <- c("Can","Brn","Rt_")
stat_trafic <- c("Lbg_Friedr","Rtl","Nck" )
NO2_trafic <-BW_stations_NO2_tbl %>% filter(name %in% stat_trafic)
Grdz_Can <- Heizungsdaten_15_19 %>%
.[["Can"]] %>% mutate(name = "Can")
Grdz_Brn <- Heizungsdaten_15_19 %>%
.[["Brn"]] %>% mutate(name = "Brn")
Grdz_Rt <- Heizungsdaten_15_19 %>%
.[["Rt"]] %>% mutate(name = "Rt_")
Grdz_data <- bind_rows(Grdz_Brn,Grdz_Can,Grdz_Rt)
summary(Grdz_data)
NO2_Grdz_data <- inner_join(NO2_data,Grdz_data, by = c("name","datetime"))
NO2_Grdz_data$NO2 <- NO2_Grdz_data$NO2 %>% na.locf()
NO2_Grdz_data$Gradzahl<- NO2_Grdz_data$Gradzahl %>% na.locf()
cor(NO2_Grdz_data$NO2,NO2_Grdz_data$Gradzahl)
NO2_Grdz_data$name <- as_factor(NO2_Grdz_data$name)
NO2_Grdz_data$station <- as_factor(NO2_Grdz_data$station)
summary(NO2_Grdz_data)
# Vergleich Verteilung mit Normalverteilung
NO2_Grdz_data %>%
ggplot() +
stat_qq(aes(sample=NO2),distribution = stats::qnorm) +
facet_wrap(~name)
NO2_Grdz_data %>%
ggplot() +
geom_histogram(aes(NO2), bins = 40) +
ggtitle("Haeufigkeitsverteilung der NO2 Messungen
Poissonverteilung")+
facet_wrap(~name)
NO2_Grdz_data %>%
ggplot() +
geom_freqpoly(aes(NO2), bins = 40) +
ggtitle("Haeufigkeitsverteilung der NO2 Messungen
Poissonverteilung")+
facet_wrap(~name)
range(NO2_Grdz_data$NO2) # 0 bis 176
NO2_Grdz_data %>% filter (Gradzahl > 5) %>%
ggplot() +
geom_freqpoly(aes(Gradzahl), bins = 40) +
ggtitle("Haeufigkeitsverteilung des Heizbedarfs
dargestellt als Gradzahlen")+
facet_wrap(~name)
range(NO2_Grdz_data$datetime)#"2015-01-01 01:00:00 UTC" "2019-10-14 23:00:00 UTC"
NO2_Grdz_data %>% filter (Gradzahl > 5) %>%
ggplot() +
geom_point(aes(Gradzahl,NO2),size = 0.1, alpha = 0.5) +
geom_smooth(method = "lm", aes(Gradzahl,NO2),col = "red")+
ggtitle("NO2 als Funktion des Heizbedarfs
im städtischen Hintergrund",
subtitle = "Regressionsgerade (rot)" )+
labs(x= "Heizbedarf (Gradzahl)", y = " NO2 [μg/m3]")+
facet_wrap(~name)
NO2_Grdz_data %>% ggplot(aes(datetime,Gradzahl))+
geom_smooth(method = "lm",aes(datetime,Gradzahl))+
geom_point(aes(datetime,Gradzahl),size = 0.01,alpha = 0.01)+
ggtitle ("Heizwaermebedarf als Gradzahl
einschl. Energie f. WW-bedarf")+
facet_wrap(~ name)
#======================================
# Verkehrsnahe Stationen
NO2_trafic <- NO2_trafic %>% filter (datetime > ymd("2015-01-01"))
NO2_trafic %>% ggplot()+
geom_point(aes(datetime,NO2),size = 0.01,alpha= 0.1) +
geom_smooth(method = "lm",aes(datetime,NO2),col = "red")+
labs( x = "",y = "NO2 [μg/m3]")+
ggtitle(" NO2 1-h Messwerte verkehrsnahe Stationen
Ludwigsburg Friedr.Str,Stgt Am Neckartor,
Reutlingen Lederstrasse",
subtitle = "Trend als Regressionsgerade (rot)")+
facet_wrap(~name)
NO2_trafic$station <- NO2_trafic$station %>% as_factor()
#====================
#Umrechnung Jahreswerte
td <- difftime(ymd_h("1970-01-01 00"),ymd_h("2015-01-01 00"),
units="secs")
jahrs <- 365*24*60*60
# use broom to visualise model
NO2_Grdz_data %>% group_by(name) %>% as.tbl() %>%
do(tidy(lm(NO2 ~ datetime,data =.),conf.int= TRUE)) %>%
filter(term == "datetime") %>%
select(name,estimate,conf.low,conf.high)%>%
mutate(estimate = estimate*jahrs,
conf.low = conf.low*jahrs,
conf.high= conf.high*jahrs) %>%
ggplot(aes(name,y = estimate,ymin =conf.low,ymax= conf.high))+
geom_errorbar()+
geom_point()+
ggtitle("Jahresrate Veränderung NO2 - Immissionen",
subtitle = "Hintergrund: Bernhausen,Cannstatt,Reutlingen")+
labs(y = "NO2 [μg/m3]pro Jahr", x = "Messstation")
NO2_Grdz_data %>% group_by(name) %>% as.tbl() %>%
do(tidy(lm(Gradzahl ~ datetime,data =.),conf.int= TRUE)) %>%
filter(term == "datetime") %>%
select(name,estimate,conf.low,conf.high)%>%
mutate(estimate = estimate*jahrs,
conf.low = conf.low*jahrs,
conf.high= conf.high*jahrs) %>%
ggplot(aes(name,y = estimate,ymin =conf.low,ymax= conf.high))+
geom_errorbar()+
geom_point()+
ggtitle("Jahresrate 15/19 Veränderung Heizbedarf",
subtitle = "Hintergrund: Bernhausen,Cannstatt,Reutlingen")+
labs(y = " Gradzahl[°C ] pro Jahr", x = "Messstation")
# NO2 Trafic regression
NO2_trafic$NO2 <- NO2_trafic$NO2 %>% na.locf()
NO2_trafic %>% group_by(name) %>% as.tbl() %>%
do(tidy(lm(NO2 ~ datetime,data =.),conf.int= TRUE)) %>%
filter(term == "datetime") %>%
select(name,estimate,conf.low,conf.high)%>%
mutate(estimate = estimate*jahrs,
conf.low = conf.low*jahrs,
conf.high= conf.high*jahrs) %>%
ggplot(aes(name,y = estimate,ymin =conf.low,ymax= conf.high))+
geom_errorbar()+
geom_point()
df <-bind_rows(NO2_trafic %>% group_by(name) %>% as.tbl() %>%
do(tidy(lm(NO2 ~ datetime,data =.),conf.int= TRUE)) %>%
filter(term == "datetime") %>%
select(name,estimate,conf.low,conf.high),NO2_Grdz_data %>% group_by(name) %>% as.tbl() %>%
do(tidy(lm(NO2 ~ datetime,data =.),conf.int= TRUE)) %>%
filter(term == "datetime") %>%
select(name,estimate,conf.low,conf.high))%>%
mutate(estimate = -estimate*jahrs,
conf.low = -conf.low*jahrs,
conf.high= -conf.high*jahrs) %>% as_tibble()
df <-arrange(df, desc(estimate))
df[,1] <- c("1 Nck","2 Lbg_friedr","3 Can", "4 Brn")
df %>% ggplot(aes(name,y = estimate,ymin =conf.low,ymax= conf.high))+
geom_errorbar()+
geom_point()+
ggtitle("Jahresrate 15/19 Red.NO2 Immissionen",
subtitle = "Städt. Hintergrund/ verkehrsnah")+
labs(y = "NO2 [μg/m3]pro Jahr", x = "Messstation") +
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
df %>% ggplot()
ggplot()
grid.table(df)
levels(BW_stations_NO2_tbl$name) <-BW_stations_NO2_tbl$name %>% levels()%>%
recode("Lbg_4" ="Lbg_Weimar","Rt_" = "Rt_pomol","Rtl" ="Rt_leder",
"Friedri"= "Friedrichshafen","Heid" ="Heidelbg","Heil"="Heilbron")
BW_stations_NO2_tbl$datetime %>% range()
df_15 <-BW_stations_NO2_tbl %>% group_by(name) %>%
do(tidy(lm(NO2 ~ datetime,data =.),conf.int= TRUE)) %>%
filter(term == "datetime") %>%
select(name,estimate,conf.low,conf.high)%>%
mutate(estimate = estimate*jahrs,
conf.low = conf.low*jahrs,
conf.high= conf.high*jahrs)%>%
as_tibble() %>%
arrange(., estimate)
# for visualisation compare datasciecebook page 185
df_15 %>% ggplot(aes(x = reorder(name,desc(estimate)),
y = estimate,ymin =conf.high,ymax= conf.low),label = name)+
geom_errorbar()+
geom_point(show.legend = TRUE)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
#coord_flip()+
ggtitle("Jahresrate NO2 -Immissionen
Reduzierung 2000 bis 2020")+
labs(y = "NO2 [μg/m3]pro Jahr", x = "Messstation")
BW_stations_NO2_tbl <- BW_stations_NO2_tbl %>%
mutate(NO2 = na.locf(NO2),station =as_factor(station)) %>%
select(name = "name",NO2 = "NO2",datetime,station)
# Verkehrsnahe Stationen
stat_trafic <- c("Nck","Rtl","Lbg_Friedr")
stat_rural <- c("Alb","Odw","Sws")
BW_stations_NO2_tbl %>%
filter(name %in% stat_trafic) %>%
ggplot(aes(x=datetime,y = NO2))+
geom_smooth(col = "red")+
geom_smooth(method= "lm",linetype = 2, col = "darkred")+
#geom_point(size = 0.00001, alpha = 0.5)+
facet_wrap(~name)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs( x = "", y= "NO2[μg/m3]")+
ggtitle( " NO2- Immissionen
Trends und Mittelwerte",
subtitle = " verkehrsnahe Messungen
Lbg.Friedrichstr.,Am Neckartor,Rtl.Lederstr. ")
BW_stations_NO2_tbl %>%
filter(name %in% c("Can","Brn","Rt_")) %>%
ggplot(aes(x=datetime,y = NO2))+
geom_smooth(col = "red")+
geom_smooth(method= "lm",linetype = 2, col = "darkred")+
#geom_point(size = 0.00001, alpha = 0.5)+
facet_wrap(~name)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs( x = "", y= "NO2[μg/m3]")+
ggtitle( " NO2- Immissionen
Trends und Mittelwerte",
subtitle = " städt. Hintergrund-Messungen
Stgt Bad Cannstatt,Bernhausen,Rtl.Pomologie ")
levels(BW_stations_NO2_tbl$name)
BW_stations_NO2_tbl %>%
filter(name %in% stat_urban) %>%
ggplot(aes(x=datetime,y = NO2))+
geom_smooth(col = "red")+
geom_smooth(method= "lm",linetype = 2, col = "darkred")+
#geom_point(size = 0.00001, alpha = 0.5)+
facet_wrap(~name)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs( x = "", y= "NO2[μg/m3]")+
ggtitle( " NO2- Immissionen
Trends und Mittelwerte",
subtitle = " städt. Hintergrund-Messungen
Stgt Bad Cannstatt,Bernhausen,Rtl.Pomologie ")
levels(BW_stations_NO2_tbl$name)
BW_stations_NO2_tbl %>%
filter(name %in% c("Alb","Sws")) %>%
ggplot(aes(x=datetime,y = NO2))+
geom_smooth(col = "red")+
geom_smooth(method= "lm",linetype = 2, col = "darkred")+
#geom_point(size = 0.00001, alpha = 0.5)+
facet_wrap(~name)+
theme(axis.text.x = element_text(angle = 90, hjust = 1))+
labs( x = "", y= "NO2[μg/m3]")+
ggtitle( " NO2- Immissionen
Trends und Mittelwerte",
subtitle = " ländl. Hintergrund-Messungen
Schw.-Alb,Schwarzwald Süd")