Tvorba DMR a DMP¶

Tato kapitola shromažďuje informace, jak pracovat s daty digitálního modelu reliéfu a povrchu poskytovaných Českým úřadem zeměměřičským a katastrálním v systému GRASS GIS.

Tyto produkty poskytuje ČÚZK ve dvou formátech:

textovém XYZ
binárním LAS, resp. komprimované formě LAZ

Souřadnice X,Y jsou referencovány v souřadnicovém systému S-JTSK (EPSG:5514), souřadnice H (nadmořská výška) ve výškovém referenčním systému Balt po vyrovnání (Bpv).

Poznámka

Tato data nejsou poskytována v režimu otevřených dat, ČÚZK nicméně zveřejnil vzorová data volně ke stažení zde. Ukázkový dataset obsahuje data pouze v textovém formátu XYZ.

Postup importu lidarových dat je popsán v předcházející kapitole.

Digitální model reliéfu¶

Data pro digitální model reliéfu, tj. digitální reprezentaci modelu reliéfu bez umělých a přírodních objektů (např. vegetace nebo budovy) poskytuje ČÚZK ve dvou verzích:

DMR4G - diskrétní body v pravidelné síti (5 x 5 m) bodů o souřadnicích X,Y,H s úplnou střední chybou výšky 0,3 m v odkrytém terénu a 1 m v zalesněném terénu. Další informace zde.
DMR5G - diskrétní body v nepravidelné trojúhelníkové síti (TIN) bodů o souřadnicích X,Y,H s úplnou střední chybou výšky 0,18 m v odkrytém terénu a 0,3 m v zalesněném terénu. Další informace zde.

Digitální model povrchu¶

Data pro digitální model povrchu (DMP), tj. digitální reprezentaci modelu reliéfu včetně umělých a přírodních objektů poskytuje ČÚZK v současnosti v jedné verzi, a to jako:

DMP1G - diskrétní body v nepravidelné sítě výškových bodů (TIN) s úplnou střední chybou výšky 0,4 m pro přesně vymezené objekty (budovy) a 0,7 m pro objekty přesně neohraničené (lesy a další prvky rostlinného pokryvu). Další informace zde.

Následující kapitoly popisují typické postupy tvorby rastrové či TIN reprezentace DMR/DMP.

Tvorba DMR/DMP reprezentace¶

Digitální model reliéfu či povrchu lze v systému GRASS vytvořit několika způsoby. První uvedený postup je vhodný především pro DMR4G. Další tři postupy jsou aplikovatelné jak na DMR5G tak DMP1G.

1. Vytvoření rastrové reprezentace převzetím hodnot vstupních bodových dat

Import vstupních dat do rastrové mapy s prostorových rozlišením odvozeným ze vstupu. Jednoduchý postup vhodný především pro vstupní bodová data v pravidelné mřížce. Příklad v kapitole Vytvoření DMR převzetím hodnot pravidelné mřížky bodových dat DMR4G.

2. Vytvoření rastrové reprezentace s vyplněním „děr“

Import vstupních dat do rastrové mapy s prostorových rozlišením odvozeným ze vstupu. Interpolace hodnot pouze v místech, kde chybí vstupní data. Příklad v kapitole Vytvoření DMR kombinací převzetí hodnot DMR5G a interpolací chybějících hodnot.

3. Spline interpolace z importovaných vektorových bodových dat

Import vstupních dat do vektorové bodové mapy, interpolace výsledného povrchu pomocí modulu v.surf.rst. Časově náročné. Vhodné pro vytvoření DMT s vysokým prostorovým rozlišením. Příklad v kapitole Vytvoření DMR metodou spline na základě DMR5G.

4. Vytvoření TIN reprezentace

Vytvoření TIN (Triangulated irregular network) reprezentace na základě vstupních bodových dat. Většina navazujících nástrojů pro analýzu topografického povrchu v systému GRASS podporuje pouze rastrovou reprezentaci digitálního modelu reliéfu či povrchu. Vhodné při exportu TIN a navazující operace již mimo systém GRASS GIS. Příklad v kapitole Vytvoření DMR v podobě TIN reprezentace.

