# Opis kodu charakterystyk emisji (speed-dependent emission factors)
lokalizacja:`/mnt/array0/lgawuc/struktura/emisje/obliczenia_e01/liniowka/yanosik/przetwarzanie/Yanosik_odnowa_lech_zima2020/multicores/emisje_licz/model_v1.0`
zbiór: `input_wielomiany_v8_charakterystyki_za2021__modified_apr23.R`
**UWAGA:** w poniższym kodzie zależność pomiędzy rodzajami pojazdów a klasą drogi jest na sztywno wpisana (zaadoptowana z dawnego modelu IMKR z POlitechniki Warszawskiej/GDDKiA) - obiekt `flo`. Jeżeli model wykorzystywany do określenia aktywności akceptowalnie poprawnie estymuje flotę to nie należy korzystać ze sztywnych danych z `flo`.
---
Na początku zdefiniowana jest funkcja która zawiera charaksterystyki/wielomiany. Opisuje jej cześć, gdyż to za dużo tekstu, który sie powtarza dla kolejnych zanieczyszczeń:
```r
wielom <- function( Vso, Vsd, Vsc, Vbus,
ilso, ilsd, ilsc, ilbus, nazwy , druk) {
ef <- as.data.frame(array(,dim=c(length(Vso),1)))
e <- as.data.frame(array(,dim=c(length(Vso),1)))
ef$COso <- ( 7.9414180424E-10*Vso^5 - 2.8881963157E-07*Vso^4 + 3.8421033952E-05*Vso^3 - 1.9438978031E-03*Vso^2 + 6.8801144142E-03*Vso+ 1.7125454624E+00)*ilso
ef$COsd <- ( -5.4310511658E-07*Vsd^3 + 3.3991919879E-04*Vsd^2 - 4.0592774345E-02*Vsd+ 1.5252174178E+00)*ilsd
ef$COsc <- ( 8.1393534101E-12*Vsc^6 - 3.9468784895E-09*Vsc^5 + 7.7123271141E-07*Vsc^4 - 7.8094272260E-05*Vsc^3 + 4.3996139970E-03*Vsc^2 - 1.3666391884E-01*Vsc+ 2.3256016166E+00)*ilsc
ef$CObus <- ( 1.8353675551E-11*Vbus^6 - 9.0958820933E-09*Vbus^5 + 1.8297780893E-06*Vbus^4 - 1.9239584911E-04*Vbus^3 + 1.1355480736E-02*Vbus^2 - 3.7204202968E-01*Vbus+ 6.3692345594E+00)*ilbus
ef$NMVOCso <- (9.313297136813470E-11*Vso^5 - 3.624515022943500E-08*Vso^4 + 4.748740375941800E-06*Vso^3 - 1.915353825019680E-04*Vso^2 - 5.319750398590240E-03*Vso + 4.098243750833050E-01)*ilso
# ( -2.8402608461E-09*Vso^4 + 4.0310795026E-07*Vso^3 + 5.5255293678E-05*Vso^2 - 1.1069224421E-02*Vso+ 4.5035366922E-01)*ilso - to daje ujemne waetosci powyzej 140 kmh
ef$NMVOCsd <- ( 6.6500835014E-13*Vsd^6 - 2.6135617236E-10*Vsd^5 + 3.6491970597E-08*Vsd^4 - 2.0820566034E-06*Vsd^3 + 4.6475430219E-05*Vsd^2 - 1.7841427148E-03*Vsd+ 1.1858358263E-01)*ilsd
ef$NMVOCsc <- ( 1.0932238916E-12*Vsc^6 - 5.4784258489E-10*Vsc^5 + 1.1030990064E-07*Vsc^4 - 1.1452033742E-05*Vsc^3 + 6.5672673876E-04*Vsc^2 - 2.0609189137E-02*Vsc+ 3.2781645831E-01)*ilsc
ef$NMVOCbus <- ( 2.6704553445E-13*Vbus^6 - 3.1460065007E-10*Vbus^5 + 1.0675432005E-07*Vbus^4 - 1.6458346556E-05*Vbus^3 + 1.3022807331E-03*Vbus^2 - 5.2979274162E-02*Vbus+ 1.0237450946E+00)*ilbus
```
..i tak dalej. Vs* to parametry prędkości a il* to parametry liczby poszczególnych rodzajów pojazdów.
ale są "kwiatki" wiec dopisuje funkcję EmBus, która koryguje wartości dla autobusów:
```r
ef$NH3so <- ( -4.5135998406E-08*Vso^3 + 1.1408295390E-05*Vso^2 - 6.0141154405E-04*Vso+ 2.0942289533E-02)*ilso
ef$NH3sd <- ( -6.8879834154E-09*Vsd^3 + 1.8338964164E-06*Vsd^2 - 1.1844229458E-04*Vsd+ 6.9170237755E-03)*ilsd
ef$NH3sc <- ( 1.6940658945E-20*Vsc^2 - 2.0735366549E-18*Vsc+ 8.6488706314E-03)*ilsc
EmBus <- function(x,il){ if(x<=50){( 2.7105054312E-20*x^2 - 3.1441863002E-18*x+ 4.6866300412E-03)*il} else {(2.7105054312E-20*x^2 - 5.8546917314E-18*x+ 4.6866300412E-03)*il}}
ef$NH3bus <- mapply(Vbus,ilbus,FUN=EmBus)
```
...i tak dalej
dla BaP sa fukcje prostsze:
```r
ef$BaPso <- ( -5.4256735506E-24*Vso+ 8.3839652317E-07)*ilso
ef$BaPsd <- ( 3.2554041304E-23*Vsd+ 1.6375085470E-06)*ilsd
ef$BaPsc <- ( -5.4256735506E-24*Vsc+ 5.2282233413E-07)*ilsc
ef$BaPbus <- ( -5.4256735506E-24*Vbus+ 5.2282233413E-07)*ilbus
```
dalej zsumowanie emisji dla zanieczyszczen:
```r
e$CO <- ef$COso + ef$COsd + ef$COsc + ef$CObus
e$CO2 <- ef$CO2so + ef$CO2sd + ef$CO2sc + ef$CO2bus
e$SO2 <- ef$SO2so + ef$SO2sd + ef$SO2sc + ef$SO2bus
e$NOx <- ef$NOxso + ef$NOxsd + ef$NOxsc + ef$NOxbus
e$PM10 <- ef$PM10so + ef$PM10sd + ef$PM10sc + ef$PM10bus
e$PM25 <- ef$PM25so + ef$PM25sd + ef$PM25sc + ef$PM25bus
e$PMEX <- ef$PMEXso + ef$PMEXsd + ef$PMEXsc + ef$PMEXbus
e$NMVOC <- ef$NMVOCso + ef$NMVOCsd + ef$NMVOCsc + ef$NMVOCbus
e$BaP <- ef$BaPso + ef$BaPsd + ef$BaPsc + ef$BaPbus
e$NH3 <- ef$NH3so + ef$NH3sd + ef$NH3sc + ef$NH3bus
e$CH4 <- ef$CH4so + ef$CH4sd + ef$CH4sc + ef$CH4bus
e$Hg <- ef$Hgso + ef$Hgsd + ef$Hgsc + ef$Hgbus
e$Ni <- ef$Niso + ef$Nisd + ef$Nisc + ef$Nibus
e$Cd <- ef$Cdso + ef$Cdsd + ef$Cdsc + ef$Cdbus
e$Pb <- ef$Pbso + ef$Pbsd + ef$Pbsc + ef$Pbbus
e$As <- ef$Asso + ef$Assd + ef$Assc + ef$Asbus
e <- dplyr::select(e,-(V1))
ef <- dplyr::select(ef,-(V1))
```
Udostępenienie zmiennej e w tle:
```robotframework
assign('e' ,e, envir = .GlobalEnv )
```
Tutaj mozna zrobic wykresy za pomocą zmiany wartości zmiennej `druk`. Tzn. tutaj jest kod do wykresow a zmiane wartości trzeba zrobić przy wywoływaniu funkcji `wielom :`
```r
if (druk == TRUE ) {
ef$Vr <- Vso
for ( k in 1:ncol(ef)){
png(paste0('wykres_',names(ef)[k],'.png' ))
plot(ef[,'Vr'],ef[,k],ylab=paste0(names(ef)[k],' [g/kilometr]'),xlab='Vr')
abline(v=0)
dev.off()
}
for ( k in 1:ncol(e)){
png(paste0('wykres_zbiorczy',names(e)[k],'.png' ))
plot(ef[,'Vr'],e[,k],ylab=paste0(names(e)[k],' [g/kilometr]'),xlab='Vr')
abline(v=0)
dev.off()
}
}
} #koniec funkcji {wielom}
```
Tutaj jest program zasadniczy.
Stworzenie wektorów prędkości dla różnych pojazdów:
```r
wVso <- as.vector(array(1,dim<-c(17,1)))
wVsd <- as.vector(array(1,dim<-c(17,1)))
wVsc <- as.vector(array(1,dim<-c(17,1)))
wVbus <- as.vector(array(1,dim<-c(17,1)))
fso <- as.vector(array(1,dim<-c(17,1)))
fsd <- as.vector(array(1,dim<-c(17,1)))
fsc <- as.vector(array(1,dim<-c(17,1)))
fbus <- as.vector(array(1,dim<-c(17,1)))
```
Wczytanie danych dotyczących występowania rodzajów pojazdów w zależności od klasy drogi:
```r
flo <- read.table('/home/lgawuc/yanosik/multicores/emisje_licz/model_v1.0/flota_z_punktow_imkr.csv',sep=' ',head=T)
```
Ten zbiór jest utworzony na podstawie danych z modelu ruchu. Ale to już stary zbiór a model ruchu był z Politechniki Warszawskiej a nie z CUPT. Mozna zaktualizować. ale... skoro model F1B i F1C i F2a maja dane o flocie o można to usunąć... lub wręcz trzeba. Moje niedopatrzenie.
zbior wygląda tak:
```r
rcl so sd scp sc bus
1 0.68847 0.08188 0 0.22546 0.0042
2 0.75634 0.08302 0 0.15515 0.00548
3 0.75886 0.07935 0 0.15185 0.00993
4 0.83955 0.07906 0 0.07178 0.00928
5 0.85928 0.07357 0 0.05042 0.01661
6 0.79738 0.09326 0 0.09888 0.01048
7 0.84765 0.0887 0 0.05464 0.00893
8 0.87934 0.07665 0 0.03582 0.00772
9 0.79897 0.07812 0 0.1141 0.00864
```
dalej kod dla pozostałych klas dróg (10-17):
```r
flo[10:14,2:6] <- flo[1:5,2:6]
flo[10:14,1] <- 10:14
flo[15:17,2] <- 1
flo[15:17,3:6] <- 0
flo[15:17,1] <- 15:17
for (i in 1:17) {
fso[i] <- flo[i,'so']
fsd[i] <- flo[i,'sd']
fsc[i] <- flo[i,'sc']
fbus[i] <- flo[i,'bus']}
```
Dalej następuje stworzenie wektorów służących do obliczenia prędkości poszczególnych rodzajów pojazdów na podstawie średniej prędkości wszystkich użytkowników. Wartości zapisane poniżej pochodzą z materiałów znajdujących się na chomiku autorstwa prof. Chłopka
```r
\\chomik2\CBE\pozyskane materiały\GDDKiA\prof_Chlopek
```
Znajduje się tam tabela:
która posłużyła do obliczenia poniżych wektorów dla przelicznia prędkości.
wyniki są dostępne w lokalizacji: `/mnt/array0/lgawuc/struktura/emisje/obliczenia_e01/liniowka/yanosik/przetwarzanie/Yanosik_odnowa_lech_zima2020/multicores/emisje_licz`
zbiór `srednia-.xlsx`
**UWAGA** Ze względu na wiek tego przeliczenia należałoby zrobić aktualizację, gdyż wektory "wVso" są przyjęte subiektywnie i należałoby je zastąpić lepszymi, np. opartymi o dane ZMR CUPT.
```r
# rcl <- 1 ; 10 (A autostrady ; OSM: motorway)
wVso[1] <- 1.285855589 ; wVsd[1] <- 1.107814045 ; wVsc[1] <- 0.81372549
wVbus[1] <- 0.820969337 ; wVso[10] <- 1.285855589 ; wVsd[10] <- 1.107814045
wVsc[10] <- 0.813725491 ; wVbus[10]<- 0.820969337
# rcl <- 2 ; 11 (S ; OSM: trunk)
wVso[2] <- 1.242236025 ; wVsd[2] <- 1.118012422 ; wVsc[2] <- 0.827160494 ; wVbus[2] <- 0.832298137
wVso[11] <- 1.242236025 ; wVsd[11] <- 1.118012422 ; wVsc[11] <- 0.827160494 ; wVbus[11] <- 0.832298137
# rcl <- 3 ; 12 (GP ; OSM: primary)
wVso[3] <- 1.223021583 ; wVsd[3] <- 1.079136691 ; wVsc[3] <- 0.857142857 ; wVbus[3] <- 0.863309353
wVso[12] <- 1.223021583 ; wVsd[12] <- 1.079136691 ; wVsc[12] <- 0.857142857 ; wVbus[12] <- 0.863309353
# rcl <- 4 ; 13 (G ; OSM: secondary)
wVso[4] <- 1.222410866 ; wVsd[4] <- 1.103565365 ; wVsc[4] <- 0.8438818 ; wVbus[4] <- 0.848896435
wVso[13] <- 1.222410866 ; wVsd[13] <- 1.103565365 ; wVsc[13] <- 0.8438818 ; wVbus[13] <- 0.848896435
# rcl <- 5 ; 14 (Z ; OSM: terriary)
wVso[5] <- 1.226993865 ; wVsd[5] <- 1.083844581 ; wVsc[5] <- 0.85279187 ; wVbus[5] <- 0.858895706
wVso[14] <- 1.226993865 ; wVsd[14] <- 1.083844581 ; wVsc[14] <- 0.85279187 ; wVbus[14] <- 0.858895706
# rcl <- 6 (D ; OSM: unclassified)
wVso[6] <- 1.192546 ; wVsd[6] <- 1.068322 ; wVsc[6] <- 0.8696 ; wVbus[6] <- 0.8695652
# rcl <- 7 (L ; OSM: residental)
wVso[7] <- 1.1023622 ; wVsd[7] <- 1.0078740 ; wVsc[7] <- 0.944881 ; wVbus[7] <- 0.94488189
```
Powiązanie z danymi o prędkości `vr` wszystkich użytkowników:
```r
Vso <- iYan$vr*wVso[iYan$rcl]
Vsd <- iYan$vr*wVsd[iYan$rcl]
Vsc <- iYan$vr*wVsc[iYan$rcl]
Vbus <- iYan$vr*wVbus[iYan$rcl]
```
Napisanie danych o flocie (prawdopodobnie do wyrzucenia, bo jak pisałem wyżej dane z F1C f1B i F2a zawierają już te dane..):
```r
iYan$so <- fso[iYan$rcl]
iYan$sd <- fsd[iYan$rcl]
iYan$sc <- fsc[iYan$rcl]
iYan$bus <- fbus[iYan$rcl]
```
Przepisanie predkosci pojazdow do tablicy `iYan` i stworzenie tablic potrzebnych do funkcji `wielom`
```r
iYan$Vso <- Vso ; ilso <- users * iYan$so * r
iYan$Vsd <- Vsd ; ilsd <- users * iYan$sd * r
iYan$Vsc <- Vsc ; ilsc <- users * iYan$sc * r
iYan$Vbus <- Vbus ; ilbus <- users * iYan$bus * r
```
Odpalenie funkcji zdefiniowanej wyżej:
```r
wielom( Vso, Vsd, Vsc, Vbus, ilso, ilsd, ilsc, ilbus, nazwy=n, druk=FALSE )
```
Mozna tutaj wpisac sztuczne dane aby zobaczyc jak wygladaja wykresy emisji w zalenosci od predkosci (ale jest to wykomentowane obecnie):
```r
# zrobienie wykresow kontrolnych dla charakterystyk emisji
# Vso <- 1:140 ; ilso <- 1
# Vsd <- 1:140 ; ilsd <- 1
# Vsc <- 1:140 ; ilsc <- 1
# Vbus <- 1:140 ; ilbus <- 1
#
# wielom( Vso, Vsd, Vsc, Vbus, ilso, ilsd, ilsc, ilbus, nazwy=n, druk=T )
```
Sign in with Wallet
Connect another wallet