Python

Základní verze skriptu

Začneme verzí skriptu bez vstupní parametrů, názvy vstupních a výstupních rastrových map jsou uvedeny na řádcích 7-12.

Na řádku 3 importuje z knihovny PyGRASS třídu Module, která nám umožní z prostředí jazyka Python spouštět moduly systému GRASS jako je např. g.mapsets (viz řádek 7) a další.

Poznámka

Spouštět moduly systému GRASS z jazyka Python umožňuje také knihovna GRASS Python Scripting Library. V našech případech budeme ale použít pro tento účel PyGRASS a jeho třídu Module.

 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
#!/usr/bin/env python

from grass.pygrass.modules import Module
from subprocess import PIPE

# vstup
vis="LC81920252013215LGN00_B4"
nir="LC81920252013215LGN00_B5"
mapset="landsat"
# vysledek
ndvi="ndvi"
r_ndvi= "r_{}".format(ndvi)

# pridat mapset do vyhledavaci cesty
Module('g.mapsets', mapset=mapset, operation='add', quiet=True)

# nastavit vypocetni region
Module('g.region', raster=vis)

# vypocet NDVI
print ("VIS: {0} ; NIR: {1}".format(vis, nir))
Module('r.mapcalc',
       expression="{o} = float({n} - {v}) / ({n} + {v})".format(o=ndvi, v=vis, n=nir),
       overwrite=True)

# reklasifikace (1,2,3)
print ("Reklasifikuji...")
# r.reclass umi reklasifikovat pouze celociselne rastry, proto pouzime
# r.recode
recode = """
 -1:0.05:1 
 0.05:0.35:2 
 0.35:1:3
"""
Module('r.recode', input=ndvi, output=r_ndvi,
       rules='-', overwrite=True, stdin_=recode)

# popisky
labels = """
1:bez vegetace, vodni plochy
2:plochy s minimalni vegetaci
3:plochy pokryte vegetaci
"""
Module('r.category', map=r_ndvi,
       separator=':', rules='-', stdin_=labels)

# tabulka barev
colors = """
 1 red
 2 yellow
 3 0 136 26
"""
Module('r.colors', map=r_ndvi,
       rules='-', quiet=True, stdin_=colors)

# vypsat vysledek
print ("Generuji report...")
report = Module('r.stats', flags='pl', input=r_ndvi, separator=':', stdout_=PIPE)

print ('-' * 80)
for trida, label, procento in map(lambda x: x.split(':'), report.outputs.stdout.splitlines()):
    print ("Trida {0} ({1:28s}): {2:>7}".format(trida, label, procento))
print ('-' * 80)

print ("Hotovo!")

Skript je ke stažení zde.

../_images/wxgui-ndvi-v1.png

Obr. 7 Příklad spuštění skriptu v GUI a vizualizace výsledku v mapovém okně.

Poznámky k volání modulů

Moduly systému GRASS se volají ve skriptech se stejnými parametry jako z příkazové řádky. Například pro volání na řádku 15:

Module('g.mapsets', mapset=mapset, operation='add', quiet=True)

by korespondující zápis pro příkazovou řádku vypadal následovně:

g.mapsets mapset=landsat operation=add --quiet

Jednotlivé parametry modulu se zadávaní jako argumenty třídy Module. Vyjímkou jsou globální přepínače (tj. ty, které jsou uvozeny dvěma pomlčkami) jako je --quiet, --overwrite a další. Ty se zadávají jako argument s hodnotou True. V tomto případě tedy quiet=True. Běžné přepínače (uvozeny jednou pomlčkou) se předávají jako hodnota argumentu flags.

Poznámka

Zkracování názvů parametrů

Při volání modulů z příkazové řádky lze názvy parametrů libovolně zkracovat, pouze s tou podmínkou, aby byly jednoznačné. V níže uvedeném případě bude následnující volání v pořádku i když méně čitelné.

g.mapsets landsat o=add --q

Podobné zkracování názvů parametrů není při použití třídy Module z knihovny PyGRASS možné.

Poznámka

Shortcuts

PyGRASS umožňuje emulovat způsob volání modulů podobně jako POSIX přes tzv. „shortcuts“. Příklad volání modulu g.mapsets (řádek 15):

from grass.pygrass.modules.shortcuts import general as g

