Poligon típusú shapefile rétegek hálózattá alakítása R nyelven

Csakúgy mint számos egyéb esetben, ezt a blogposztot is egy aktuálisan felmerült kutatási probléma inspirálta. Én mostanában a hitbizományokkal foglalkozom. Ezeknek a nagybirtokoknak az volt a sajátossága, hogy egy meghatározott rokoni körön belül öröklődtek, vagyis nem volt szabad – sem egészben, sem részben – idegeneknek eladni azokat. A 19-20. század fordulójának közéleti diskurzusában gyakori téma volt az a feltételezés, hogy a hitbizományok leuralnak bizonyos vidékeket az országban, ahonnan a kivándorlásba, egykézésbe, agrárszocialista szervezkedésbe menekül a földszerzés lehetőségétől elzárt, kilátástalan helyzetben lévő paraszti népesség.

De hogyan lehetne meghatározni a hitbizományokat által belakott vidékeket? A közigazgatási egységek – vármegyék, járások – nem feltétlenül jelentenek releváns kereteket egy társadalomtörténeti vizsgálathoz. Ezek túlságosan nagyok és a rájuk kiszámolt középértékek elfedik az olyankor igen jelentős belső különbségeket. Földrajzi távolságokkal szintén nem dolgozhatunk. Hiszen egy-egy település (olykor igen kiterjedt) határán belül nem tudható, hogy pontosan hol feküdtek a hitbizományként lekötött parcellák. Így arra gondoltam, hogy megpróbálkozom az érintett települések egymáshoz viszonyított helyzetét figyelembe venni. Olyan területeket keresek az országban, ahol összefüggő tömböt vagy láncolatot alkottak a hitbizományi birtokok, vagyis ahol sok, egymással szomszédos településen találunk ilyeneket.

Mivel itt egy nagyobb adathalmaz vizsgálatáról van szó, ezért szükség van egy olyan módszerre, amellyel automatizáltan megragadható a “szomszédos település” fogalma. Erre kínál lehetőséget a hálózatos megközelítés.

A felépíteni tervezett hálózatban Magyarország települései jelentik a csomópontokat. A linkek pedig azokat a helységeket kötik össze, amelyeknek a határa fizikailag is érintkezik egymással. Az így létrehozott hálózatban egy konkrét település Ego-networkjének lekérdezésével megkaphatjuk azokat a további településeket, amelyek közvetlenül szomszédosak azzal. De lehetőség kínálkozik például a kisebb-nagyobb összefüggő részhálózatok (komponensek) megfigyelésére, vagy a hálózatos értelemben vett távolságok kiszámítására is. Szóval egy olyan eszközköz juthatunk, amelynek segítségével az egymáshoz viszonyított helyzetük alapján tudjuk megragadni és jellemezni a települések egy bizonyos, számunkra éppen érdekes körét.

A hálózat alapanyaga

Az 1910-es népszámlálás szerint az akkori Magyarország – Fiume és Horvátország nélküli – területén 12.541 település volt található. Azt hiszem magától értetődő, hogy ekkora adattömegnél a hálózat manuális létrehozása fel sem merülhet. Szükség van tehát egy olyan alapanyagra, ami valamilyen módon tartalmazza a hálózat elkészítéséhez szükséges alapvető információkat (csomópontok, kapcsolatok). És egyébként a számítógép által is értelmezhető formátumban van. Erre kiválóan alkalmasnak látszik a Demeter Gábor által vezetett GISta Hungarorum OTKA projekt (2015–2017) keretei között létrehozott községhatáros digitális alaptérkép.

Amit itt látunk, az egy úgynevezett poligon típusú shapefile réteg. Ez azt jelenti, hogy az ábrázolás alapegysége az önmagába visszatérő vonal által határolt felület. Jelen esetben 12.541 településünk van, amelyeket ugyanennyi poligon szimbolizál. (A shapefile egy vektoros adatformátum, amit térinformatikai alkalmazásokhoz fejlesztettek ki.)

Egy poligon típusú térképi réteg létrehozásának menete általában az, hogy a térinformatikai szoftverbe betöltenek egy nagy felbontású szkennelt térképet, majd az ábrázolni kívánt térképi egységeket körbe rajzolják egy másik rétegen. Nagyjából úgy, mintha egy átlátszó fóliát helyeznének az eredeti térképre.

Ha jobban belenagyítunk a térképbe, akkor megfigyelhetjük, hogy a poligonok körvonala töréspontokból és a köztük lévő egyenes vonalakból áll. Az alábbi ábrán a töréspontokat csillagokkal ábrázoltam. (Ezek normál esetben természetesen nem ilyen hangsúlyosak.)

