# 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`
Sign in with Wallet
Connect another wallet