g.mapsets(mapset=mapset, operation='add', quiet=True)

Vstup

Některé moduly přijímají vstup ze souboru, např. r.recode s parametrem rules (řádky 30-36). Místo fyzického vytvoření vstupního souboru na disku lze použít standardní vstup, konktrétně argument stdin_ s hodnotou řetězce, který má být na vstupu. V tomto případě musí parameter modulu rules nabývat hodnoty -.

recode = """
 -1:0.05:1 
 0.05:0.35:2 
 0.35:1:3
"""
Module('r.recode', input=ndvi, output=r_ndvi,
       rules='-', overwrite=True, stdin_=recode)

Zpracování výstupu

U modulů, které svůj výstup zapisují na standardní výstup, lze jejich výstup zachytit přes argument stdout_=PIPE. Obsah výstupu je potom uložen jako řetězec v atributu třídy Module outputs.stdout, viz řádek 4,58,61-62.

from subprocess import PIPE
report = Module('r.stats', flags='pl', input=r_ndvi, separator=':', stdout_=PIPE)
for trida, label, procento in map(lambda x: x.split(':'), report.outputs.stdout.splitlines()):
    print ("Trida {0} ({1:28s}): {2:>7}".format(trida, label, procento))

Pokročilejší verze skriptu

Pokročilejší verze skriptu je rozšířena o:

  • uživatelské rozhraní, řádky 3-19
  • vstupní parametry (mapset, output_postfix a classes), viz řádky 6-19
  • uživatelské rozhraní je zpracováno funkcí parse() (řádek 102), která je součástí balíčku grass.script (řádek 24)
  • hodnoty parametrů jsou na řádku 102 uloženy do proměnné options, přepínače do proměnné flags, ty jsou dále použity na řádcích 30-31,57
  • celý kód je vložen do funkce main() (řádek 29)
  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
#!/usr/bin/env python

#%module
#% description: Creates reclassified NDVI.
#%end
#%option G_OPT_M_MAPSET
#% required: yes
#% answer: landsat
#%end
#%option
#% key: output_postfix
#% type: string
#% description: Postfix for output maps
#% answer: ndvi
#%end
#%option G_OPT_F_INPUT
#% key: classes
#% required: no
#%end

import sys
from subprocess import PIPE

import grass.script as grass

from grass.pygrass.modules import Module
from grass.pygrass.gis import Mapset

def main():
    mapset = options['mapset']
    ndvi = options['output_postfix']
    
    # pridat mapset do vyhledavaci cesty
    Module('g.mapsets', mapset=mapset, operation='add', quiet=True)
    
    try:
        vis = Mapset(mapset).glist('raster', pattern='*B4')[0]
        nir = Mapset(mapset).glist('raster', pattern='*B5')[0]
    except IndexError:
        grass.fatal("Nelze najit vstupni kanaly")
    r_ndvi= "r_{}".format(ndvi)
    
    # nastavit vypocetni region
    Module('g.region', raster=vis)
    
    # vypocet NDIV
    grass.message("VIS: {0} ; NIR: {1}".format(vis, nir))
    Module('r.mapcalc',
           expression="{o} = float({n} - {v}) / ({n} + {v})".format(o=ndvi, v=vis, n=nir),
           overwrite=True)
    
    # reklasifikace (1,2,3)
    grass.message("Reklasifikuji...")
    # r.reclass umi reklasifikovat pouze celociselne rastry, proto pouzime
    # r.recode
    args = {}
    if options['classes']:
        args['rules'] = options['classes']
    else:
        recode = """
-1:0.05:1 
0.05:0.35:2 
0.35:1:3
        """
        args['rules'] = '-'
        args['stdin_'] = recode
    Module('r.recode', input=ndvi, output=r_ndvi,
           overwrite=True, **args)
    
    # popisky
    labels = """
1:bez vegetace, vodni plochy
2:plochy s minimalni vegetaci
3:plochy pokryte vegetaci
"""
    Module('r.category', map=r_ndvi,
           separator=':', rules='-', stdin_=labels)
    
    # tabulka barev
    colors = """
1 red
2 yellow
3 0 136 26
"""
    Module('r.colors', map=r_ndvi,
           rules='-', quiet=True, stdin_=colors)
    
    # vypsat vysledek
    grass.message("Generuji report...")
    report = Module('r.stats', flags='pl', input=r_ndvi, separator=':', stdout_=PIPE)
    
    print('-' * 80)
    for trida, label, procento in map(lambda x: x.split(':'), report.outputs.stdout.splitlines()):
        print("Trida {0} ({1:28s}): {2:>7}".format(trida, label, procento))
    print('-' * 80)

    grass.message("Hotovo!")

    return 0