Mivel a földrajzi tér minden pontjának tartoznia kell valahová, ezért a településeket szimbolizáló poligonok szorosan érintkeznek egymással, közös határvonalaik, azokon pedig közös töréspontok vannak.

A szomszédos települések hálózatának létrehozásához a fenti térkép azon tulajdonságát fogom kihasználni, hogy egy-egy töréspont több településhez is tartozhat, vagyis amely településeknek közös töréspontjaik vannak, azok biztosan szomszédosak egymással.

Az elképzelésem az, hogy lekérdezem az egyes települések határvonalain lévő töréspontokat. Ezt követően minden települést összehasonlítok minden másik településsel, és ahol közös töréspontokat találok, ott szomszédosnak tekintem azokat.

A hálózat előállítása

A feladat végrehajtásához elsőként be kell szereznünk a kérdéses községhatáros térképet. Ezt a GISta Hungarorum oldalán az ide belinkelt virtuális mappából tudjuk megtenni. Ebben több különböző, más-más célra szolgáló shapefile réteg található. Nekünk az MO_Telepules nevű fájlokra(!) van szükségünk. Egy shapefile réteg ugyanis mindig több fájlból áll, melyekből minimálisan az *.SHP, *.SHX és *.DBF kiterjesztésűek kellenek ahhoz, hogy egy térinformatikai rendszer értelmezni tudja őket. Hozzunk létre egy munkakönyvtárat a számítógépünkön és másoljuk bele ebbe az említett fájlokat!

Ezt követően az sf csomag segítségével beolvassuk a shapefile tartalmát az R-be. A továbbiakban felhasználom még a listr és az igraph csomagokat is, ezért egyúttal ezeket is betöltöm a rendszerbe. (Egy adott függvény származási helyének egyértelmű beazonosításához a csomag::függvény() formulát használom a kódban. Az automatikusan betöltődő függvényeket nem jelölöm külön.)

Az alábbiakban R nyelven (v4.2.2) írt kódot használok a feladat végrehajtásához. A magyarázó szövegek közé ékelt fekete kódblokkok tartalmát az RStudio-ban egymás alá illesztve elvileg bárki által reprodukálható az itt bemutatott műveletsor. A kódblokkok # kezdetű sorai pusztán magyarázó funkcióval bírnak, ezekre a program futtatásakor nincs szükség.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# A szükséges csomagok betöltése.
# Ezeket a legelső használat előtt az install.packages("csomag_neve") utasítással telepíteni kell.
library(sf)
library(listr)
library(igraph)

# Ide a saját munkakönyvtárunk elérését kell beírni!
setwd("C:/Munkakönyvtár")

# A shapefile réteg beolvasása
telepulesek <- sf::st_read("MO_Telepules.shp")
sf::st_crs(telepulesek) <- 3857

A kód futtatását követően az MO_Telepules shapefájl tartalma bekerült a telepulesek nevű adatkeretbe. Vessünk egy pillantást ennek első 10 sorára!

1

print(telepulesek)

## Simple feature collection with 12622 features and 9 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 1601630 ymin: 5539229 xmax: 2950561 ymax: 6381733
## Projected CRS: WGS 84 / Pseudo-Mercator
## First 10 features:
##            telepulesn status JaraszSzek MegyeSzekh  IDJaras IDJGeo IDMegye
## 1                Pécs   tjv.          T          T M0100001  M0104     M01
## 2          Albertfalu   <NA>          F          F    M0101  M0101     M01
## 3          Baranyabán   <NA>          F          F    M0101  M0101     M01
## 4  Baranyaszentistván   <NA>          F          F    M0101  M0101     M01
## 5          Baranyavár   <NA>          F          F    M0101  M0101     M01
## 6              Bellye   <NA>          F          F    M0101  M0101     M01
## 7               Benge   <NA>          F          F    M0101  M0101     M01
## 8             Bezedek   <NA>          F          F    M0101  M0101     M01
## 9             Bolmány   <NA>          F          F    M0101  M0101     M01
## 10              Csúza   <NA>          F          F    M0101  M0101     M01
##    Nev_megj IDTel1910                       geometry
## 1      <NA>  M0100001 MULTIPOLYGON (((2031533 580...
## 2      <NA>  M0101001 MULTIPOLYGON (((2087952 573...
## 3      <NA>  M0101002 MULTIPOLYGON (((2075914 575...
## 4      <NA>  M0101003 MULTIPOLYGON (((2064136 574...
## 5      <NA>  M0101004 MULTIPOLYGON (((2076479 574...
## 6      <NA>  M0101005 MULTIPOLYGON (((2084531 571...
## 7      <NA>  M0101006 MULTIPOLYGON (((2068382 574...
## 8      <NA>  M0101007 MULTIPOLYGON (((2073360 575...
## 9      <NA>  M0101008 MULTIPOLYGON (((2063308 573...
## 10     <NA>  M0101009 MULTIPOLYGON (((2088661 575...

