Korona nas je sve naterala da razmišljamo o odlasku u prirodu i o kupovini vikendica. Iako ne planiramo ni žena ni ja da kupujemo vikendicu u dogledno vreme (jel ste videli vi te cene!?), to nije razlog da ne razmišljam o tome gde bi ona bila i kako bi izgledala. Ono što bih voleo je neko fino mesto na što većoj nadmorskoj visini, ali na žalost – u okolini Beograda teško da tako nešto postoji (možda Genex). Počeo sam da razmišljam kako analitički da nađem mesto za vikendicu tako što ću da izračunam gde bi to tačno bilo!

TL;DR

Neću da smaram nikoga da ovo sve čita ako ga ovo ne zanima samo da sazna pobednike, pa evo kratkog izdanja.

  • Ukoliko se gleda apsolutno sve, pobednik je Košutnjak.
  • Ukoliko se gleda malo van Beograda, a bilo koje mesto, čak iako ne može da se gradi, onda je pobednik vrh na Avali.
  • Ukoliko se gleda nešto iole realno, prvo mesto dele Lipovica i Kopaonik, a drugo Divčibare i vikend naselja na Avali.

Ideja

OK, zadatak je prost. Prvo mi treba mapa Srbije sa visinama i zatim mapa Srbije sa udaljenostima od Beograda. Onda, da bih našao što veću nadmorsku visinu, a što bliže Beogradu, podelim ta dva broja za svaki piksel na karti i eto – nađem najoptimalnije mesto! Time dobijam mapu “metar/minut” gde svaki piksel predstavlja ovu imaginarnu metriku, ali nije ni bitna tačna jedinica mere, već su bitni relativni odnosi ovih vrednost za celu Srbiju. Žena mi je rekla da sam zaludan i zašto da samo prosto ne kupimo nešto na Divčibarama, ali sam osećao da mogu bolje od toga! Da li sam uspeo u tome – nastavite da čitate da biste saznali! (a ako nešto možda naučite o QGIS-u ili se zainteresujete za kartografiju, tim još bolje:)

Topo mapa Srbije

Ovo je bio lak posao. Samo treba naći najbolju topografsku mapu Srbije. Najbolje besplatno što postoji danas su satelitski snimci Kopernikus satelita (pa dobro, nije skroz besplatno, poreski obveznici u EU ovo omogućuju:), gde se nude DEM podloge sa 1 arc-second rezolucijama (20-tak metara) koji mogu da se skinu sa https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1. Kad se skine GeoTIFF od 5GB odavde, bolje da se odseče samo na Srbiju (za odsecanje skinite granice Srbije odavde) čime dobijemo nešto ovako (klik za uvećanje):

Jel samo ja vidim guštera u Vojvodini:)

Za one koji ne znaju – ovo je rasterizovana karta Srbije, gde svaki piksel nosi neku vrednost od 0m do 2169m (Midžor). Da se razumemo, ovaj osnovni crno-beli GeoTIFF prikaz jeste ružan, ali ako se doda bolji color ramp gradijent i stavi se “Blending mode” na “Darken”, dobijamo ovako malo lepše:

Što južnije, to… svežije?

OK, prvi deo posla gotov, idemo na teži deo – udaljenosti od Beograda.

Mapa udaljenosti od Beograda

Na početku sam rekao da nam treba mapa udaljenosti od Beograda, ali sam tu mislio na izohrone. Postoje izohipse, izobate, izobare… šta su sad izohrone? Izohrone (iso = isto, chronos = vreme) su linije koje na mapi povezuju sve tačke do kojih se dolazi u isto vreme. Prepuštam da za detalje pogledate vikipediju. Jako zanimljiva oblast urbanističkog planiranja, taksija, špedicija…, ali i dosta zajebana (da ne kažem “tehnički izazovna”:D). Mislio sam da ću prosto da nađem neku gotovu mapu izohrona, ili da ću samo da iskoristim neki postojeći servis za ovo baziran na OSM-u, ali ne lezi vraže. Svi servisi koje sam probao, a bilo ih je mnogo (ORS, GraphHopper, Bing, HERE, iso4app, …) ili nisu mogli da obrade celu Srbiju ili ne daju granularnosti koji mi trebaju! Ograničenja su uglavnom u kratkim vremenskim intervalima ili površini koja može da se obuhvati. Postoji ovaj sajt koji čak može da ovo izračuna, ali nažalost, meni je trebalo vreme dolaska za SVAKU tačku Srbije, a ne samo da dobijem poligone za 30/60/90/120… minuta. BTW, evo šta ovaj sajt k(l)aže dokle sve može da se stigne iz Beograda za 5h:

