Pre izvesnog vremena smo uzeli Mladenovac i izmapirali sve kuće u njemu. Postavilo se pitanje koju bismo sledeću opštinu mogli da mapiramo. Bilo je logično uzeti neku gde je mapiranost zgrada najgora, ali kako to naći?

U ovom postu ću pričati malo kako provaliti gde je mapiranost mala, kako Python-om procesirati geojson i Shapefile-ove i pokušati da dođem do odgovora na pitanje koja je opština u Srbiji najgore mapirana.

Šac metoda

Najprostiji i najočigledniji način je da se ode na sajt disaster.ninja, da se odabere odnos zgrada i naseljenosti i da se vidi „šac” metodom šta će da ispadne. Evo kako izgleda cela Srbija (primetiti kako se Mladenovac zeleni😀):

Nažalost, „šac” metoda se razlikuje od nivoa zuma i neki gradovi deluju loše mapirani, ali onda zumiraš i vidiš da se ipak donekle OK urađeni. Evo npr. kako izgleda Jagodina iz daljine i kada se zumira:

Daleko da je savršena, ali vide se zeleni (mapirani) delovi. To nije ni blizu kao kad se pogleda Mladenovac, gde kada smo počeli da ga radimo, bilo je ukupno 2 (slovima: dve) zgrade na celoj opštini.

Obilazeći Srbiju ovako na ovom sajtu dolazimo do sledećih „kandidata”:

• Grad Sombor (ne cela opština)
• Opština Bosilegrad
• Opština Osečina
• Opština Bela Crkva
• Opština Despotovac
• Opština Ćićevac
• Opština Surdulica
• Naselje Bečej

Da li su ovo najbolji kandidati ili ne – teško je reći. Ali zato možemo da koristimo analitiku😉.

Analitika

Evo šta mi je palo na pamet. Ovaj disaster.ninja prikazuje odnos zgrada iz OSM-a i gustine naseljenosti i to po nekim šestaugaonicima, ali još bolje bi bilo da možemo da imamo odnos stvarnog broja zgrada i zgrada iz OSM-a, i to po opštinama – tako bi znali koliko i gde još stvarno fali. Međutim, nemamo prave zgrade (da imamo – problem možda ne bi ni postojao😀), ali imamo nešto što je možda blizu tome, a to su zgrade koje je izgenerisao Microsoft uz pomoć mašinskog učenja i dao na otvoreno i besplatno korišćenje. Zvaćemo ih ML zgrade na dalje. Koji mi je plan – 1. uzmemo zgrade iz OSM-a, 2. uzmemo ove ML zgrade, 3. nađemo njihov odnos za svaku opštinu, sortiramo to od najboljeg do najgoreg i tako nađemo sledeću opštinu za mapiranje! Ako vas ne zanima programiranje, preskočite na deo „Rezultati” da vidite šta se dobija.

Da se bacimo na posao – od sastojaka nam trebaju:

• ML zgrade koje uzmemo odavde (nađite Srbiju tu i skinućete ZIP u kome su zgrade u „geojsonl” formatu),
• OSM zgrade ćemo da uzmemo odavde (skinite .shp.zip u kome je jedan od Shapefile-ova sa zgradama iz OSM-a),
• Administrativne granice koje je najlakše skinuti odavde (kliknite „Download” i odaberite Srbiju sa svim granicama od nivoa 2 do 9 i dobićete geojson sa svim granicama)

Šta nam sad treba – da otvorimo ovaj geojson sa granicama (iz poslednjeg koraka) i da dodamo nove kolone, npr. „osm_buildings” i „ml_buildings” u kojima ćemo da upišemo koliko smo zgrada našli za svaku opštinu. Da otvorimo zgrade, dovoljno je u Python-u da uradimo ovo:

1
2

import geopandas as gpd
boundaries = gpd.read_file('serbia-boundaries.geojson')

Tako dobijamo sve granice u GeoPandas data frame-e. Da bismo otvorili Shapefile sa OSM zgradama, koristićemo Fiona biblioteku.

