UAV és GIS találkozása: Hova érdemes kilátót vagy rádiótornyot építeni?

Előszeretettel szokás bemutatni a GIS-ben rejlő adatelemző potenciált az optimális térrészek leválogatásával. Jelenlegi cikkünk valójában erről fog szólni, tehát megkeressük, hogy egy adott területen belül, meghatározott kritériumok szerint, mely pontok felelnek meg egy kilátó vagy rádiótorony építésére.

Hogy jön ide a drón? A GIS szoftverek adatok nélkül féllábú óriások! A GIS által használt alapadatok – mint például egy DEM, településhatár vagy úthálózat térkép – légi felvételezéssel is előállítható. Ha érdekel, hogy egy vektoros úthálózat vagy településhatár térképet hogyan állíthatsz elő ortofotóból, akkor olvasd el a korábbi cikkünket.

A kilátó/rádiótorony optimális helyének kiválasztásához a következő kritériumokat határoztuk meg:

1. A környék legmagasabb pontja és az attól 150 méterrel alacsonyabb térszínek jöhetnek szóba, hogy még elég magasan legyen a kilátó/ rádiótorony.
2. A településektől minimum 500 méterre található helyeket keressük, hogy ne zavarjuk a lakosságot.
3. Az úthálózattól maximum 100 méterre fekvő területek alkalmasak, hogy autóval könnyen megközelíthető legyen.
4. A kilátó/rádiótorony 15 méter magas lesz.
5. A rendelkezésre álló helyszínek közül a legnagyobb területet belátó pont (az aktuális terepfelszín felett 15 méter magasságban) legyen kiválasztva. Ez azért fontos, mert ha kilátóról beszélünk, akkor a látogatóknak minél nagyobb területet szeretnénk bemutatni, ha rádiótoronyról van szó, akkor a lehető legkevesebb árnyékoló, zavaró hatást kell elérnünk.

Megjegyzés: A műveletsort a legújabb QGIS Desktop 3.24.0 szoftverben végeztük el. Ha régebbi QGIS-t használsz, akkor felhívom a figyelmedet, hogy az itt leírtak csak a QGIS 3.x verziókkal fog működni.
Továbbá, ha külön indítható a GRASS GIS függvényekkel ellátott változat (pl. QGIS Desktop 3.14.0 with GRASS 7.8.3), akkor a GRASS funkciókkal bővített változatot indítsd el. Az új QGIS szoftverbe már alapból benne vannak a GRASS függvények.

LEGMAGASABB PONTOK LEVÁLOGATÁSA GIS-BEN

Annak érdekében, hogy a vizsgált terület legmagasabb pontját és az attól 150 méterrel alacsonyabb részeket ki tudjuk jelölni, szükségünk lesz egy DEM-re. Ha megnyitjuk a domborzatmodellt QGIS alatt, akkor már a rétegfában láthatóvá válik, hogy hány méter a legmagasabb pontunk. 

A „Processing Toolbox” alatt kikereshető a „Reclassify by Table” függvény. Ezt megnyitva a DEM réteget kell kiválasztani, majd ezen belül a „Reclassification table” opció alatt állíthatod be a magassági értékek intervallumát. A példámban a „Value” értéket 1-re állítottam be. Ez csak annyit jelent, hogy a raszter értékét 1-re állítja a keresett cellákban, minden más „NoData” értéket kap.

Javaslom, hogy a „Use no data when no range matches value” opciót kapcsold be. A következő képen láthatod, hogy én milyen beállításokon futtattam le a függvényt.

Ha minden jól ment, akkor valami fekete paca fog megjelenni a térképen. Ez egy raszteres állomány és a továbbiakban vektorokkal fogunk dolgozni. Van olyan eset, amikor könnyebb, ha egységes adatmodellt használunk és véleményem szerint ez is ilyen. Konvertáljuk át a kijelölt területet vektorrá! Írd be a „Processing Toolbox” keresőjébe, hogy „Polygonize (raster to vector)”. Az előbukkanó GDAL függvényt nyisd meg, majd válaszd ki az előbb előállított rasztert az „Input layer” mezőnél. 

Az eredményed egy vektoros állomány lesz, ami jól illeszkedik a raszter határaihoz.

TELEPÜLÉSEK ÖVEZETÉNEK SZŰRÉSE