Brže u Zagreb, nego u Bosilegrad

Možda sam mogao i da startujem OSRM lokalno, ali mislim da bi mi trebalo mnogo vremena da ga podesim (pa bar mnogo više nego što sam bio spreman da potrošim na ovo:) A onda sam skontao da OSMnx (softver za rad sa ulicama iz OSM-a) ima mogućnost da računa najkraći put između dve tačke! A to je upravo ono što mi je trebalo – da učitam OSM ulice Srbije, da ih pretvorim u graf (jer šta je rutiranje nego nalaženje puta kroz graf:D), da granama u grafu dodelim težinu (vreme potrebno za prolaz kroz ulicu, tj. granu) i da za svaki čvor (spoj dve ulice) izračunam vreme dolaženja do njega. Deluje kao nešto komplikovano, ali OSMnx rešava većinu ovih stvari, a ima čak i primer sa izohronama. Ceo posao se sveo na sledeći kod:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
# Load whole street network in Serbia
# Driving only, this is lot of data already
graph = ox.graph_from_place('Serbia', network_type='drive')
 
# Choose node from which we will calculate travel times from (Novi Beograd)
origin_node = ox.nearest_nodes(graph, 20,42321, 44,81834)
 
# Add time it takes to cross each street/edge
graph = ox.project_graph(graph)
graph = ox.add_edge_speeds(graph, fallback=50)
graph = ox.add_edge_travel_times(graph)
 
# Since basic graphs do not know anything about traffic jams, traffic lights etc.
# I am increasing these for 1.5 time (empirically chosen constant:D)
for _, _, _, data in graph.edges(data=True, keys=True):
    data["travel_time"] = data["travel_time"] * 1.5
 
# Calculate how much time it takes to all other nodes from origin_node
# and write that to all nodes
res, _ = nx.single_source_dijkstra(graph, origin_node, weight='travel_time')
for node_id, data in graph.nodes(data=True):
    data["travel_time"] = res[node_id] if node_id in res else -1
 
ox.save_graph_geopackage(graph, 'serbia.gpkg')

Par pojašnjenja oko ovog koda. OSMnx je jako zahtevan memorijski i za Srbiju treba minimum 64GB da se iskonstruište ceo graf. Takođe, OSMnx nema mogućnost da računa “turn restriction-e”, tj. ne zna baš-baš tačan put da izračuna, ali je ovo što je izračunao dovoljno dobro u praksi. Takođe, sve dobijene vrednosti sam pomnožio sa 1.5 da kompenzujem stvarni svet – OSMnx računa teorijska vremena, bez semafora, gužvi, čekanja… Za divno čudo, lupio sam konstantu 1.5, ali se ispostavila dosta precizna! No, ovo i ne služi da bude savršena rasterska karta izohrona, već dobra aproksimacija; čak i da je sve pomereno za 2h – nije bitno, pošto su nam bitne samo relativne vrednosti.

Sa ovim smo dobili samo tačke na mapi gde se ulice spajaju i vremena dolaženja do tih ulica. To kada se učita u QGIS izgleda nekako ovako prosto:

Sve tačke u grafu obojene po vremenu dolaska iz Novog Beograda

