Atgal...

QGIS: Izolinijų kūrimas

Įvadas

Ar kuriant topografinį, ar orientacinį, o kartais ir tiesiog turistinį žemėlapį, mums reikia vietovės aukščio duomenų. Tam paprastai naudojamos izolinijos (izohipses, horizontalės). Tokie aukščio duomenys leidžia suprasti, koks yra vietovės reljefas, geriau suplanuoti žygį ar tiesiog tiksliau susiorientuoti vietovėje (pvz. miške, kur kitų orientyrų be aukščio gali ir nebūti).

Izolinijos

Lietuvoje izolinijas galime parsisiųsti iš TOP50LKS duomenų rinkinio. Jame, be visų kitų, rasite ir izohipsių (bei bergštrichų) sluoksnius:

Izolinijos TOP50LKS

Tai puikus duomenų rinkinys. Deja, jis paskaičiuotas būtent masteliui 1:50000. Taigi jei norite aukščio informacijos stambesniam masteliui, šito duomenų rinkinio neužteks. Ką tada daryti? Ogi pasidaryti savo izohipses naudojant QGIS ir kitus atviro kodo įrankius!

Pradinių duomenų atsisiuntimas

Tam, kad pasigamintumėme savo izohipses, mums reikia aukščio duomenų. Lietuvos teritorijai šiuo metu tinkamiausias duomenų rinkinys - „SEŽP_0.5LT – LR teritorijos skaitmeniniai erdviniai žemės paviršiaus lazerinio skenavimo taškų duomenys“. Tai 2009 metais sukurtas rinkinys. Šiuo metu vyksta naujas Lietuvos paviršiaus lazerinis skenavimas, bet rezultatų palaukti teks bent iki 2023 metų.

Šiuos duomenis galima atsisiųsti iš Lietuvos erdvinių duomenų portalo - geoportalo.

SEŽP duomenys suskirstyti į atskirus lapus:

SEŽP pasirinkimas. 1 žingsnis

Naudodami pasirinkimo žemėlapyje įrankius, nurodykite, kurios konkrečiai vietos aukščio duomenų pageidaujate:

SEŽP pasirinkimas. 2 žingsnis

Spauskite mygtuką „Išsaugoti pasirinkimą“. Patvirtinkite užsakymo detales (nesibaiminkite, kaina bus 0), patvirtinkite asmenybę, patvirtinkite užsakymo detales ir šiek tiek palaukite, kol jums bus paruoštas užsakytas duomenų rinkinys (paprastai tai užtrunka tik 5-10 minučių).

Pastaba Detalesnes instrukcijas, kaip atsisiųsti duomenų rinkinius iš Geoportalo rasite šiame pradžiamokslyje.

Duomenis gausite ASCII Simple Point Cloud (*.xyz) formatu. T.y. jūsų gautas failas yra paprastas tekstinis failas, kurio kiekvienoje eilutėje yra vieno taško informacija - platuma, ilguma ir aukštis:

 590002.840  6100563.830 155.258
 590000.610  6100525.350 156.018
 590000.320  6100497.610 156.968
 590004.380  6100542.570 155.628
 590005.440  6100554.230 155.258
 590004.300  6100540.910 155.778
 590001.820  6100499.800 156.878
 ...

Neprivaloma pradinių duomenų vizualizacija

Mums realiai reikia šituos duomenis konvertuoti į aukščio gardelę (rastrą), bet jei labai įdomu pažiūrėti, kas ten per duomenis, jei norite geriau suprasti, ką mes darysime apačioj - galite atsidaryti šį failą QGIS'e naudodami „Atskirto teksto atidarymo“ funkcionalumą (dar žinomą kaip CSV importą).

Pakartosiu: šito daryti nėra būtina, izohipsėm šito žingsnio nereikia. Bet jei vis tiek norite pabandyti, tai atskirtų duomenų importe nurodykite, kad:

XYZ duomenų importas

Viskas, spauskite „Pridėti“ ir po kokios minutės gausite didžiulę aukščio taškų aibę. Šiek tiek patvarkę vizualizavimą, gauname tokį vaizdą:

Aukščio taškai

Ką mes čia matome? Darant lazerinį skenavimą, pagal požymius, kurie gerokai išeina už šio straipsnio temos, nustatoma, ar lazeris pataikė į žemę, medį, vandenį ar pastatą. SEŽP faile yra tik tie taškai, kurie pataikė į žemę. Taigi paveikslėlyje matote, kad (beveik) nėra taškų vandenyje, ant pastatų, o štai miške taškų yra gerokai mažiau - nes tik dalis lazerio spindulių sugebėjo pasiekti žemę - didžioji dauguma atsitrenkė į lapus ir grįžo.

Ir jei vanduo mus nelabai duomina, nes izohipsių ten mums ir nereikia, pastatai irgi neesminė detalė, tai faktas, kad miške yra gerokai mažiau taškų mums yra labai svarbus, nes tai įtakos mūsų tolimesnius sprendimus konfigūruojant aukščio rastro kūrimą.

Patarimas nenurimstantiems: simbolizuoti šituos aukščio duomenis galite graduotai, kaip reikšmę nurodykite „$z“ - aukštį - vizualizuosite ne tik taškų pasisirstymą, bet ir aukštį.

Aukščio rastro kūrimas

Tam, kad sukurtumėm kontūrus, mums reikia turėti aukščio rastrą. T.y. tokios paties dydžio elementų/celių matricą su aukščio reikšme (paprastai tokia informacija įrašoma į GeoTIFF failus). Tam mes naudosime apdorojimo įrankinėje randamus Rastro analizės -> Tinklelio įrankius arba tiesiog per meniu Rastras -> Analize -> Tinklelis ***.

Pradžioje panaudokime patį paprasčiausią - „artimiausio kaimyno“ (angl. nearest neighbour) algoritmą. Kiekvienai celei šis algoritmas tiesiog suras patį artimiausią tašką ir iš jo paims aukštį. Gauname tokį vaizdelį (čia priartintas vaizdas, mastelis 1:500)

Aukščio rastras 1

Rezultatas lyg ir neblogai atrodo, bet pažiūrėkime į sluoksnio savybes:

Sluoksnio savybės

Kaip matome, iš viso taškų yra 278 horizontaliai ir 232 vertikaliai (šitų skaičių mums prireiks šiek tiek vėliau), o vieno taško dydis yra ~3x3m. Tai gan stambi gardelė, jei jums tokia tinka, galite iš karto eiti į kitą skyrių. O jei norite gaminti izohipses stambesniam masteliui, tai tikriausiai norėsite gerokai mažesnio taško dydžio.

Taško dydžio Tinklelio įrankiuose tiesiogiai nurodyti negalima. Bet tose pačiose sluoksnio (jųsų pradžioje sugeneruoto artimiausio kaimyno rastro) savybėse pažiūrėkite, koks yra rastro plotis ir aukštis. Kadangi plotis ir aukštis buvo 278 232 (juos mes įžiūrėjome sluoksnio savybėse aukščiau), tai jei norime vietoj 3x3m gardelės turėti keturi kartus mažesnę (pagal ilgį, ir 16 kartų pagal plotį), tiesiog padauginsime plotį bei aukštį iš 4 ir rezultatą nurodysime Tinklelio įrankių laukelyje „Papildomi komandinės eilutės parametrai [pasirinktinai]“ kaip parametrą „-outsize 1112 928“:

Algoritmo savybės 1

Vėl paleidžiame artimiausio kaimyno skaičiavimą ir jau gauname tokį vaizdą:

Aukščio rastras 2