A fentiekből kiolvasható, hogy itt egy georeferált térképpel van dolgunk. Ez a lényegét tekintve azt jelenti, hogy a digitális térkép minden pontja a valós földrajzi tér egy konkrét pontjának felel meg. A georeferálást a WGS 84 / Pseudo-Mercator (EPSG:3857) rendszerben készítették el a térkép alkotói. A töréspontok koordinátáihoz ennek megfelelő formátumban juthatunk majd hozzá.

A táblázat minden sora egy önálló térképi objektum, jelen esetben egy-egy település adatait tartalmazza. Ezeknek megkapjuk például a nevét (telepulesn), az egyedi azonosítóját (IDTel1910), valamint az őket szimbolizáló poligonok határvonalán lévő töréspontok koordinátáit (geometry).

A konzolra kinyomtatott információk szerint összesen 12.622 térképi objektumunk van. Mivel a korábban elmondottak alapján már tudjuk, hogy 12.541 település volt a kérdéses területen, ezért felmerülhet a gyanú, hogy valami hibával találkoztunk. Az adatkeretünknek tehát több sora van, mint ahány településnek lennie kell.

Némi szemrevételezéssel elég könnyen kiszúrható, hogy mi itt a gond: a táblázatban lévő térképi objektumok egy része nem kapcsolódik egyetlen településhez sem. Míg az utóbbiak egyedi azonosítója rendre M betűvel kezdődik és további 7 számot tartalmaz (pl. M2507000), addig van 80 olyan sor, ahol az azonosítók egy X betűből és 2 számjegyből állnak (pl. X12). Ezekhez név sincs megadva. Én azt feltételezem, hogy olyan exklávékkal van dolgunk, amelyek hovatartozását a térkép készítőinek nem sikerült megállapítani.

A problémára jó megoldás nincs. A települések határai közé beékelődő apró, sehova nem tartozó foltok bizonyos esetekben megzavarhatják a szomszédsági kapcsolatok detektálását. Ugyanakkor ha megtartjuk őket, akkor ezen túlmenően még azzal is számolnunk kell, hogy kvázi településként jelennek meg a rendszerben. Számomra az tűnik a kisebbik rossznak, ha ezeket a – három karakterből álló azonosítóval rendelkező – sorokat eltávolítjuk a táblázatból. Egyúttal hasonlóan járok el a szűkebb Magyarországtól földrajzilag különálló Fiume esetében is. Megjegyzem továbbá, hogy valami miatt a térképről hiányzik három település neve (Kistapolca, Hernádfa, Szepeshely). Ezeket a rend kedvéért pótolom.

1
2
3
4
5
6
7

# A három karakterből álló azonosítóval rendelkező sorok, valamint Fiume eltávolítása.
telepulesek <- telepulesek[- c(which(nchar(telepulesek$IDTel1910) == 3), which(telepulesek$telepulesn == "Fiume")),]

# A hiányzó településnevek pótolása.
telepulesek[telepulesek$IDTel1910 == "M0106028", "telepulesn"] <- "Kistapolca"
telepulesek[telepulesek$IDTel1910 == "M0107024", "telepulesn"] <- "Hernádfa"
telepulesek[telepulesek$IDTel1910 == "M3308022", "telepulesn"] <- "Szepeshely"

A továbbiakban a célunk az, hogy a shapefile réteget reprezentáló adatkeretből kinyerjük a hálózat létrehozásához szükséges információkat, majd feldolgozzuk azokat.

Az összegyűjtött adatokat a telepulesek2 nevű listában fogom tárolni. Ennek annyi eleme lesz, ahány településünk van.

Első körben létrehozom magát a listát, majd elhelyezem benne a települések nevét, egyedi azonosítóját, valamint a települést szimbolizáló poligon töréspontjainak koordinátáit. Itt az érdemi feladatot ez utóbbiak előállítása jelenti a telepulesek adatkeret geometry oszlopának tartalmából. A művelet során tekintettel kell lenni arra is, hogy az exklávék vagy enklávék esetleges jelenléte megbonyolítja azt az adatstruktúrát, amelyből ki akarjuk nyerni a koordinátákat. Nézzünk egy-egy vizuális példát ezekre!

