# FAZA 1B_OR (dla OR) - opis kodu (*program drugi - pozostałe 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/quantiled_Yan_APRs_data`
zbiór: `rprog_SmartYan_faza_1b_minor_roads_v2.R`
Jest to program, który jest wykonywany jako **drugi** w kolejności, działa on dla "pozostałych dróg", tj. tych o klasach OSM 3 i niższych. Klasy 1 i 2 są przetwarzane przez inny program - [link](https://hackmd.io/@yanosik/BkB0MNALn).
poniżej opis komend w kodzie.
---
Podział klas dróg i odpowiadającym im liczbom wg danych Yanosik
```r
# klasy drog OSM oraz odpowiadajace im numery
# MOTORWAY; 1 autostrady
# TRUNK; 2 drogi ekspresowe
# PRIMARY; 3 drogi krajowe
# SECONDARY; 4 wojewodzkie? ponizej podzial nie jest juz taki scisly
# TERTIARY; 5
# UNCLASSIFIED; 6
# RESIDENTIAL; 7
# LIVING_STREET; 8
# SERVICE; 9
# MOTORWAY_LINK; 10
# TRUNK_LINK; 11
# PRIMARY_LINK; 12
# SECONDARY_LINK; 13
# TERTIARY_LINK; 14
# ROAD; 15
# TRACK; 16
```
wyczyszczenie obszaru roboczego:
```r
rm(list=ls())
require(sf)
require(dplyr)
```
definicja przydatnych później funkcji:
```r
fnoDATA <- function (x) {
x[is.na(x)] <- 0
return(x)}
`%ni%` <- Negate(`%in%`)
```
wczytanie podstawowego zbioru danych Yanosik dla calej PL i ujednolicenie nazewnictwa klas dróg:
```r
yb <- st_read('/mnt/array0/lgawuc/struktura/emisje/obliczenia_e01/liniowka/yanosik/przetwarzanie/Yanosik_odnowa_lech_zima2020/gooddata/base_yanosik/2017_yanosik_baseData_UUhrSUM_SYGhrSUM.shp')
rcl <- yb$road_class
rcl <- ifelse(rcl=='10','1',rcl)
rcl <- ifelse(rcl=='11','2',rcl)
rcl <- ifelse(rcl=='12','3',rcl)
rcl <- ifelse(rcl=='13','4',rcl)
rcl <- ifelse(rcl=='13','5',rcl)
yb$rcl <- rcl
```
Ograniczenie danych yanosik - wyrzucenie klas 1 i 2 (major roads) bo one sa zrobione osobno. Major roads 1 i 2 sa mocno zmienne w ostatnich latach dlatego jest to zrobione osobno nie na danycn Yanosik lecz na geometrii OSM wprost (link)
```r
ys <- filter(yb,yb$rcl %ni% c('1','2')) %>% select(-road_class)
```
```r
#
# wyrzucenie ze zbioru APRow tych ktore maja klase 1 i 2 (jw)
aprPL <- st_read('/home/lgawuc/yanosik/runs/run_43_major_roads_gdal_grid/apry_all_w_PL.shp')
aprPL <- filter(aprPL,aprPL$OSMclass %ni% c('1','2'))
aprPL$trf <- aprPL$Psumaro + aprPL$Rsumaro
ap <- aprPL
ap$so <- (ap$Pso + ap$Rso)/2
ap$sd <- (ap$Psd + ap$Rsd)/2
ap$sc <- (ap$Psc + ap$Rsc)/2
ap$bus <- (ap$Pbus + ap$Rbus)/2
```
ta wersja modelu (SmartYan faza 1b) jest oparta o statystyczna kwantyzacje (grupowanie) danych wedlug decyli
tutaj tworze wektor tych decyli
```r
qs <- seq(0,1,0.1)
```
petla sluzaca do stworzenia 10 zbiorow APR dla 10 decyli:
`przAPR` - dla danych CUPT/PPR/APR
`przYan` - dla danych Yanosik
```r
for (i in 1:10){
# for (i in 1:2){
przYan <- filter(ys, (ys$uu_hr > as.numeric(quantile((ys$uu_hr ),qs[i])) & ys$uu_hr < as.numeric(quantile((ys$uu_hr ),qs[i+1]))) )
yanNam <- paste0('2017_yanosik_baseData_UUhrSUM__Qu',qs[i],'_do_',qs[i+1])
st_write(przYan,layer=yanNam ,dsn='.',driver='ESRI Shapefile',append=F)
aprNam <-paste0('APRs_allPL__Qu',qs[i],'_do_',qs[i+1])
przAPR <- filter(ap, (ap$trf > as.numeric(quantile((ap$trf ),qs[i])) & ap$trf < as.numeric(quantile((ap$trf ),qs[i+1]))) )
st_write(przAPR,layer=aprNam,dsn='.',driver='ESRI Shapefile',append=F)
}
```
Stworzenie wektorow 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')
```
Definicja nazwy ostatecznego zbioru wynikowego i usunięcie już istniejacego (bo potem jest `append=T`)
```r
osmOutNam <-p aste0('yanosik_by_SmartYan_faza_1b__minor_roads_interpolGGrid')
unlink(paste0(osmOutNam,'.*'))
```
za pomoca funkcji `gdal_grid` stworzenie rastrow z przeinterpolowanymi wartosciami z APRow.
Tutaj są dwie pętle po `i` i po `g`. Pętla po `i` jest dla decyli (dziesięć kroków). Pętla po `g` jest dla rocznego natężenia ruchu `trf` oraz frakcji floty.
Niektóre części kodu są wykomentowane (np. zapis do pliku) gdyż służą sprawdzaniu poprawności wyników. Model działa prawidłowo więc nie trzeba zapisywać wszystkich kroków iteracji na dysk. Zostawiam wykomentowane bo może się przydać w dalszych pracach z programem.
```r
require(raster)
t <-Sys.time()
for (k in 1:10){
# for (k in 1:1){
aprNam <- paste0('APRs_allPL__Qu',qs[k],'_do_',qs[k+1])
aprNamShp <- paste0('APRs_allPL__Qu',qs[k],'_do_',qs[k+1],'.shp')
yanNam <- paste0('2017_yanosik_baseData_UUhrSUM__Qu',qs[k],'_do_',qs[k+1])
yanNamShp <- paste0('2017_yanosik_baseData_UUhrSUM__Qu',qs[k],'_do_',qs[k+1],'.shp')
# wczytanie danych Yanosik dla poszczegolnych decyli
osm <- st_read(yanNamShp)
osm <- tibble::rowid_to_column(osm,'rowID')
# zamiana na centroidy bo funkcja raster::extract szybciej dziala jak sie wezmie punkty
osmCt <- st_centroid(osm)
# wyrzucenie pustych geometrii
osmCt <- osmCt[!st_is_empty(osmCt),]
osmCt <- dplyr::select(osmCt,c('rowID','geometry'))
message(' ')
print(paste0('new quantile: ',k,' ',qs[k],'_do_',qs[k+1] ))
message(' ')
# for (g in 1:1){
for (g in 1:5){
rstNam <- paste0( aprNam,'__gdalGridded_',rstNamVc[g],'.tif')
# komenda do odpalenia gdal_grid - zewnetrzna
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 ',aprNamShp,' ',rstNam)
system( cmd, intern = FALSE,ignore.stdout = FALSE, ignore.stderr = TRUE,wait = T, input = NULL)
#
# po stworzeniu rastrow nalezy przepisac ich wartosci do pustej geometrii obiektu osm
rst <- raster(rstNam)
rst <- flip(rst, direction='y') # z jakiegos powodu rastry wychodza jakies odwrocone
icrst <- fnoDATA(raster::extract(rst,osmCt)) #zastosowanie wcześnie zdefniowanej funkcji fnoDATA
osmCt$icrst <- round(icrst,2)
#st_write(osmCt,layer=paste0('osmCt_',rstNamVcFull[r],'.shp'),dsn='.',driver='ESRI Shapefile',append=F) - mozna zapisac albo nie ale to sluzy tylko do sprawdzenia
df_osmCt <- osmCt %>% st_drop_geometry()
names(df_osmCt) <- c('rowID',rstNamVc[g])
# tutaj jest problem z klasa obiektu ale da sie to zrobic przez st_as_sf
osm <- st_as_sf(left_join(osm,df_osmCt,by='rowID'))
} #{g}
```
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 <- round(osm$trf,1)
osm$trf24h <- round(osm$trf/300/2,1)
osm$trf <- round(osm$trf/2,1)
st_write(osm,layer=osmOutNam,dsn='.',driver='ESRI Shapefile',append=T)
# osmAllQu <- bind_rows(osmAllQu,osm)
}#{k}
```
na tym etapie bilbioteka raster nie jest potrzebna (a zawiera duplujące się funkcje z dplyr) więc ją wyłączam
```r
detach(package:raster,unload=TRUE)
print(Sys.time()-t) # obliczenie czasu jaki został skonsumowany na obliczenia w pętli po `i` oraz `g`
```
Na tym etapie drogi pozostałe (poza rcl 1 i 2) są już skalibrowane, więc można zrobić polaczenie danych z przetworzonymi wcześniej (program pierwszy) major roads policzonymi wczesniej.
gdzie `major` oznacza wynikony zbiór danych wygenerowany przez program pierwszy [link](https://hackmd.io/@lgawuc/BkB0MNALn).
```r
major <- st_read('/mnt/array0/lgawuc/struktura/emisje/obliczenia_e01/liniowka/yanosik/przetwarzanie/Yanosik_odnowa_lech_zima2020/runs/run_43_major_roads_gdal_grid/yanosik_by_SmartYan_faza_1b__major_roads_rcl1i2_interpolGGrid.shp') %>% dplyr::select(-rowID)
names(major) <- c('rcl','maxspeed','tunnel','trf','so','sd','sc','bus','trfBulk','trf24h','geometry')
major <- major %>% select(rcl,maxspeed,trf24h,so,sd,sc,bus,geometry)
major$uu_hr <- -999
major$v_sr_rok <- -999
major$nr_segm <- -999
osm <- st_read(paste0(osmOutNam,'.shp'))
osm$maxspeed <- -999
osmMerge <- osm %>% select(uu_hr,v_sr_rok,nr_segm,rcl,trf24h,maxspeed, so,sd,sc,bus,geometry)
fout <- bind_rows(osmMerge,major)
unlink('yanosik_by_SmartYan_faza_1b__minor_AND_major_roads.*')
st_write(fout,layer="yanosik_by_SmartYan_faza_1b__minor_AND_major_roads",dsn='.',driver='ESRI Shapefile',append=F)
```
Ostatni etap to uruchomienie programu który dzieli powstałe dane na poszczególne województwa:
```r
cmd <- paste0('nohup R CMD BATCH rprog_cut_po_woj_v1.R')
system( cmd, intern = FALSE,ignore.stdout = FALSE, ignore.stderr = TRUE,wait = FALSE, input = NULL)
```
___
koniec zbioru `rprog_SmartYan_faza_1b_minor_roads_v2.R`
Sign in with Wallet
Connect another wallet