---
# System prepended metadata

title: Opis kodu charakterystyk emisji (speed-dependent emission factors)

---

# 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:
![](https://hackmd.io/_uploads/BylGW10wn.png)

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 )
```