Az exklávék halmazából álló Majsamiklósvár (felső kép) és az enklávét tartalmazó Somogyszil (alsó kép) ugyan speciális és ritka eseteket jelenítenek meg, de az algoritmusnak természetesen tudni kell kezelni az ilyeneket is.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

# Az üres lista létrehozása
telepulesek2 <- vector(mode = "list", length = nrow(telepulesek))

# Egyenként végigmegyünk a lista elemein és feltöltjük azokat az alapadatokkal.
for (i in 1:length(telepulesek2)) {
  # A haladást jelző számláló.
  print(i)
  # A település ID-jének rögzítése.
  telepulesek2[[i]][[1]] <- telepulesek$IDTel1910[i]
  # A település nevének rögzítése.
  telepulesek2[[i]][[2]] <- telepulesek$telepulesn[i]
  # A töréspontok koordinátáinak kinyerése és rögzítése.
  x <- do.call(rbind.data.frame, listr::list_flatten(telepulesek$geometry[[i]]))
  telepulesek2[[i]][[3]] <- paste(as.character.numeric_version(x$V1), as.character.numeric_version(x$V2))
  # A következő körben ide kerülnek majd a szomszédos települések ID-jei.
  telepulesek2[[i]][[4]] <- character()
  # A következő körben ide kerül majd a település Ego-hálózata.
  telepulesek2[[i]][[5]] <- data.frame()
}

Ha idáig eljutottunk, akkor már minden település esetében rendelkezésünkre állnak azok a koordináták, amelyek a határvonalukon lévő töréspontokhoz tartoznak. Nézzünk erre egy példát!

1

telepulesek2[[998]][[3]]

##  [1] "2029600 5907222" "2029674 5908292" "2030042 5909085" "2030413 5909777"
##  [5] "2031369 5910516" "2031983 5910799" "2032611 5910735" "2033097 5910551"
##  [9] "2033028 5910466" "2033016 5910445" "2032955 5910298" "2032825 5909992"
## [13] "2032630 5909526" "2032379 5908928" "2032173 5908435" "2032034 5908102"
## [17] "2031961 5907928" "2031950 5907901" "2031916 5907891" "2031622 5907807"
## [21] "2030989 5907624" "2030754 5907556" "2030024 5907345" "2029600 5907222"

Ezek a koordináták az összes településre lekérdezhetők. Amennyiben ugyanaz a töréspont-koordináta több helységnél is előfordul, akkor az érintett települések szomszédosnak tekinthetők egymással. A következő lépésben minden települést összehasonlítok minden másik településsel, és ennek során megnézem, hogy a töréspont-koordinátáik halmazában van-e közös rész.

Az alábbi kód nagyon lassan fut le, hiszen 12.541*12.540 = 157.264.140 összehasonlítást kell elvégeznie a számítógépnek. Az algoritmuson bizonyára lehetne finomítani a sebesség növelése érdekében. Mivel azonban csak egyetlenegy alkalommal kell lefuttatni, ezért nem foglalkozom ezzel a kérdéssel.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# Egyenként végigmegyünk a lista elemein, azaz a településeken.
for (i in 1:length(telepulesek2)) {
  # A haladást jelző számláló.
  print(i)
  # Az aktuálisan vizsgált települést összehasonlítjuk az összes többivel.
  for (j in setdiff(1:length(telepulesek2), i)) {
    # Ha van közös része a töréspontok koordinátáinak, akkor feljegyezzük a szomszédosnak minősített település azonosítóját.
    if (length(intersect(telepulesek2[[i]][[3]], telepulesek2[[j]][[3]])) != 0) {
      telepulesek2[[i]][[4]] <- c(telepulesek2[[i]][[4]], telepulesek2[[j]][[1]])
    }
  }
  # A település Ego-hálózatát alkotó kapcsolatok legenerálása
  telepulesek2[[i]][[5]] <- expand.grid(telepulesek2[[i]][[1]], telepulesek2[[i]][[4]])
}

A kód végrehajtása során mindegyik településhez feljegyzésre kerülnek a vele szomszédos helységek ID-jei, valamint a település Ego-hálózatát alkotó kapcsolatpárok is legenerálásra kerülnek egy adatkeretben. Az alábbiakban Eger (M2507000) példáját láthatjuk.

1

telepulesek2[[5202]][[4]]

## [1] "M2501003" "M2501005" "M2501008" "M2501009" "M2505007" "M3003015" "M3003020"
## [8] "M3003022" "M3003026"