Vytvoření DMR převzetím hodnot pravidelné mřížky bodových dat DMR4G¶

V prvním kroku zjistíme hraniční souřadnice importovaných dat, viz kapitola import dat.

r.in.xyz -sg input=HLIN04_4g.xyz separator=space

Na základě toho nastavíme výpočetní region a to tak, aby střed buňky odpovídal vstupní pravidelné mřížce bodů. Prostorové rozlišení nastavíme na 5m, což odpovídá vzdálenosti bodů DMR4G.

g.region n=-1088005 s=-1090000 e=-625005 w=-627500 b=461.5 t=554.31
g.region n=n+2.5 s=s-2.5 w=w-2.5 e=e+2.5 res=5

Následně na to data do nastaveného výpočetního regionu naimportujeme.

r.in.xyz input=HLIN04_4g.xyz separator=space output=HLIN04_4g

Poznámka

Podobný postup by bylo možné aplikovat na binární vstupní data ve formátu LAS/LAZ a modul r.in.lidar, viz kapitola import dat.

Obr. 52 Ukázka výsledného produktu digitálního modelu reliéfu vytvořeno na bázi DMR4G převzetím hodnot vstupních bodů. Na obrázku jsou pro ilustraci vykreslena vstupní pravidelná síť bodů DMR4G a výstupní mřížka rastrové mapy.

Tip

Zájmové území by nemělo obsahovat místa bez dat. To můžeme zkontrolovat pomocí modulu r.univar.

r.univar HLIN04_4g

...
total null cells: 0
...

Obr. 53 Vizualizace terénu na bázi DMR4G ve 3D pohledu.

Vytvoření DMR kombinací převzetí hodnot DMR5G a interpolací chybějících hodnot¶

Začneme obdobně jako v předchozí kapitole importem vstupních dat do rastrové reprezentace podle toho, zda jsou vstupní data dostupná v textovém nebo binárním formátu. Příklad níže ukazuje postup pro textová data.

r.in.xyz -sg input=HLIN04_5g.xyz separator=space
g.region n=-1088000.076 s=-1090000.059 e=-624999.829 w=-627499.828 b=461.312 t=554.334

Poznámka

ZÚ poskytuje výškopisná data po dlaždicích SM5 (2x2.5 km). V tomto kontextu by bylo vhodnější načíst data do výpočetního regionu zarovnaného dle kladu SM5, tj. na celá čísla:

g.region n=-1088000 s=-1090000 e=-625000 w=-627500 b=461.312 t=554.334

V našem případě nastavíme prostorové rozlišení na 3 metry (vycházíme z průměrné hustoty vstupních bodů, viz poznámky k importu vektorových dat). Předtím je samozřejmě nutné importovat data do vektorové reprezentace, postup je uveden v kapitole Vytvoření DMR metodou spline na základě DMR5G.

Estimated point density: 0.1102
Estimated mean distance between points: 3.012

Poté data naimportujeme:

r.in.xyz input=HLIN04_5g.xyz separator=space output=HLIN04_5g

Rastrová mapa bude vzhledem k povaze dat DMR5G obsahovat místa bez dat (no-data). V našem případě jde téměř o polovinu buněk(!) Tato místa se bude snažit v dalším kroku vyplnit na základně interpolovaných hodnot.

r.univar map=HLIN04_5g

total null and non-null cells: 557112
total null cells: 249326

Obr. 54 Importovaná dlaždice DMR5G do rastrové reprezentace.

Pro vyplnění chybějících hodnot DMR použijeme modul r.fillnulls. Modul podporuje tři různé metody interpolace: rst (spline, tj. metoda použitá v kapitole Vytvoření DMR metodou spline na základě DMR5G), bilinear (bilinear), bicubic (bicubic). Pro DMR je vhodná metoda spline (rst). Nicnémě vzhledem k její časové náročnosti (desítky minut na uvedené dlaždici) zvolíme v tomto případě akceptovatelnou bikubickou interpolaci (bicubic).

r.fillnulls input=HLIN04_5g output=HLIN04_5g_filled method=bicubic

