Automatizace v GRASS GIS, PyGRASS

Cvičení

Vstupní data

• Zeměměřičský úřad
  • DMR5G

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

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.

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:

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 - KRAV59_dmp"

Příklad CMH v tabulce barev differences:

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ě:

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ů:

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:

• GRASS Script Package jako vstupní brána do skriptovacích možností jazyka Python v GRASS GIS
• pyGRASS moderní objektově orientované Python API
• Temporal Framework Python API pro práci s časoprostorovými daty

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.

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ě .

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).