A kritérium szerint elő kell állítani a települések 500 méteres övezetét, majd ezeket a területeket ki kell törölnünk a potenciális térrészekből. Ha így teszünk, akkor a fennmaradó helyek (légvonalban) távolabb fognak esni a lakott zónáktól, mint 500 méter. Először is nyissuk meg a települések vektoros réteget. Ezután, minden objektum köré egy övezet generálható a felső menüszalagban található „Vector/Geoprocessing Tools/Buffer…” függvénnyel.

Válaszd ki a települések vektoros réteget, állíts be 500 méter távolságot.

Olyan az eredmény, mintha elhízott települések lennének szétfolyva a térképen. ☺ Látható, hogy több olyan terület is van, ami átfedésben áll az előző fejezetben leválogatott dombsági helyekkel.

A metszetet képező területeket kell törölnünk, hiszen azok nem felelnek meg a kritériumrendszerünknek. A legegyszerűbben úgy tudod kivágni az alkalmatlan területeket, ha a „Vector/Geoprocessing Tools/Difference…” függvényt használod.

Itt az „Input layer” helyén add meg a vektoros magasságilag optimális helyeket, majd az „Overlay layer” esetében az övezettel ellátott településeket.

Látható, hogy az új részeredménnyel már csökkent az optimális terület mérete.

ÚTHÁLÓZATHOZ KÖZELI TERÜLET KIJELÖLÉSE GIS-BEN

Ebben a fejezetben a harmadik kritérium szerint válogatjuk le a területeket. Ennek értelmében az úthálózattól 100 méterre található helyeket kell kijelölni. Az előzőkben bemutatott „Buffer” függvény segítségével készítsük el az úthálózat réteg 100 méteres övezetét.

Az eredmény olyan, mint egy giliszta térkép. Vajon létezik giliszta fajtákat bemutató térkép…? 🙂 Szemben az előző fejezetben bemutatott, törléssel járó művelettel, itt a közös pontokra vagyunk kíváncsiak, ezért használjuk a „Vector/Geoprocessing Tools/Clip…” függvényt. 

Legegyszerűbb úgy lehet elképzelni a „Clip” működését, hogy van egy tésztánk (az előző fejezetben leválogatott terület) és van egy tésztaformánk (úthálózat övezete). A tésztaformával ki akarjuk szaggatni a releváns területeket, hiszen ez lesz értékes számunkra. Az „Input layer” amiből ki akarunk vágni, tehát az előző fejezet végére előállt „Difference”, míg az „Overlay layer” a forma, amit ki akarunk vágni, azaz az úthálózat övezete.

Az eredményen látszódik, hogy tovább szűkítettük a megfelelő helyeket.

A BELÁTHATÓ TERÜLETEK MÉRTÉKE

A belátható területek elemzésére („Viewshed”) kettő függvény használható QGIS-ben: GRASS és GDAL. Mindkettő „Viewshed” függvény csak pontokkal képes kiszámolni a belátható területeket, de nekünk poligon áll rendelkezésünkre. Ez a fejezet egy kicsit bonyolultabb lesz, mint az előzők, mivel több részfeladat megoldásából jutunk el a célunkig:

1. Első lépésben a poligon területére pontokat kell illeszteni.
2. Második lépésként az összes szóba jöhető pontra le kell futtatni a „Viewshed” függvényt, majd ki kell számolni mindegyik térkép belátható területeit hektárban.
3. Végezetül, a harmadik lépés, hogy kiválasszuk azt a pontot, amelyről a legnagyobb terület látható be.

Nem mondom… szép kihívás! Fokozatokban közelítsük meg a problémát!

AZ OPTIMÁLIS TERÜLETEN BELÜLI PONTOK KIOSZTÁSA

A pontok illesztését megteheted, ha a „Vector/Research Tools/Random Points Inside Polygons…” funkciót kiválasztod.

A függvény segítségével a leválogatott területünkön belül szabályos mintázatban elhelyezhetünk pontokat. Válasszuk ki az „Input layer” esetében a kijelölt területet (clipped), majd állítsd be, hogy mennyi pontot szeretnél elhelyezni a területen belül „Point count or density”. Itt érdemes nem irreálisan magas számot megadni, hiszen ne feledd, hogy mindegyik pontra különálló térképeket is létre kell hozni! Én beállítottam a minimális távolságot is a pontok között („Minimum distance between points”), hogy egyenletesebb legyen a térbeli kiosztás.