Obr. 55 Importovaná dlaždice DMR5G po interpolaci chybějících hodnot. Barevná tabulka nastavena na elevation.

Obr. 56 Vizualizace terénu na bázi DMR5G ve 3D pohledu.

Velmi dobrou alternativou k výše uvedenému modulu je r.fill.stats. Ten je určen k rychlému doplnění chybějících hodnot na základě statistiky vstupních dat. Výchozí metoda (wmean) v podstatě odpovídá interpolační metodě IDW, mean poté aplikaci „low-pass filru“. Dále jsou dostupné metody median, mode. Modul je ideální pro zpracovaní velkého objemu dat ve velmi vysokém rozlišení. Na rozdíl od r.fillnulls interpoluje chybějící hodnoty pouze v daném okolí definovaném parametrem cells (výchozí je okolí 8 buněk). Tento postup je vhodný v případě, že místa s chybějícími hodnotami tvoří menší celky. Což není náš případ, kdy počet chybějících hodnot dosahuje téměř 50%.

Vytvoření DMR metodou spline na základě DMR5G¶

Lidarová data importujeme do vektorové reprezentace, postup pro textový a binární formát v kapitole Import lidarových dat. Příklad níže ukazuje import dat v textovém formátu.

v.in.ascii input=HLIN04_5g.xyz output=HLIN04_5g separator=space z=3 -tbz

V našem případě zvolíme prostorové rozlišení 3 metry (odvozené z průměrné hustoty vstupních bodů, viz předchozí kapitola).

g.region vector=HLIN04_5g res=3 -a

Poté spustíme proces interpolace:

v.surf.rst input=HLIN04_5g elevation=HLIN04_5g

Důležité

Modul v.surf.rst patří mezi extrémně výpočetně náročné nástroje. Na testovacím PC trvala interpolace pro výše zmíněná data 30 min.

Od verze GRASS 7.4 podporuje modul paralelizaci výpočtu, což může vést k signifikantnímu zrychlení výpočtu. V našem případě rozložení výpočtu na 12 jader CPU (parametr nprocs=12) vedlo ke snížení výpočetního času na 10 min.

v.surf.rst input=HLIN04_5g elevation=HLIN04_5g nprocs=12

Obr. 57 Ukázka výsledného produktu digitálního modelu reliéfu vytvořeno na bázi DMR5G metodou spline s rozlišením 3m.

Obr. 58 Porovnání vytvořených DMR převzetím hodnot DMR4G při rozlišení 5m a interpolací spline v rozlišení 3m.

Obr. 59 Ukázka výsledného produktu digitálního modelu povrchu vytvořeného spline interpolací v prostorovém rozlišení 2.5m.

Todo

porovnat modely spline vs. bicubic

Vytvoření DMR v podobě TIN reprezentace¶

Vstupní data importujeme do vektorové mapy podobně jako v předcházejícím případě. Pomocí modulu v.delaunay vytvoříme na základě Delaunayho triangulace nepravidelnou trojuhelníkovou síť (TIN). Import vstupních bodů provedeme s menší úpravou, a to bez přepínače -b. Modul v.delaunay totiž vyžaduje pro svůj běh topologické informace vstupních dat.

v.in.ascii input=HLIN04_5g.xyz output=HLIN04_5g separator=space z=3 -tz
v.delaunay input=HLIN04_5g output=HLIN04_5g_TIN

Obr. 60 TIN reprezentace DMR.

Dávkové zpracování dlaždic DMR/DMP a vytvoření výsledné mozaiky¶

ČÚZK poskytuje typicky data ve formě dlaždic, v případě testovacího datasetu se jedná o:

HLIN04,
HLIN05,
HLIN14,
HLIN15.