if __name__ == "__main__":
    options, flags = grass.parser()
    sys.exit(main())

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

../_images/wxgui-ndvi-v2-0.png

Obr. 8 Příklad spuštění pokročilé verze skriptu v GUI, výběr vstupních parametru v dialogu nástroje.

../_images/wxgui-ndvi-v2-1.png

Obr. 9 Výsledek je vypsán do záložky Command output v dialogu nástroje.

Poznámky k uživatelskému rozhraní

Příklad parametru output_postfix

který definuje jeho

  • název (key)
  • popisek (description)
  • výchozí hodnotu (answer)
  • a typ parametru (type)

U dalších parametrů jsou použity tzv. standardizované volby, např. G_OPT_M_MAPSET definuje parametr pro volbu mapsetu. V našem případě nastavíme parametr skriptu jako povinný (required: yes) a doplníme výchozí volbu (mapset landsat), viz Obr. 8.

Ve výsledku se skript chová jako standardní modul systému GRASS, přepínačem --help obdržíme informace o jeho syntaxi.

ndvi-v2.py --help
Description:
 Creates reclassified NDVI.

Usage:
 ndvi-v2.py mapset=name [output_postfix=string] [classes=name] [--help]
   [--verbose] [--quiet] [--ui]

Flags:
 --h   Print usage summary
 --v   Verbose module output
 --q   Quiet module output
 --ui  Force launching GUI dialog

Parameters:
          mapset   Name of mapset (default: current search path)
                    '.' for current mapset
                   default: landsat
  output_postfix   Postfix for output maps
                   default: ndvi
         classes   Name of input file

Poznámky k vypisování informačních zpráv

Nahradili jsme funkci print() pro vypisování zpráv o průběhu funkcí message() z balíčku grass.script.

print ("Reklasifikuji...")

přepsáno na

    grass.message("Reklasifikuji...")

Díky tomu budou fungovat globální přepínače --quiet a --verbose pro tichý, resp. upovídaný mód. Např. při použítí volby --quiet se vypíše pouze výsledný report, ostatní zprávy o průběhu výpočtu budou skryty.

ndvi-v2.py mapset=landsat --q
--------------------------------------------------------------------------------
Trida 1 (bez vegetace, vodni plochy  ):   0.28%
Trida 2 (plochy s minimalni vegetaci ):  30.24%
Trida 3 (plochy pokryte vegetaci     ):  21.00%
Trida * (no data                     ):  48.49%
--------------------------------------------------------------------------------

Poznámky k hledání vstupních rastrových dat

Pro nalezení rastrových map končících na B4 a B5 použijeme funkci glist třídy Mapset. Třída Mapset je součástí balíčku pygrass.gis.

Poznámka

Blokem try/except zachytíme chybu v případě, že rastrové mapy nebudou nalezeny. Potom zavoláme funkci fatal() z knihovny grass.script, která skript ukončí.

Poznámky ke spuštění modulu

Pokud skript spustíme bez parametrů mělo by vyskočit grafické okno podobné ostatním modulům systému GRASS.

../_images/wxgui-ndvi-v2-0.png

Obr. 10 Vygenerovaný grafický dialog skriptu.

Poznámka

Spuštění skriptu s parametrem classes

Mějme soubor classes.txt s odlišným rozdělením tříd:

-1:0.1:1 
0.1:0.40:2 
0.40:1:3

Soubor ke stažení zde.

Spuštění skriptu bude vypadat následovně

ndvi-v2.py map=landsat classes=classes.txt

s výsledkem:

 --------------------------------------------------------------------------------
 Trida 1 (bez vegetace, vodni plochy  ):   0.68%
 Trida 2 (plochy s minimalni vegetaci ):  38.43%
 Trida 3 (plochy pokryte vegetaci     ):  12.39%
 Trida * (no data                     ):  48.49%
--------------------------------------------------------------------------------