Přeskočit obsah

06. Automatizace v GRASS GIS, PyGRASS

Cvičení

Vstupní data

Tip pro Bash

Data DMR5G můžeme dekomprimovat dávkově:

for f in *.zip; do unzip $f;done

Lidarová data ve formátu LAS/LAZ lze importovat pomocí dvou nástrojů založených na knihovně PDAL:

  • r.in.pdal vytvoří rastrovou vrtsvu
  • v.in.pdal vytvoří vektorovou bodovou vrstvu

Poznámka

V systému GRASS najdeme nástroje r.in.lidar a v.in.lidar, které nabízí velmi podobnou funkcionalitu. Tyto nástroje jsou založené na zastaralé a v dnešní době již neudržované knihovně libLAS. Proto je nebudeme používat.

Tvorba rastrového DMR z lidarových dat

Vstupní data importujme pomocí r.in.pdal s nastaveným prostorovým rozlišením (resolution) výsledné rastrové vrstvy. Vstupní LAZ soubory neobsahují informaci o souřadnicovém systému, proto použijme přepínač -o. Důležitý je přepínač -e, který rozšíří výpočetní region na základě prostorového rozsahu vstupních lidarových dat.

r.in.pdal -oe input=155FGIS/06/dmr5g/KRAV59.laz output=KRAV59 resolution=5

Základní statistiku ve formě reportu vypišme pomocí nástroje r.report.

Důležité

Před výpočtem je nutné nastavit výpočetní region:

g.region raster=KRAV59

Rastrová vrstva obsahuje 7% buněk bez informace:

r.report map=KRAV59 nsteps=10
...
| *|no data. . . . . . . . . . . . . . . . . | 14231|  7.12|
...

Tyto chybějící informace doplňme pomocí nástroje r.fill.stats, který ve výchozím nastavení používá IDW:

r.fill.stats -k input=KRAV59 output=KRAV59_filled

Vyzkoušejme různé nastavení parametru distance (počet buňek pro interpolaci). Příklad pro distance=10:

Úkol

Počet buněk bez informace zjistěte pomocí r.univar nebo r.stats.

Úkol

Podobně stáhněte (stejné dlaždice) a naimportujte produkt DMP1G.

Výpočet model výšky porostu (CHM)

Canopy Height Model (CHM) je rozdíl výšek mezi digitálním model povrchem a terénu.

CHM vypočítejme pomocí nástroje mapové algrebry r.mapcalc:

r.mapcalc expression="KRAV59_chm = KRAV59_dmp - KRAV59"

Příklad CMH v tabulce barev differences:

Úkol

Stáhněte vrstvu budov z RÚIAN (např. pomocí QGIS pluginu). Tuto vrstvu s vhodně zvolenou obalovou zónou kolem budov použijte pro odstranění velkých rozdílu detekovaných CHM.

Tvorba modelu pro výpočet CHM

Výpočetní proces popsaný výše převedeme do formy modelu. Grafický modelář spustíme a postupně přidáváme výpočetní nástroje .

Výsledný model může vypadat následovně:

Tip

Vytvořené vrstvy lze přidat automaticky do mapového okna z kontextovému menu nad datovou položkou (Display).

Parametrizace modelu

Uživatelské vstupy lze definovat dvěma způsoby:

  • parametrizací jednotlivých parametrů použitých nástrojů v modelu
  • použitím proměnných

Začněme parametrizací nástrojů. U každého parametru je možnost zvolit "Parametrized in model". Pokud parametr vyplníme a přesto políčko zašrtneme, parametr bude parametrizován s definovanou výchozí hodnotou. Vyzkoušejme nastavit parametrizaci u r.in.pdal, a to u input a resolution.

Podobně parametrizujme i ostatní uživatelské vstupy. Po spuštění modelu se objeví dialog umožňující nastavení zvolených parametrů:

Tento způsob parametrizace modelu není v našem případě vhodný. Některé parametry budou duplikovany, jako např. resolution v případě r.in.pdal.

Upravme model s využitím proměnných. V založce Variables přidáme čtyři proměnné:

  • path Cesta k adresáři se vstupními daty
  • tile Označení dlaždice
  • resolution Cílové prostorové rozlišení
  • fill Počet buněk pro výplnění děr

V parametrech se proměnné používají s prefixem %, v případě parametru input nástroje r.in.pdal půjde o %path/dmr5g/%tile.laz. Po spuštění modelu se opět objeví dialog pro nastavení vstupních parametrů:

Úkol

Vyzkoušejte postupně zpracovat všechny stažené dlaždice.

Python API

GRASS nabízí hned několik Python balíčků, podrobnosti naleznete v dokumentaci. Zaměřme se na tři hlavní knihovny související s obsahem tohoto kurzu:

PyGRASS je navržen jako objektově orientované API jazyka Python pro GRASS GIS. To je zásadní rozdíl oproti GRASS Script Package, který se skládá z procedur - funkcí. Je důležité zdůraznit, že PyGRASS není náhradou za GRASS Script Package. Oba balíky existují vedle sebe. Záleží na uživateli, který balíček ve svých skriptech použije. Je také možné oba balíčky zkombinovat v jednom skriptu (nedoporučuje se).

První kroky v jazyce Python provedeme v záložce Python. Syntax do příkazové řádky, např.:

r.fill.stats -k input=KRAV59 output=KRAV59_filled distance=10 --o