Jau kaip ir geriau, bet matome, kad atsirado plėmai - artimiausio kaimyno algoritmas kiekvieno rastro taško paskaičiavimui naudoja tik vieną (artimiausią) kaimyną, todėl rezultatas niekaip neglotninamas, gauname „laiptelius“, kurie reikš prastenę gautų izohipsių kokybę. Mums reikia naudoti kitą metodą. Visų metodų neaprašinėsiu, apie juos galite pasiskaityti GDAL tinklelių dokumentacijoje. Pagrindinė mintis - mes kiekvienam gardelės elementui skaičiuoti norime naudoti ne vieną, o daugiau netoliese esančių taškų. Norime skaičiuoti svertinį vidurkį, t.y. arčiau esantys taškai turėtų daryti didesnę įtaką rezultatui, nei toliau esantys taškai. Tam mums tiks „atvirkštinio atstumo iki laipsnio“ tinklelio algoritmas (jį rasite tame pačiame meniu Rastras -> Analizė -> Tinklelis (atvirkštinis atstumas iki laipsnio). Paskaičiuokime tinklelį iš naujo, tik jau su atvirkštinio atstumo algoritmu ir didesniu outsize:

Aukščio rastras 3

Kaip matome, atsirado kažkokių šiukšlių. Problema ta, kad plotuose su medžiais stipriai trūksta duomenų, todėl reikia šiek tiek padidinti maksimalų atstumą, kuriuo šis algoritmas ieško gretimų taškų vidurkio skaičiavimui (nieko nenurodžius gretimų taškų ieškoma ~2-3m atstumu). Realiai jums reikėtų pasižiūrėti, kokio dydžio yra patys didžiausi tarpai ir tada nurodyti dvigubai mažesnį atstumą, kad taškas vidury taškų „vakuumo“ visgi rastų realių taškų skaičiavimams. Pabadykime nurodyti paieškos spindulį 10m (algoritmo parametruose tą patį skaičių įveskite abiems spinduliams: tiek pirmam - X ašiai, tiek antram - Y ašiai).

Algoritmo savybės 2

Aukščio rastras 4

Na štai, jau turime pakankamai detalų ir pakankamai glotnų aukštio rastrą, kuris tiks izohipsių kūrimui. Priminsiu, kad yra daugiau nei du aukštio rastro kūrimo algoritmai, kiekvienas iš jų turi nemažai parametrų. Taigi yra krūva būdų dar labiau pasitinkinti aukščio rastro ir tokiu būdu pagerinti galutinių kontūrų kokybę.

Izohipsių kūrimas

Kai jau turime aukščio rastro GeoTIFF failą, galima kurti kontūrus.

Kontūrai kuriami naudojant algoritmą „Kontūras“. Šiame algoritme reikia nurodyti mums reikiamą „Intervalą tarp kontūro linijų“. Tebūnie tai 1 (metras). Paskaičiuojam kontūrus:

Kontūrai 1

Dabar reikia sutvarkyti gautus kontūrus. Visų pirma mums reikia išimti mažytes detales, kurios prideda triukšmo, bet neduoda realios naudos. Tam galime naudoti algoritmą „Ištraukti pagal išraišką“. Palikime tik tuos kontūrus, kurių ilgis didesnis už 50 metrų - tam naudokime išraišką „$length > 50“:

Ištraukti pagal išraišką

Kontūrai 2

Toliau mums reikia panaikinti mažiukus vingius, tam naudosime paprastinimo algoritmą. Kadangi mano gardelė buvo kiek mažiau už metrą, tai ir paprastinsime mažesnius nei 1 kvadratinio metro ploto vingius naudodami Visvalingam paprastinimo algoritmą. Supaprastinti kontūrai atrodo taip:

Kontūrai 3

Na ir dabar beliko suglotninti rezultatą. Tam naudosime Glotninimo algoritmą su 5 iteracijomis:

Kontūrai 4

Štai ir viskas. Nepamirškite, kad izohipsių užrašai visada orientuojami į įkalnę, t.y. kai kuriose žemėlapio vietose aukščio užrašai matomi aukštyn kojomis ir taip ir turi būti. Žinoma galite rezultatą dar kartą paprastinti (didesne tolerancija) ir tada dar kartą glotninti. Taip pat galite ir paprastinimo bei glotninimo parametrus pabandyti keisti, kad gautumėte jums labiau tinkamą rezultatą. Jei tokių izohipsių jums reikia suskaičiuoti daug, tai visus šiuos veiksmus galite sudėti į QGIS Modelį.

2021-09-27