1
2

import fiona
osm_buildings = fiona.open("gis_osm_buildings_a_free_1.shp", "r")

E, sad iteriramo za svaku OSM zgradu i nađemo sve preseke sa administrativnim granicama koje smo uzeli kao u prethodnom koraku (hoćemo da dobijemo sva naselja i opštine koje zgrada geometrijski seče). GeoPandas ima podršku za preseke geometrija, pa ni ovo nije teško:

1
2
3
4
5

from shapely.geometry import shape
 
for osm_building in osm_buildings:
    building_shape = shape(osm_building['geometry'])
    filtered_boundaries = boundaries.sindex.query(building_shape, predicate='intersects')

Primetite da moramo da konvertujemo ulaznu zgradu u „shape”, i da koristimo „sindex” (tj. „spatial index” da bismo ubrzali operaciju presecanja zgrade sa granicama). E, sad, operacija „query” vraća indekse granica (tj. indekse pozicija tih objekata u data frame-u). Tehnički, mogli bismo da odmah napravimo novu kolonu „osm_buildings” i da to odmah vraćamo nazad u data frame tako što inkrementiramo tu kolonu na indeksu svake presečene zgrade, ali bi to bilo jako sporo, pa sam se odlučio da sve te vraćene indekse čuvam u privremenom dictionary-ju (običnom Python) koji ćemo tek na kraju da prebacimo u data frame:

1
2

for idx in filtered_boundaries:
    counter[idx] += 1

A na kraju uradimo ovako nešto i snimimo novi geojson:

1
2
3

import pandas as pd
boundaries['osm_buildings'] = pd.Series(counter)
boundaries.to_file('boundaries_with_osm.geojson', driver='GeoJSON', encoding='utf-8')

Što se tiče ML zgrada, situacija je slična – iteriramo za sve zgrade i dodamo ih kao novu kolonu. Razlika je ovog puta što čitamo „geojsonl” format (to je geojson gde je jedan objekat u jednoj liniji u fajlu) i što je fajl oooogroman, pa sam odlučio da ne koristim ove fensi Python biblioteke (GeoPandas) jer su još nezrele i ne mogu da strimuju fajl, nego moraju da ga učitaju celog, a ove ML zgrade su ogromne i ovo ne bi radilo ko nema par gigabajta RAM-a. Iskoristio sam činjenicu da je samo jedna zgrada u liniji ulaznog fajla i čitao sam onda liniju po liniju (najgluplje moguće, na standardni Python način) i tako obrađivao zgrade. Sve u svemu, ceo kod izgleda nekako ovako:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51

import geopandas as gpd
import pandas as pd
import fiona
import json
from shapely.geometry import shape
 
 
def append_ml_buildings(boundaries):
    counter = {}
    for idx in boundaries.index:
        counter[idx] = 0
 
    with open('serbia-ml-buildings.geojsonl', 'r') as building_file:
        for line in building_file:
            building_shape = shape(json.loads(line)['geometry'])
            filtered_boundaries = boundaries.sindex.query(building_shape, predicate='intersects')
            for idx in filtered_boundaries:
                counter[idx] += 1
 
    return counter
 
 
def append_osm_buildings(boundaries):
    counter = {}
    for idx in boundaries.index:
        counter[idx] = 0
 
    with fiona.open("gis_osm_buildings_a_free_1.shp", "r") as osm_buildings:
        for osm_building in osm_buildings:
            building_shape = shape(osm_building['geometry'])
            filtered_boundaries = boundaries.sindex.query(building_shape, predicate='intersects')
            for idx in filtered_boundaries:
                counter[idx] += 1
 
    return counter
 
 
def main():
    boundaries = gpd.read_file('serbia-boundaries.geojson')
    _ = boundaries.sindex
 
    counter = append_ml_buildings(boundaries)
    boundaries['ml_buildings'] = pd.Series(counter)
 
    counter = append_osm_buildings(boundaries)
    boundaries['osm_buildings'] = pd.Series(counter)
    boundaries.to_file('output.geojson', driver='GeoJSON', encoding='utf-8')
 
 