Lze jednoduše přepsat do syntaxe Python (zde v syntaxi PyGRASS):

from grass.pygrass.modules import Module

Module('r.fill.stats', flags='k', input='KRAV59', output='KRAV59_filled', distance=10, overwrite=True)

Vylepšení modelu ve formě Python skriptu

Model lze exportovat do formy skriptu v jazyku Python. Před exportem ještě v nastavení modeláře zvolíme PyGRASS API. Výchozím Python API je totiž GRASS Script Package.

Poznámka

Volba nastavení Python API je dostupná od verze GRASS 8.4.

Po přepnutí do záložky Python editor se model převede do Python syntaxe. Tu lze upravovat přímo v této záložce anebo uložit do souboru na disk Save as.

Takto uložený skript lze upravovat ve vašem oblíbeném IDE. V našem případě ale použijme jednoduchý vestavěný Python editor dostupný v nástrojové liště .

Úkol

Upravte skript, aby končil rozumně formulovanou chybou v případě, že vstupní soubor neexsituje.

Upravme skript tak, aby zpracoval všechny dostupné dlaždice najednou. Můžeme volit mezi dvě možnostmi:

  1. využít parametr file nástroje r.in.pdal
  2. pomocí cyklu projít postupně vstupní LAZ soubory, každý zpracovat odděleně a posléze vytvořit mozaiku pomocí r.series (tato cesta by umožňovala zpracování dlaždic paralelizovat pomocí ParallelModuleQueue. Příklad najdete v GRASS Jena Workshop.

Zvolme první variantu. Skript upravíme následovně:

  1. odstraňme parametr skriptu tile
  2. vytvořme funkci collect_files, která pro daný adresář vrátí seznam LAZ souborů.
  3. seznam LAZ souborů předáme do parametru file jako soubor uložený na disku anebo použijeme standardní vstup jako v ukázce níže.
  4. odstraňme duplicitu kódu - vytvořme novou funkci process_product()
  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
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
#!/usr/bin/env python3
#
##############################################################################
#
# MODULE:       model
#
# AUTHOR(S):    martin
#
# PURPOSE:      Script generated by wxGUI Graphical Modeler.
#
# DATE:         Wed Apr  3 22:01:42 2024
#
##############################################################################

# %module
# % description: Script generated by wxGUI Graphical Modeler.
# %end
# %option G_OPT_M_DIR
# % key: path
# % description: Path to input LAZ data
# % required: yes
# % answer: /home/martin/geodata/vyuka/155FGIS/06
# %end
# %option
# % key: resolution
# % description: Target resolution
# % required: yes
# % type: integer
# % answer: 5
# %end
# %option
# % key: fill
# % description: Fill distance (in pixels)
# % required: yes
# % type: integer
# % answer: 10
# %end

import sys
import os
import atexit
import glob

from grass.script import parser, tempfile
from grass.pygrass.modules import Module

def cleanup():
    pass

def collect_files(path):
    fn = tempfile(create=False)
    with open(fn, "w") as fd:
        fd.write("\n".join(glob.glob(os.path.join(path, "*.laz"))))

    return fn

def main(options, flags):
    process_product("dmr5g")
    process_product("dmp1g")

    Module("r.mapcalc.simple",
           overwrite=True,
           expression="A-B",
           a="dmp1g",
           b="dmr5g",
           output="chm")

    Module("r.colors",
           map="chm",
           color="differences",
           offset=0,
           scale=1)

    return 0

def process_product(name):
    Module("r.in.pdal",
           flags='eo',
           overwrite=True,
           file=collect_files(f"{options['path']}/{name}/"),
           output=name,
           method="mean",
           type="FCELL",
           zscale=1.0,
           iscale=1.0,
           resolution=f"{options['resolution']}",
           dimension="z")

    Module("g.region",
           overwrite=True,
           raster=name)

    Module("r.fill.stats",
           flags='k',
           overwrite=True,
           input=name,
           output=name+"_f",
           distance=f"{options['fill']}",
           mode="wmean",
           power=2.0,
           cells=8)

    Module("r.colors",
           map=name+"_f",
           color="elevation",
           offset=0,
           scale=1)

if __name__ == "__main__":
    options, flags = parser()
    atexit.register(cleanup)
    os.environ["GRASS_OVERWRITE"] = "1"
    sys.exit(main(options, flags))

Ukázka CHM (na pozadí ČÚZK WMS Ortofoto):

Přístup k datům

Pro přístup k rastrovým datům lze použít dvě třídy:

  • RasterRow pro čtení a zápis po řádcích
  • RasterSegment pro čtení a zápis po segmentech (dlaždicích)

Nejprve vypišme základní metadata:

from grass.pygrass.raster import RasterRow

with RasterRow('chm') as data:
    print(data.info.cols, data.info.rows)
    min, max = data.info.range
    print(min, max)
    print(max - min)

Dále zkusme reklasifikovat raster pomocí NumPy:

import numpy as np
from grass.pygrass.raster import RasterRow

with RasterRow('chm') as data:
    array = np.array(data)

    array[array < 2] = 1
    array[(array > 2) & (array < 7) ] = 2
    array[array > 7] = 3

    print(np.unique(array, return_counts=True))

Podobně můžeme přístupovat i k vektorovým vrstvám pomocí Vector (bez topologie) a VectorTopo (včetně topologie).

Tip

Další ukázky najdete v materiálech GRASS Jena Workshop: část 13 a část 14