Zpracování dlaždic můžeme urychlit paralelizací výpočtu, k tomu použijeme nástroje frameworku PyGRASS, viz kapitola PyGRASS. Moduly systému GRASS umožňuje spouštět třída Module, viz příklad v kapitole Úvod do skriptování. Paralelizaci výpočtu je možné poměrně jednoduše implementovat pomocí třídy ParallelModuleQueue, viz dokumentace. Příklad použití si ukážeme na jednoduché operaci importu dat pomocí modulu v.in.ascii:

    import_module = Module('v.in.ascii', separator='space', z=3, flags='tbz',
                           overwrite=overwrite(), quiet=True, run_=False)
    try:
        for f in files:
            basename = os.path.basename(f)
            message("Importing <{}>...".format(f))
            import_task = deepcopy(import_module)
            queue.put(import_task(input=f, output=mapname))
        queue.wait()

Komentáře:

Na řádcích 1-2 vytvoříme objekt modulu v.in.ascii obsahující společné parametry.
Na řádku 7 vytvoříme pro každý importovaný soubor kopii objektu pro import.
Této kopii nastavíme na souboru závislé parametry - cestu k importovanému souboru (input) a název výstupní vektorové mapy (output). Spustíme instanci takto upraveného objektu a zařadíme do fronty (queue jako instance třídy ParallelModuleQueue), viz řádek 8.
```
    queue = ParallelModuleQueue(int(options['nprocs']))
```
Poté necháme všechny paralelně běžící procesy doběhnout, viz řádek 9.

Podobně lze paralelně volat i interpolační modul v.surf.rst pouze s tím rozdílem, že je třeba pro každý proces nastavit příslušný region na základě vstupních dat (dlaždice). Tuto operaci nám zásadně usnadní třída MultiModule, která umožňuje registrovat navazující moduly jako jeden proces ve frontě. V prvním kroku zaregistrujeme modul g.region, pomocí kterého nastavíme výpočetní region dlaždice (11) a poté přidáme do procesu interpolační modul v.surf.rst (12).

    region_module = Module('g.region', n='n+{}'.format(offset), s='s-{}'.format(offset),
                           e='e+{}'.format(offset), w='w-{}'.format(offset),
                           quiet=True)
    rst_module = Module('v.surf.rst', nprocs=rst_nprocs,
                        overwrite=overwrite(), quiet=True, run_=False)
    try:
        for mapname in maps:
            message("Interpolating <{}>...".format(mapname))
            region_task = deepcopy(region_module)
            rst_task = deepcopy(rst_module)
            mm = MultiModule([region_task(vector=mapname),
                              rst_task(input=mapname, elevation=mapname)],
                              sync=False, set_temp_region=True)
            queue.put(mm)
        queue.wait()

Komentáře:

Ve skriptu nastavujeme výpočetní region dlaždice o něco větší tak, aby u výstupních dlaždic DMT vznikly pásy překryvu a bylo možno vytvořit výslednou mozaiku DMT bez ostrých přechodů na hranicích dlaždic, viz offset jako desetinásobek zadaného rozlišení na řádce 1-2.

Poznámka

Třída MultiModule je dostupná od verze GRASS 7.4.

Výsledná mozaika DMT jednotlivých dlaždic může být vytvořena modulem r.series a vhodnou statistickou metodou, viz 109. Výsledný skript může vypadat následovně:

#!/usr/bin/env python

#%module
#% description: Creates DMT from input XYZ data.
#%end
# overwrite: yes
#%option G_OPT_M_DIR
#% required: yes
#%end
#%option
#% key: pattern
#% type: string
#% description: File name pattern to filter files in input directory
#%end
#%option G_OPT_R_ELEV
#% description: Name for output elevation raster map mosaics
#%end
#%option
#% key: resolution
#% description: Output resolution
#% type: double
#%end
#%option
#% key: nprocs
#% description: Number of processes
#% answer: 1
#% type: integer
#%end
#%option
#% key: rst_nprocs
#% description: Number of v.surf.rst processes
#% answer: 1
#% type: integer
#%end

import os
import sys
import time
from copy import deepcopy

from grass.script.core import parser, message, fatal, overwrite

from grass.pygrass.modules import Module, MultiModule, ParallelModuleQueue
from grass.exceptions import CalledModuleError

