NDVI - Pokročilejší verze skriptu

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

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

#%module
#% description: Creates reclassified NDVI based on given AOI.
#%end
#%option G_OPT_R_MAP
#% key: red
#% description: Name of red channel
#% answer: LC81920252013215LGN00_B4@landsat
#%end
#%option G_OPT_R_MAP
#% key: nir
#% description: Name of nir channel
#% answer: LC81920252013215LGN00_B5@landsat
#%end
#%option G_OPT_V_MAP
#% key: aoi
#% answer: obce@ruian_praha
#%end
#%option G_OPT_F_INPUT
#% key: classes
#% required: no
#%end

import sys
from subprocess import PIPE

from grass.script import parser, fatal, message

from grass.pygrass.modules import Module

def main():
    vis = options['red']
    nir = options['nir']
    aoi = options['aoi']
    if '@' in aoi:
        aoi_name = aoi.split('@')[0]
    else:
        aoi_name = aoi
    
    ndvi = "ndvi_{}".format(aoi_name)
    r_ndvi= "r_ndvi_{}".format(aoi_name)

    # 0. nastavit vypocetni region
    Module('g.region', align=vis, vector=aoi)
    
    # 1. vypocet NDVI
    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)
    
    # 2. reklasifikace (1,2,3)
    message("Reklasifikuji...")
    # r.reclass umi reklasifikovat pouze celociselne rastry, proto pouzime
    # r.recode
    if not options['classes']:
        recode = """
        -1:0.05:1 
        0.05:0.35:2 
        0.35:1:3
        """
        kwargs = {
            'rules' : '-',
            'stdin_' : recode
        }
    else:
        kwargs = {
            'rules' : options['classes']
        }
    Module('r.recode', input=ndvi, output=r_ndvi,
           overwrite=True, **kwargs)
    
    # 3. tabulka barev
    if options['classes']:
        colors = """
        1 red
        2 yellow
        3 0 136 26
        """
        Module('r.colors', map=r_ndvi,
               rules='-', quiet=True, stdin_=colors)
    
    # 4. vypsat vysledek
    message("Generuji report...")
    report = Module('r.stats', flags='pl', input=r_ndvi, separator=':', stdout_=PIPE)
    
    for line in report.outputs.stdout.splitlines():
        trida, popisek, procento = line.split(':')
        print("Trida {}: {}".format(trida, procento))
        
    return 0

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

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

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

Obr. 11 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. 12 Výsledek je vypsán do záložky Command output v dialogu nástroje.

Poznámky k uživatelskému rozhraní

Parametry jsou definovány pomocí tzv. standardizovaných voleb, např. G_OPT_R_MAP definující parametr pro volbu rastrové mapy. V našem případě změníme potřebné vlastnosti (key, description) a zvolíme výchozí hodnotu parametru pro snažší testování (answer).

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

Creates reclassified NDVI based on given AOI.
Usage:
 ndvi-v2.py red=name nir=name aoi=name [classes=name]
[--help]
   [--verbose] [--quiet] [--ui]
Parameters:
      red   Name of red channel
            default: LC81920252013215LGN00_B4@landsat
      nir   Name of nir channel
            default: LC81920252013215LGN00_B5@landsat
      aoi   Name of vector map
            default: obce@ruian_praha
  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("VIS: {0} ; NIR: {1}".format(vis, nir))

přepsáno na

    message("VIS: {0} ; NIR: {1}".format(vis, nir))

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 red=LC81920252013215LGN00_B4@landsat nir=LC81920252013215LGN00_B5@landsat aoi=obce@ruian_praha --q

Trida 1: 1.30%
Trida 2: 72.33%
Trida 3: 26.37%

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. 13 Vygenerovaný grafický dialog skriptu.

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 red=LC81920252013215LGN00_B4@landsat nir=LC81920252013215LGN00_B5@landsat aoi=obce@ruian_praha classes=classes.txt

s výsledkem:

Trida 1: 4.70%
Trida 2: 79.37%
Trida 3: 15.93%