Na ovome mora još da se radi jer su neke deonice jako duge (npr. autoputevi) i tačke su onda retke i dalje računanje izohrona postaje jako nelinearno i isprekidano. Takođe ovako ima dosta rupa (pogledajte Vojvodinu, to su sve putevi, ali nema čvorova jer nema presecanja puteva). Zato sam odlučio da iz QGIS-a isparčam sve linije na manje celine od po 200m, uz pomoć funkcije “Points along geometry”, a onda da tako dobijenim tačkama nađem nova vremena za vreme dolaženja tako što interpoliram vreme dolaženja od početne do kranje tačke, proporcionalno daljini nove tačke od početne tačke, u odnosu na dužinu celog puta. Uf, na primeru će biti možda lakše. Npr. ako je u početnu tačku puta dolazak u 1200s, a u krajnju tačku u 1800s, ceo put je dužine 200m, a nova tačka je na 50m od početne, onda će vreme dolaska u novu tačku biti 1350s (50m:200m=(1350s-1200s):(1800s-1200s) proporcija). Ovo sam uradio tako što sam napravio novi atribut koji se izračunava sa “scale_linear” funkcijom (najjača funkcija ikada, ne znam kako je nema u standardnim bibliotekama u svim jezicima):

1
2
3
4
5
6
7
scale_linear(
	attribute($currentfeature, 'distance'), // distanca
	0, // početna dužina
	attribute($currentfeature, 'length'), // krajnja dužina
        attribute(get_feature('layer', 'osmid', attribute($currentfeature, 'from')), 'travel_time'),
	attribute(get_feature('layer', 'osmid', attribute($currentfeature, 'to')), 'travel_time')
)

Pomoću ove magije, dobio sam nove čvorove sa interpoliranim vremenima dolazaka. Da bi vizuelno razumeli o čemu se radi, evo slika pre i posle dodavanja novih čvorova:

Prvo se vide čvorovi i vremena samo tamo gde se ulice spajaju u grafu, a posle se vidi i gomila usputnih čvorova. Oznake na čvorovima su sekunde, a ova lokacija je inače ovde.

Sad, kada je ovo odrađeno, treba samo ove tačke prebaciti u rastersku formu. Za to sam iskoristio QGIS funkciju “Raster→Analysis→Grid (Inverse distance to a power)” sa podrazumevanim vrednostima i voilà – evo rasterske mape izohrona u Srbiji:

Jel nazirete “žile” od autoputeva?

Huh, nije bilo naivno, ali sam naprčkao nešto. Ima ovde još mesta za poboljšanje. Na primer:

  • prvo je što sam koristio “Inverse distance to a power” što znači da tamo gde nema tačaka, da se tamo vrednosti smanjuju, a treba da se povećavaju! Teši me što je takvih mesta malo, pa nije problem. Rešenje je verovatno da se daljine prebace u 1/x, da se onda rasteriziju, pa vrati opet preko 1/x funkcije
  • drugo je što imamo rupe po Srbiji tamo gde nema ulica – idealno bi bilo da sva mesta po Srbiji bez ijednog puta označimo nekako (npr. napravi se pravougaona mreža tačaka), i da se tim tačkama dodele neka vremena dolaska (npr. peške od najbližeg puta/tačke)
  • treće je nekako smanjiti probleme tačaka koje su blizu, a imaju užasno različita vremena – npr. tačka na autoputu na jednu stranu kaže da je vreme dolaska 1h, a u suprotnom smeru (zbog okretanja na najbližoj petlji) bude 1.5h. Za to nemam neko dobro rešenje

No, sve u svemu, zadovoljan sam rezultatom – nije savršeno, ali će raditi posao za traženje vikendice:)

Čisto da se vidi koliko je ovo ispalo dobro – pogledajte koliko malo izolovanih “ostrva” ima (izolovana ostrva ne smeju da postoje na ovakvim mapama jer bi to značilo da u neko ostrvo može da se stigne pre nego u okolinu tog ostrva!) i koliko su izohrone precizne (1h razlike na 6h, u odnosu na HERE izohrone pomenute iznad):

Znam da će prvi komentari da budu “e, ovo nije tačno zato što do mesta X ima Y minuta, a tebi piše Z” ili “ja brate kad spičim do X, to ne može da se prikaže na ovoj tvojoj karti”, ali opet – poenta ove karte nije da bude precizna izohrona, nego da predstavi relativna vremena dolaženja do svih krakova Srbije!

