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:

  1. textovém XYZ

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

../_images/dmr4g.png

Obr. 51 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
...
../_images/dmr4g-3d.png

Obr. 52 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

g.region res=3 -a

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
../_images/dmr5g-holes.png

Obr. 53 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
../_images/dmr5g-filled.png

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

../_images/dmr5g-3d.png

Obr. 55 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
../_images/dmr5g.png

Obr. 56 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.

../_images/dmr4g-vs-5g.png

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

../_images/dsm-cuzk.png

Obr. 58 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
../_images/tin.png

Obr. 59 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:

1    import_module = Module('v.in.ascii', separator='space', z=3, flags='tbz',
2                           overwrite=overwrite(), quiet=True, run_=False)
3    try:
4        for f in files:
5            basename = os.path.basename(f)
6            message("Importing <{}>...".format(f))
7            import_task = deepcopy(import_module)
8            queue.put(import_task(input=f, output=mapname))
9        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).

 1    region_module = Module('g.region', n='n+{}'.format(offset), s='s-{}'.format(offset),
 2                           e='e+{}'.format(offset), w='w-{}'.format(offset),
 3                           quiet=True)
 4    rst_module = Module('v.surf.rst', nprocs=rst_nprocs,
 5                        overwrite=overwrite(), quiet=True, run_=False)
 6    try:
 7        for mapname in maps:
 8            message("Interpolating <{}>...".format(mapname))
 9            region_task = deepcopy(region_module)
10            rst_task = deepcopy(rst_module)
11            mm = MultiModule([region_task(vector=mapname),
12                              rst_task(input=mapname, elevation=mapname)],
13                              sync=False, set_temp_region=True)
14            queue.put(mm)
15        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ě:

  1#!/usr/bin/env python
  2
  3#%module
  4#% description: Creates DMT from input XYZ data.
  5#%end
  6# overwrite: yes
  7#%option G_OPT_M_DIR
  8#% required: yes
  9#%end
 10#%option
 11#% key: pattern
 12#% type: string
 13#% description: File name pattern to filter files in input directory
 14#%end
 15#%option G_OPT_R_ELEV
 16#% description: Name for output elevation raster map mosaics
 17#%end
 18#%option
 19#% key: resolution
 20#% description: Output resolution
 21#% type: double
 22#%end
 23#%option
 24#% key: nprocs
 25#% description: Number of processes
 26#% answer: 1
 27#% type: integer
 28#%end
 29#%option
 30#% key: rst_nprocs
 31#% description: Number of v.surf.rst processes
 32#% answer: 1
 33#% type: integer
 34#%end
 35
 36import os
 37import sys
 38import time
 39from copy import deepcopy
 40
 41from grass.script.core import parser, message, fatal, overwrite
 42
 43from grass.pygrass.modules import Module, MultiModule, ParallelModuleQueue
 44from grass.exceptions import CalledModuleError
 45
 46def import_files(directory, pattern):
 47    maps = []
 48    if pattern:
 49        from glob import glob
 50        files = glob('{dir}{sep}{pat}'.format(
 51            dir=directory, sep=os.path.sep, pat=pattern)
 52        )
 53    else:
 54        files = map(lambda x: os.path.join(directory, x),
 55                    os.listdir(directory)
 56        )
 57
 58    start = time.time()
 59
 60    import_module = Module('v.in.ascii', separator='space', z=3, flags='tbz',
 61                           overwrite=overwrite(), quiet=True, run_=False)
 62    try:
 63        for f in files:
 64            basename = os.path.basename(f)
 65            mapname = os.path.splitext(basename)[0]
 66            maps.append(mapname)
 67            message("Importing <{}>...".format(f))
 68            import_task = deepcopy(import_module)
 69            queue.put(import_task(input=f, output=mapname))
 70        queue.wait()
 71    except CalledModuleError:
 72        return sys.exit(1)
 73
 74    if not maps:
 75        fatal("No input files found")
 76
 77    message("Import finished in {:.0f} sec".format(time.time() - start))
 78
 79    return maps
 80
 81def create_dmt_tiles(maps, res, rst_nprocs, offset_multiplier=10):
 82    offset=res * offset_multiplier
 83
 84    start = time.time()
 85
 86    region_module = Module('g.region', n='n+{}'.format(offset), s='s-{}'.format(offset),
 87                           e='e+{}'.format(offset), w='w-{}'.format(offset),
 88                           quiet=True)
 89    rst_module = Module('v.surf.rst', nprocs=rst_nprocs,
 90                        overwrite=overwrite(), quiet=True, run_=False)
 91    try:
 92        for mapname in maps:
 93            message("Interpolating <{}>...".format(mapname))
 94            region_task = deepcopy(region_module)
 95            rst_task = deepcopy(rst_module)
 96            mm = MultiModule([region_task(vector=mapname),
 97                              rst_task(input=mapname, elevation=mapname)],
 98                              sync=False, set_temp_region=True)
 99            queue.put(mm)
100        queue.wait()
101    except CalledModuleError:
102        return sys.exit(1)
103
104    message("Interpolation finished in {:.0f} min".format((time.time() - start) / 60.))
105    
106def patch_tiles(maps, output, resolution):
107    message("Patching tiles <{}>...".format(','.join(maps)))
108    Module('g.region', raster=maps, res=resolution)
109    Module('r.series', input=maps, output=output, method='average')
110    Module('r.colors', map=output, color='elevation')
111
112def main():
113    start = time.time()
114
115    maps = import_files(options['input'], options['pattern'])
116    create_dmt_tiles(maps, float(options['resolution']), int(options['rst_nprocs']))
117    patch_tiles(maps, options['elevation'], options['resolution'])
118
119    message("Done in {:.0f} min".format((time.time() - start) / 60.))
120    
121    return 0
122
123if __name__ == "__main__":
124    options, flags = parser()
125
126    # queue for parallel jobs
127    queue = ParallelModuleQueue(int(options['nprocs']))
128
129    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

../_images/dmr5g_ortofoto_3d.png

Obr. 60 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.