---
# System prepended metadata

title: FAZA 1B_OR (dla OR) - opis kodu (*program pierwszy - główne drogi*)

---

# FAZA 1B_OR (dla OR) - opis kodu (*program pierwszy - główne drogi*)
## wersja dla OR, która uwzględniająca nowopowstałe drogi główne

lokalizacja programu na klastrze:`/mnt/array0/lgawuc/struktura/emisje/obliczenia_e01/liniowka/yanosik/przetwarzanie/Yanosik_odnowa_lech_zima2020/runs/run_43_major_roads_gdal_grid`

skrypt: `rprog_SmartYan_faza_1b_major_roads_v1.R` 


Jest to program, który jest wykonywany jako **pierwszy** w kolejności, działa on dla "dróg głównych", tj. tych o klasach OSM 1 i 2. Klasy 3 i niższe są przetwarzane przez inny program - [link](https://hackmd.io/@yanosik/H1EABxRU2).
poniżej opis komend w kodzie.

---

wczytanie bibliotek i czyszczenie obszaru roboczego:

```r
	require(sf)
	require(dplyr)

	rm(list=ls())
#
# funkcja pomoznicza do zamiany NA na zera
	  fnoDATA <- function (x) {
            x[is.na(x)] <- 0
            return(x)}
```

wczytanie aktualnych danych OSM. Dane muszą zostać wcześniej pobrane z repozytorium. Dane musza pokrywac całą Polskę, więc trzeba ściągnąć wszystkie województwa i przetworzyć do postaci tożsamej z poniższym zbiorem `osm_m_v2_roads.shp`

```r
	osm	<- st_read('/mnt/array0/lgawuc/struktura/emisje/dane_e03/OSM_05sty2023/osm_m_v2_roads.shp')
```

Selekcja autostrad i drog ekspresowych klasy 1 i 2. Tylko te drogi podlegają przetwarzaniu w niniejszym programie. Zamiana nazw klas dróg ze stringów na cyfry i zapis jako `osm_major_roads_for_model_rcl1i2.shp`

```r
        osm    <- filter(osm,osm$fclass %in% c('motorway','motorway_link','trunk','trunk_link')) %>% select(fclass,maxspeed,tunnel)
	osm$fclass	<- gsub('motorway_link',1,osm$fclass)
	osm$fclass	<- gsub('motorway',1,osm$fclass)
	osm$fclass	<- gsub('trunk_link',2,osm$fclass)
	osm$fclass	<- gsub('trunk',2,osm$fclass)
        st_write(osm,'osm_major_roads_for_model_rcl1i2',dsn='.',driver='ESRI Shapefile',append=F)
```
Wczytanie wszystkich istniejacych APRow w PL (w petli po wojewodztwach) i stworzenie jednolitej tablicy dla calego kraju:

```r
	la	<- Sys.glob('/home/lgawuc/yanosik/gooddata/Do_poprawy_geometrii_podzial_wojewodztwa_siatka_ZMAIK/aprs/*shp')
	apr	<- st_read(la[1])	
	for (i in 2:length(la)){
		ab	<- st_read(la[i])	
		apr	<- rbind(apr,ab)
	}

	st_write(apr,layer='apry_all_w_PL',dsn='.',driver='ESRI Shapefile',append=F)
```
Tutaj trzeba wybrac podzbiór APRów - tylko takie punkty które leżą w małej odległości od osi dróg obiektu osm. (ale jest wykomentowane bo trwa zbyt długo)
```r
	#apr_w_buf	<- st_intersection(apr,obf)  - to trwa okropnie dlugo wiec zrobione w Qgis - gdzie jest znacznie szybciej. 
```


Tutaj można wczytać APRy na drogach głównych powstałych z selekcji w QGIS (odleglość od geometrii OSM: 0.0005 deg czyli ok 50 m) i obliczenie trf (rocznego sumarycznego natężenia ruchu w obie strony) oraz frakcji floty (so sc sd bus) jako sumy. Po przerzuceniu tych danych do geometrii OSM z aktualną  siecią dróg będzie założenie ze ruch w obie strony jest taki sam. Flota w obu kierunkach też jest taka sama tak samo jak flota. 
 	obf	<- st_read('apry_all_w_buf_rclAiS_0.0005.shp')
    
        
... ostatecznie to powyższe nie jest zrobione. Prościej jest wczytać APRy z całej PL powstałe wcześniej i zrobić selekcję tylko po OSMclass (nie wykonuję  buforowania w QGIs albo intersekcji w R):

```r	
	obf	<- filter(apr,apr$OSMclass %in% c('1','2'))
```

Obliczenie usrednionych wartosci natezenia ruchu i floty (usrednionych dla dwoch kierunkow)
```r
	obf$trf		<- obf$Psumaro + obf$Rsumaro
	obf$so		<- (obf$Pso    + obf$Rso)/2	
	obf$sd		<- (obf$Psd    + obf$Rsd)/2
	obf$sc		<- (obf$Psc    + obf$Rsc)/2	
	obf$bus		<- (obf$Pbus   + obf$Rbus)/2	
	
	obf	<- select(obf,c('geometry','nrpkt','trf','so','sd','sc','bus'))
```
Zapis punktów APR leżących przy drogach głównych wraz z obliczonymi `trf so sd sc bus`
```r
        st_write( obf , dsn='.',layer='apry_rcl1i2_trf_so_sd_sc_bus',driver='ESRI Shapefile',append=F)
```
Stworzenie wektorów zawierających nazwy oraz indeksy dla floty.
`trf` - sumara pojazdow jakie przejechały dany punkt w ciągu roku (roczne natężenie ruchu)
`so` - samochody osobowe
`sd` - samochody dostawcze
`sc` - samochody ciężarowe
`bus` - autobusy dalekobieżne
Zmienne `so, sd, sc, bus` powinny się sumować do 1 gdyż oznaczają udział poszczególnych rodzajów pojazdów w sumarycznym natężeniu ruchu.

```r
	rstNamVc	<- c('trf','so','sd','sc','bus')
	rstNamVcFull	<- c('1','2','3','4','5')
```
Za pomoca zewnętrznego narzędzia `gdal_grid` stworzenie 5ciu rastrów z przeinterpolowanymi wartościami z APRów. 

**UWAGA**: wynikowa rozdzielczość rastra jest tutaj ustawiona na 1000px. Wyniki modelu były zadowalające (w ocenie wizualnej oraz po obliczeniu sum emisji NOx dla całej Polski, które zostały porównane z sumami EMEP), więc nie widziałem potrzeby aby zwiększać rozdzielczość. Jej wpływ na dokładność wersji SmartYan w fazie 1B oraz 1C jest przedmiotem analizy w ewaluacji wyników modelu SmartYan, która jest przedmiotem powstającej publikacji naukowej. Ewaluacja jest oparta o referecyjne pomiary ZDM APR oraz GDDKiA SCPR. Ewaluacja jest przeprowadzona tylko dla woj. mazowieckiego. 

Poniższa pętla jest realizowana po indeksie `g` dla pięciu elementów `trf so sd sc bus`:

```r
	t <-Sys.time()	
#	for (g in 1:1){
	for (g in 1:5){

        # obiekt z wynikową nazwą zbioru
		rstNam		<- paste0('apry_rcl1i2_trf_so_sd_sc_bus__',rstNamVc[g],'.tif')
		rstNamVcFull[g]	<- rstNam
	
        # stworzenie stringa z komendą dla zewnętrznej funkcji gdal_grid
	        cmd     <- paste0('gdal_grid -zfield ',rstNamVc[g],' -a invdist:power=6.0:smoothing=0.0 -outsize 1000 1000 -a_srs EPSG:4326 -of GTiff  -txe 13.5 24.5 -tye 48.5 55.5 apry_all_w_buf_rclAiS_0.0005_trf_so_sd_sc_bus.shp ',rstNam)
	       	system( cmd, intern = FALSE,ignore.stdout = FALSE, ignore.stderr = TRUE,wait = T, input = NULL) #uruchomienie komendy zewnętrznej
	}#{g}	
        print(Sys.time()-t)
```
Po stworzeniu rastrów należy przepisać ich wartości do geometrii obiektu osm.
```r
    
	osm	<- tibble::rowid_to_column(osm,'rowID')
	osmCt	<- st_centroid(osm) #stworzenie centroidów segmentów dróg
	osmCt	<- osmCt[!st_is_empty(osmCt),] # wyrzucenie nieprawidłowych obiektów
        osmCt   <- dplyr::select(osmCt,c('rowID','geometry'))

	require(raster)
#	for (r in 1:1){
	for (r in 1:5){

        	rst     	<- raster(rstNamVcFull[r])
		rst		<- flip(rst, direction='y') # z jakiegoś powodu dane wychodzą obrucone
	        icrst		<- fnoDATA(raster::extract(rst,osmCt)) # znalezienie wartości pikseli rastra odpowiadające
	        osmCt$icrst	<- round(icrst,2) # zaokrąglenie 
		st_write(osmCt,layer=paste0('osmCt_',rstNamVcFull[r],'.shp'),dsn='.',driver='ESRI Shapefile',append=F)

        # połączenie danych wyekstrahowanych z rastra z danymi przestrzennymi - obiekt "osm" gdzie wcześniej dopisana została kolumna "rowID"
		df_osmCt	<- osmCt %>% st_drop_geometry()
		names(df_osmCt)	<- c('rowID',rstNamVc[r])
		osm		<- st_as_sf(left_join(osm,df_osmCt,by='rowID'))
	} #{r}
```
Przeliczenie natezenia ruchu na polowe (jeden segment Yanosik to ruch w jedna strone). W tej wersji programu liczone sa tylko drogi najwyzszego rzedu a one sa jednokierunkowe wiec trf musi byc podzielone przez 2

 **UWAGA**: dane sciagniete wprost z repozytorium OSM sa jednokierunkowe, mimo ze rzeczywistosci ruch odbywa się  w dwóch kierunkach. informacja o tym znajduje sie w kolumnie bodaj 'oneway' (w danych OSM ściągniętych z repozytorium, w danych Yanosik tej kolumny nie ma). Dane Yanosik dostarczone przez firme Neptis maja rowniez segmenty jednokierunkowe (tj. jeden segment drogi, stanowiący jeden wers w zbiorze danych jest zawsze jednokierunkowy). Jednakże segmenty dróg w danych Yanosik są zduplikowane jeżeli dana droga jest dwukierunkowa w rzeczywistosci.
Tj. jeden segment ma geometrie o punkt zalamania A B C D E F, a drugi kierunek ruchu ma punkty zamałania F E D C B A. Wyjaśnienie z ilustracjami tutaj [link](https://hackmd.io/@yanosik/HkVR9wmIn)

Kolejny krok to obliczenie skalibrowanego nateżenia ruchu. Wartość `300` jest typowo przyjmowaną wartością w inżynierii ruchu dla przeskalowania natężenia ruchu dobowego na sumę roczną. Jest to związane z faktem, iz wszystkie modele ruchu sumulują ruch dobowy na max 24 h do przodu i w inżynierii ruchu nie zajmują się ruchem rocznym. Wartość `300` zostala zasugerowana w drodze licznych konsulacji przez prof. Andrzeja Szaratę z Politechniki Krakowskiej.
```r

	osm$trfBulk	<- osm$trf
	osm$trf24h	<- osm$trf/300/2
	osm$trf		<- osm$trf/2
```
Zapis wynikowego zbioru danych, który jest potem brany do obliczeń w programie drugim [link](https://hackmd.io/@yanosik/H1EABxRU2).
```r
      st_write(osm,layer='yanosik_by_SmartYan_faza_1b__major_roads_rcl1i2_interpolGGrid',dsn='.',driver='ESRI Shapefile',append=F)
```
---
koniec zbioru `rprog_SmartYan_faza_1b_major_roads_v1.R`
	