def import_files(directory, pattern):
    maps = []
    if pattern:
        from glob import glob
        files = glob('{dir}{sep}{pat}'.format(
            dir=directory, sep=os.path.sep, pat=pattern)
        )
    else:
        files = map(lambda x: os.path.join(directory, x),
                    os.listdir(directory)
        )

    start = time.time()

    import_module = Module('v.in.ascii', separator='space', z=3, flags='tbz',
                           overwrite=overwrite(), quiet=True, run_=False)
    try:
        for f in files:
            basename = os.path.basename(f)
            mapname = os.path.splitext(basename)[0]
            maps.append(mapname)
            message("Importing <{}>...".format(f))
            import_task = deepcopy(import_module)
            queue.put(import_task(input=f, output=mapname))
        queue.wait()
    except CalledModuleError:
        return sys.exit(1)

    if not maps:
        fatal("No input files found")

    message("Import finished in {:.0f} sec".format(time.time() - start))

    return maps

def create_dmt_tiles(maps, res, rst_nprocs, offset_multiplier=10):
    offset=res * offset_multiplier

    start = time.time()

    region_module = Module('g.region', n='n+{}'.format(offset), s='s-{}'.format(offset),
                           e='e+{}'.format(offset), w='w-{}'.format(offset),
                           quiet=True)
    rst_module = Module('v.surf.rst', nprocs=rst_nprocs,
                        overwrite=overwrite(), quiet=True, run_=False)
    try:
        for mapname in maps:
            message("Interpolating <{}>...".format(mapname))
            region_task = deepcopy(region_module)
            rst_task = deepcopy(rst_module)
            mm = MultiModule([region_task(vector=mapname),
                              rst_task(input=mapname, elevation=mapname)],
                              sync=False, set_temp_region=True)
            queue.put(mm)
        queue.wait()
    except CalledModuleError:
        return sys.exit(1)

    message("Interpolation finished in {:.0f} min".format((time.time() - start) / 60.))
    
def patch_tiles(maps, output, resolution):
    message("Patching tiles <{}>...".format(','.join(maps)))
    Module('g.region', raster=maps, res=resolution)
    Module('r.series', input=maps, output=output, method='average')
    Module('r.colors', map=output, color='elevation')

def main():
    start = time.time()

    maps = import_files(options['input'], options['pattern'])
    create_dmt_tiles(maps, float(options['resolution']), int(options['rst_nprocs']))
    patch_tiles(maps, options['elevation'], options['resolution'])

    message("Done in {:.0f} min".format((time.time() - start) / 60.))
    
    return 0

if __name__ == "__main__":
    options, flags = parser()

    # queue for parallel jobs
    queue = ParallelModuleQueue(int(options['nprocs']))

    sys.exit(main())

Poznámka

Skript umožňuje paralelizovat jak běh modulů pro import a interpolaci (parametr nprocs), tak modulu v.surf.rst jako takového (rst_nprocs).

Ukázka volání (v tomto případě bude vytíženo při interpolaci 12 jader CPU, výpočet trval přes hodinu a 20 minut):

create-dmt.py input=VYSKOPIS elevation=HLIN_5g pattern=*5g* resolution=3 nprocs=4 rst_nprocs=3

Výsledná verze skriptu ke stažení zde.

Poznámka

aktualizovat screenshot v rozliseni 3m

Obr. 61 Ukázka vizualizace výsledného DMT složeného jako mozaika ze čtyř dlaždic vstupních bodových dat s rozlišením 1m ve 3D potaženého ortofotem (ČÚZK WMS), viz kapitola Vizualizace ve 3D.

Tvorba DMR a DMP¶

Digitální model reliéfu¶

Digitální model povrchu¶

Tvorba DMR/DMP reprezentace¶

Vytvoření DMR převzetím hodnot pravidelné mřížky bodových dat DMR4G¶

Vytvoření DMR kombinací převzetí hodnot DMR5G a interpolací chybějících hodnot¶

Vytvoření DMR metodou spline na základě DMR5G¶

Vytvoření DMR v podobě TIN reprezentace¶

Dávkové zpracování dlaždic DMR/DMP a vytvoření výsledné mozaiky¶

Table of Contents

Vyhledávání