1

telepulesek2[[5202]][[5]]

##       Var1     Var2
## 1 M2507000 M2501003
## 2 M2507000 M2501005
## 3 M2507000 M2501008
## 4 M2507000 M2501009
## 5 M2507000 M2505007
## 6 M2507000 M3003015
## 7 M2507000 M3003020
## 8 M2507000 M3003022
## 9 M2507000 M3003026

Ezt követően a telepulesek2 listából kibányásszuk és edges néven önálló adatkeretben egyesítjük az Ego-hálózatokat alkotó kapcsolatpárokat. E művelet során a duplán (oda-vissza) rögzített szomszédsági kapcsolatoktól is megszabadulunk.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

# Az üres adatkeret létrehozása.
edges <- data.frame("Var1" = as.character(), "Var2" = as.character())
# Az Ego-hálózatot alkotó kapcsolatok összegyűjtése és egyesítése.
for (i in 1:length(telepulesek2)) {
  print(i)
  edges <- rbind(edges, telepulesek2[[i]][[5]])
}
edges$Var1 <- as.character(edges$Var1)
edges$Var2 <- as.character(edges$Var2)

# A duplumok eltávolítása az adatkeretből.
edges$source <- ifelse(edges$Var1 < edges$Var2, edges$Var1, edges$Var2)
edges$target <- ifelse(edges$Var1 > edges$Var2, edges$Var1, edges$Var2)
edges <- unique.data.frame(edges[, 3:4])

Az utolsó lépés a hálózatos adatstruktúra létrehozása a fentiekben előállított adatokból. Ha ezt elmentjük, akkor a továbbiakban elég a hálózatot betölteni és nincs szükség az itt közölt algoritmus újbóli lefuttatására.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

# A hálózat létrehozása
halozat <- igraph::graph_from_data_frame(edges, directed = F)
V(halozat)$label <- telepulesek[match(V(halozat)$name, telepulesek$IDTel1910),]$telepulesn

# A hálózat elmentése a későbbi felhasználás céljából.
save(halozat, file = "halozat.RData")

# A hálózat visszatöltésének lehetősége a munkakönyvtárból.
# load("halozat.RData")

# A hálózat jellemzőinek kiíratása a konzolra.
print(halozat)

## IGRAPH 26ca71b UN-- 12541 37198 -- 
## + attr: name (v/c), label (v/c)
## + edges from 26ca71b (vertex names):
##  [1] M0100001--M0102005 M0100001--M0102043 M0100001--M0102053 M0100001--M0104014
##  [5] M0100001--M0104019 M0100001--M0104026 M0100001--M0104027 M0100001--M0104029
##  [9] M0100001--M0104030 M0100001--M0104032 M0100001--M0104037 M0100001--M0104047
## [13] M0101001--M0101012 M0101001--M0101018 M0101001--M0101026 M0101002--M0101004
## [17] M0101002--M0101007 M0101002--M0101014 M0101002--M0101017 M0101002--M0101023
## [21] M0101002--M0101032 M0101002--M0101033 M0101002--M0103003 M0101002--M0103011
## [25] M0101002--M0103016 M0101003--M0101006 M0101003--M0101008 M0101003--M0101016
## [29] M0101003--M0101028 M0101003--M0101031 M0101003--M0106004 M0101003--M0106040
## + ... omitted several edges

A hálózat tehát egy csomópontokból és kapcsolatokból álló adatstruktúra, ami alapesetben nem egy látványos valami. Megfelelően elrendezett adatok halmaza. Jelen esetben 12.541 csomópontból és köztük 37.198 kapcsolatból áll. Én a későbbiekben az igraph csomagban található függvények segítségével fogom kezelni a települések hálózatát. De nem maga a hálózat a fontos a számomra, mert ez csak egy eszköz. Ennek segítségével választom majd ki a kutatásom szempontjából releváns településeket.

Megtehetjük ugyanakkor azt is, hogy a hálózatunkat átkonvertáljuk a Gephi nevű hálózatkezelő szoftver által megkívánt formátumba, s így egy vizuális reprezentációt kölcsönzünk neki. Az alábbi ábrán a GeoLayout plugint használtam fel a csomópontok elrendezéséhez. Ennek köszönhetően a hálózatba kötött települések a földrajzi elhelyezkedésüknek megfelelő helyen látszódnak. A helységek WGS 84 formátumú koordinátáit egy régebbi posztban leírt módon állítottam elő egy másik térkép átkonvertálásával.

Íme a láthatóvá varázsolt végeredmény!