if __name__ == '__main__':
    main()

Rezultati

Interesantno, ali neki delovi Srbije uopšte nisu imali ML zgrade – celo Trgovište, jedan deo oko Pirota i jedan mali deo posle Priboja na tromeđi sa BiH i Crnom Gorom. Tamo definitivno ima zgrada na satelitskim snimcima, tako da ne znam šta je posredi. Ako uzmemo odnos OSM i ML zgrada i probamo da iscrtamo sva naselja („admin_level=9” u OSM terminologiji), dobijamo ovakvu situaciju (klik na veću sliku):

Ovo je dosta nepregledno jer gledamo jako sitne teritorije, ali mogu da se primete neke stvari. Beograd i Užice su sjajni, Vojvodina ne zaostaje mnogo, kao ni delovi oko većih gradova. Ali – mi gledamo da mapiramo teritoriju slične veličine kao opština Mladenovac, pa je onda logičnije da gledamo opštine („admin_level=8” u OSM terminologiji) i gradove, a ne naselja. Kada to vizuelizujemo, dobijamo ovako nešto (klik za veću sliku):

Kad sam radio ovu sliku, bio sam malo mudriji, pa sam stavio i „100+ %” jer ima mesta gde je broj zgrada u OSM-u veći nego broj ML zgrada (što je skroz OK i legitimno). Opet se vidi kako Mladenovac dominira, ali se vidi i par „rupa”. Ali koja je rupa najgora? Na kraju je data cela tabela sa svim opštinama, pa analizirajte! Par mojih komentara:

• Rača je nešto manja od Mladenovca, nema velikih gradova i deluje kao dobar kandidat
• Ovi moji kandidati dobijeni „šac” metodom (Sombor, Bosilegrad…) nisu ni u prvih 10 (osim Ćićevca), toliko o šac metodama😀
• U moju, i u odbranu svih šac metoda – ja sam prirodno po sajtu disaster.ninja gledao veće gradove, pa su mi oni bili na vrhu, tako da mi je analiza bila nagnuta (engl. biased) na tu stranu. Ovde su uglavnom pri vrhu opštine koje nemaju baš velike gradove u svojoj blizini (jer opštine sa velikim gradovima obično imaju nešto bar mapirano), tako da je i ova analiza nagnuta na drugu stranu. Možda nije loše ni tako odabrati opštinu!
• Kad smo kod toga i kod „impakta” – prema Vikiju, Mladenovac ima 50,000 duša, Rača 10,000, a Sombor 80,000, tako da bi definitivno bilo i teže i veći impakt da se odradi Sombor.
• Sve vizuelizacije lažu, pa tako i ova iznad. Kada je čovek vidi, pomisli da je Srbija relativno dobro mapirana, ali to je samo zato što sam tako odabrao granice opsega. Da su mi granice opsega bile zaokružene (0-20%, 20-40%…), otkrila bi se mnogo realnija situacija izmapiranosti zgrada u Srbiji, ovako nešto:

Još jedan pokazatelj kako sa podacima i vizuelizacijama može da se manipuliše ljudima.

Na kraju dana, treba odabrati:

• Sombor – slične površine kao Mladenovac, sa više stanovnika, veći impakt, ali teže će ga biti završiti. Nije centralna Srbija, pa posle Mladenovca malo šaramo i u Vojvodini.
• Kikinda – manje stanovnika nego Sombor, ali isto u Vojvodini, a najgore mapiran grad u Vojvodini
• Jagodina – najgore mapiran grad na nivou cele Srbije
• Rača – mala opština, najmanji apsolutni odnos izmapiranih zgrada u Srbiji, manji impakt, ali možemo je relativno brzo završiti.

Na kraju, evo tabele svih ML i OSM zgrada za sve opštine: