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.

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

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

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 dat v textovém či binárním formátu je popsán v předcházející kapitole.

Digitální model reliéfu či povrchu lze v systému GRASS vytvořit několika způsoby:

1. Vytvoření rastru 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. Postup vhodný především pro pravidelnou mřížku vstupních dat.

2. Spline interpolace z bodových dat

Import vstupních dat do vektorové bodové mapy, interpolace výsledného povrchu pomocí modulu v.surf.rst. Vhodné pro vytvoření DMT s vysokým prostorovým rozlišením.

3. Import to rastrové mapy a vyplnění „děr“

Todo

Doplnit

4. Vytvoření TIN reprezentace

Todo

Doplnit

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.

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

Následující dva příklady vytvoření rastrové reprezentace DMT prezentují dva rozdílné přístupy:

  • vytvoření DMT převzetím hodnot pravidelné mřížky bodových dat DMR4G, prostorové rozlišení je odvozeno ze vstupních dat (tj. 5m)
  • vytvoření vysoce podrobného DMT metodou spline na základě DMR5G

Vytvoření DMT 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
../_images/dmr4g.png

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

Vytvoření vysoce podrobného DMT metodou spline na základě DMR5G

Nastavíme výpočetní region na základě vstupní bodové vektorové mapy, prostorové rozlišení zarovnáme (přepínač -a) na požadované rozlišení, viz kapitola import dat.

Poznámka

Zvolené prostorové rozlišení může vycházet z průměrné vzdálenosti vstupních bodů, tuto statistiku spočítáme pomocí modulu v.lidar.edgedetection.

v.lidar.edgedetection -e input=HLIN04_5g
Estimated point density: 0.1102
Estimated mean distance between points: 3.012

Pokud používáte verzi GRASS 7.2.1 a nižší, je nutné přidat ještě parametr output a to i přestože modul žádný výstup nevytváří.

V našem případě zvolíme prostorové rozlišení 1m.

g.region vector=HLIN04_5g res=1 -a

Poznámka

Pokud používáte GRASS verze 7.2.0 a nižší, tak naimportujte vstupní data včetně topologie (tj. vynechte přepínač -b v případě importních modulů v.in.ascii a v.in.lidar). V opačném případě nebude výše uvedený příkaz pro nastavení výpočetního regionu fungovat.

Poté spustíme proces interpolace:

v.surf.rst input=HLIN04_5g elevation=HLIN04_5g

Důležité

Modul v.surf.rst poskytuje dobré výsledky, bohužel je ale velmi pomalý. Na testovacím PC trvala interpolace pro výše zmíněná data 3hod 14min!

Od verze GRASS 7.3 (aktuální vývojová větev) 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 8 jader CPU (parametr nprocs=8) vedlo ke snížení výpočetního času na 30 min.

v.surf.rst input=HLIN04_5g elevation=HLIN04_5g nprocs=8
../_images/dmr5g.png

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

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

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

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ů (např. vegetace nebo budovy) 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.

Souřadnice X,Y jsou podobně jako v případě DMR4/5G 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).

Produkt DMP vytvoříme obdobně jako v případě DMR interpolací spline:

v.in.ascii input=HLIN04_1g.xyz output=HLIN04_1g separator=space z=3 -tbz
g.region vector=HLIN04_1g res=1 -a
v.surf.rst input=HLIN04_1g elevation=HLIN04_1g nprocs=8

Poznámka

Pokud zpracováváte DMR a DMP současně pro stejné území, je vhodné zachovat stejný výpočetní region. Ten můžete nastavit na základě více vektorových map současně, v tomto případě DMR5G a DMP1G:

g.region vector=HLIN04_1g,HLIN04_5g res=1 -a
../_images/dsm-cuzk.png

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

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
2
3
4
5
6
7
8
9
    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).

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
    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.

Varování

Třída MultiModule je v současnosti dostupná pouze ve vývojové verzi systému GRASS 7.3, viz dokumentace.

Todo

Dopsat scénář řešení v nižších verzích bez třídy MultiModule.

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
  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
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
#!/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=1 nprocs=4 rst_nprocs=3

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

../_images/dmr5g_ortofoto_3d.png

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