Sad samo podelimo ove dve mape i…

Da rezimiramo – imamo:

  • rastersku topo mapu sa vrednostima piksela od 0m do 2169m, i
  • rastersku mapu izohrona sa vrednostima od 0h do 6.5h

Takve dve mape možemo da podelimo iz QGIS-a koristeći “Raster→Raster Calculator…” funkciju:

Evo koliko je ovo prosto u QGIS-u

I dobijemo nešto ovako:

Vrednosti idu od 0 do 16 (kada se podele metri sa minutama; 16 je npr. 1000m za 62.5 min). Već i sa ovakve slike se vide pojedini delovi koji se ističu, vidi se Beograd da se ističe, kao i Divčibare, Kopaonik, Suva i Stara Planina i Beljanica. Sva sreća pa se niko nije setio da ide ili kupuje vikendice na ovim mestima! Neće da izađe na dobro, ovoliko muke da bi na kraju žena bila u pravu. Možda ako uradimo izohipse, možda ćemo videti malo bolje i preciznije. Da bismo našli izohipse, treba odabrati “Raster→Extraction→Contour”. Stavimo da je interval na “1” (da su izohipse svakih 1 metar/minuta), isfiltriramo sve ispod 15 i posmatramo šta nam je preostalo. Pobednik je ipak:

Ma daaaaaj…

Dobra vest je da žena ne može da mi kaže “jel sam ti lepo rekla!”. Loša vest je da nešto nema vikendica u ponudi na Košutnjaku:) Priča se ovih dana o nekim finim investitorima koji bi da grade neke vikendice po Košutnjaku, malo veće vikendice doduše i za više porodica i od betona tako da liče na zgrade, ali se nadam da im neće proći:)

Ako pogledamo sve preko vrednosti 9 metara/minut, sve je i dalje u Beogradu. Ima tu i Avale, ali na ovom delu nema vikendica koliko ja znam. Primetićete da malo više stvari ima po Novom Beogradu pošto mi je početna tačka stavljena tamo:

Ispod 8 je slika ovakva nekakva:

Sa slike se ne vidi najbolje, pa evo top pobednika (ako ne računamo mesta po Beogradu):

Deo Metara/minut Visina Komentar
Avala vrh 11 511m Ovde ne može da se prave vikendice
Rudnik 8 1132m Samo Cvijićev vrh, isto mala površina
Stara planina 8 2916m Samo greben oko Midžora, isto nema gradnje
Lipovica 7 ~300m Prvo realno mesto gde ima logike kupiti vikendicu
Mali i Srednji Povlen 7 1346m Selo Taor, samo vrhovi
Beljanica 7 1339m Samo vrh
Šiljak 7 1560m Samo vrh
Trem 7 1810m Samo par vrhova na Suvoj planini
Kopaonik 7 >1800 Bogami, ceo centar Kopaonika je ovde
Avala vikend naselja 6 ~250m Vikend naselja Šuplja stena i Karagača
Divčibare 6 >950m Centar Divčibara i okolina, mogu vikendice
Kosmaj 6 626m Samo vrh
Bukulja 6 696m Samo vrh
Rudnik 6 >750m Ne znam Rudnik dobro, ali deluje da ovde možda može da bude vikendica

Sve u svemu, ako ne gledamo vrhove, prvo mesto dele Lipovaja i Kopaonik, a drugo Divčibare i vikend naselja na Avali.

Umesto kraja

Nadam se da vam je ova analiza bila zanimljiva kao i meni dok sam je pravio! Teško je držati pažnju ljudima, ali se nadam da ste bar bacili pogled na slike:) Sve fajlove možete da skinete odavde. Ostalo mi je nekih stvari koje nisam završio (npr. malo preciznije računanje izohrona ili da probam da dam veću težinu visini, pa da bude metar/minut2), ali negde crta mora da se podvuče. A ženi mogu samo da kažem “eto jednom da i ti ne budeš u pravu”:)