Az eredmény egy jól kivehető, szabályos mintázatú ponthalmaz lett. Igaz, a távolság beállítása miatt az előre kért 50 pont helyett 53 pontot generált a QGIS.

VIEWSHED TÉRKÉP GENERÁLÁSA

Talán a beláthatóság (viewshed) térkép gyártása jelenti a legnagyobb problémát. Gondolj csak bele, hogy 53 pont mindegyikére manuálisan kell lefuttatni a „Viewshed” függvényt, majd mindegyikre ki kell számolni a belátható területek mértékét. Ha ez több potenciális helyet vagy nagyobb mintavételi sűrűséget jelentene, akkor akár napok elmennének ezzel a lépéssel. De nem nekünk! A lusta ember – mint amilyen én is vagyok – megtanul folyamatokat automatizálni és hagyja, hogy a számítógép dolgozzon helyette. Ráadásul, a kóddal megismételhetővé, ellenőrizhetővé válik egy teljes folyamat. Írjunk egy Python kódot, ami elvégzi helyettünk a piszkos munkát! 🙂

Nyisd meg QGIS alatt ikonnal ellátott „Python Console”-t, majd a „Show Editor” gombot nyomd meg. A szöveges területbe másold be a következő kódot.

Megjegyzés: Mielőtt használnád a szkriptet, pár fontos információt megosztok veled. A kódban kommentelve jeleztem (a # jel utáni szövegrész), hogy milyen változókat kell átírnod. Tehát a kódsorban felülről lefelé haladva: DEM réteg neve; a pontokat tartalmazó réteg neve; a „viewshed” térképek kimentésének abszolút elérési útja a fájl nevével és „.tif” végződéssel. Ügyelj arra, hogy az elérési utat normál perjellel „/” add meg és ne fordított perjellel „\”.
A másik fontos dolog, hogy GRASS függvénnyel írtam meg a kódot. Ha erről többet szeretnél megtudni, akkor itt elolvashatod a részletes leírást. Az ‘observer_elevation’ értékét 15 méterre állítottam be a példa kedvéért. Ha te más paraméterrel rendelkező kilátót/rádiótornyot terveznél, akkor írd át ezt az értéket.

import processing
import numpy as np
from osgeo import gdal
 
# Enter the name of the DEM layer
rasterLayer = QgsProject.instance().mapLayersByName('dem')[0]
 
# Enter the name of the Points layer
pointlayer = QgsProject.instance().mapLayersByName('Random points')[0]
 
# Enter the file path of the output raster
outputFile = 'C:/Users/.../YOUR FOLDER/viewshed.tif'
 
new_column = pointlayer.dataProvider()
new_column.addAttributes([QgsField('area', QVariant.Double)])
pointlayer.updateFields()
 
iter = pointlayer.getFeatures()
 
for feature in iter:
 
    geom = feature.geometry()
    observerxy = geom.asPoint()
    strobserverxy = str(observerxy[0]) + "," + str(observerxy[1])
 
    processing.run('grass7:r.viewshed', {
    'input': rasterLayer,
    'output': outputFile,
    'coordinates': strobserverxy,
    'observer_elevation': 15,
    '--overwrite':True,
    '-b':True
    })
 
    viewshedLayer = QgsRasterLayer(outputFile)
    
    pixel_size_x = viewshedLayer.rasterUnitsPerPixelX()
    pixel_size_y = viewshedLayer.rasterUnitsPerPixelY()
 
    hectare_per_cell = pixel_size_x * pixel_size_y / 10000
 
    ds = gdal.Open(outputFile)
 
    raster = ds.ReadAsArray()
    
    all_pixels_1 = np.where((raster == 0)|(raster == 255), np.nan, raster)
    
    total_pixel_area = np.nansum(all_pixels_1)
    total_area_raster = total_pixel_area * hectare_per_cell
    atts = {pointlayer.fields().lookupField('area'):float(total_area_raster)}
    new_column.changeAttributeValues({feature.id(): atts})

Miután beállítottál mindent, nyomj rá a zöld háromszög gombra („Run Script”). Ezután nyugodtan menj és igyál meg egy kávét, megérdemled! 🙂 Ugye, hogy mennyivel jobb, hogy a gép dolgozik helyetted?

Az előzőkben bemutatott kód úgy fut le, hogy az összes pontra kiszámolja a beláthatósági térképet, de mindig felülírja a kimeneti fájlt. Ezt helytakarékossági szempont miatt írtam meg így. Ha szükséged lenne egy olyan megoldásra, ami az összes térképet külön kimenti, akkor használd a következő kódsort.

Megjegyzés: Az előző szkripthez írt megjegyzéshez képest csupán egy apró részletet kell másképp csinálnod. A „viewshed” térképek kimentésének abszolút elérési útja végén, a fájl neve után ne írd be, hogy „.tif”!

import processing
import numpy as np
from osgeo import gdal
 
i = 0
 
# Enter the name of the DEM layer
rasterLayer = QgsProject.instance().mapLayersByName('dem')[0]
# Enter the name of the Points layer
pointlayer = QgsProject.instance().mapLayersByName('Random points')[0]
# Enter the file path of the output raster
outputFile = 'C:/Users/.../YOUR FOLDER/viewshed'
 
new_column = pointlayer.dataProvider()
new_column.addAttributes([QgsField('area', QVariant.Double)])
pointlayer.updateFields()
 
iter = pointlayer.getFeatures()
 
for feature in iter:
 
    outputFileTif = outputFile + str(i) + '.tif'
    
    geom = feature.geometry()
    observerxy = geom.asPoint()
    strobserverxy = str(observerxy[0]) + "," + str(observerxy[1])
 
    processing.run('grass7:r.viewshed', {
    'input': rasterLayer,
    'output': outputFileTif,
    'coordinates': strobserverxy,
    'observer_elevation': 15,
    '--overwrite':True,
    '-b':True
    })
 
    viewshedLayer = QgsRasterLayer(outputFileTif)
    
    pixel_size_x = viewshedLayer.rasterUnitsPerPixelX()
    pixel_size_y = viewshedLayer.rasterUnitsPerPixelY()
 
    hectare_per_cell = pixel_size_x * pixel_size_y / 10000
 
    ds = gdal.Open(outputFileTif)
 
    raster = ds.ReadAsArray()
    
    all_pixels_1 = np.where((raster == 0)|(raster == 255), np.nan, raster)
    
    total_pixel_area = np.nansum(all_pixels_1)
    total_area_raster = total_pixel_area * hectare_per_cell
    atts = {pointlayer.fields().lookupField('area'):float(total_area_raster)}
    new_column.changeAttributeValues({feature.id(): atts})
    i += 1

Néhány adat a szkript futásáról: A cikk írása során én egy 5×5 méteres felbontású DEM-mel dolgoztam, ami 3028×1899 raszterből áll. 53 pontra vizsgáltam a beláthatóságot és az abból kiszámolható terület méretét. A kódot Windows 10 Pro 64 bites operációs rendszerre telepített QGIS Desktop 3.24.0 futattam. A PC specifikációi: Intel Core i9-10900F CPU 2.81 GHz; 16 GB RAM; Radeon RX 570 GPU 4GB GDDR5. A kód futási ideje 30 percet vett igénybe.

AZ OPTIMÁLIS TELEPÍTÉSI HELY KIVÁLASZTÁSA

A szkript futtatása után az eredmény könnyen leolvasható a térképről, ha megnyitjuk a pontok attribútum tábláját. Ezt úgy teheted meg, ha jobb egérgombbal kattintasz a rétegre, majd kiválasztod az „Open Attribute Table” opciót.

Kattints kétszer az új „area” nevű oszlop fejlécére, így az értékeket csökkenő sorrendbe állíthatod. Így már csak a sor elején lévő számra kell kattintani, majd kiválaszthatjuk a legmagasabb értéket. Ez azt jelenti, hogy ebből a kiválasztott pontból látható be a legnagyobb terület. A kód megírása során a területegységeket hektárban határoztam meg.

ZÁRÓ GONDOLATOK

Ezzel a cikkel az volt a szándékom, hogy népszerűsítsem a folyamatok automatizálását open-source környezetben. Remélem sokaknak hasznos lesz ez a tartalom is.

Az itt bemutatott kódokat bármilyen projekt keretében nyugodtan felhasználhatod! 🙂

Ha nagyon-nagyon tetszett, amit olvastál, akkor meg is oszthatod az ismerőseiddel. Ne fogd vissza magad! 🙂

Tetszett, amit olvastál? Akarsz hasonlókat olvasni?