NDVI - Základní verze skriptu

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

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

Poznámka

Spouštět moduly systému GRASS z jazyka Python umožňuje také knihovna GRASS Python Scripting Library. V našem případě budeme striktně používat knihovnu PyGRASS.

#!/usr/bin/env python3

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

# vstup
vis="LC81920252013215LGN00_B4@landsat"
nir="LC81920252013215LGN00_B5@landsat"
aoi="obce@ruian_praha"
# vysledek
ndvi="ndvi"
r_ndvi= "r_{}".format(ndvi)

# 0. nastavit vypocetni region
Module('g.region', align=vis, vector=aoi)

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

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

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

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

print('-' * 80)
for line in report.outputs.stdout.splitlines():
    trida, popisek, procento = line.split(':')
    print("Trida {}: {}".format(trida, procento))
print('-' * 80)

Skript je ke stažení zde.

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

Tip

Namísto funkce print(...) zkuste použít

from grass.pygrass.messages import Messenger

msgr = Messanger()
msgr.message('...')

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.region', align=vis, vector=aoi)

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

g.region align=LC81920252013215LGN00_B4@landsat vector=obce@ruian_praha

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, např. overwrite=True. Běžné přepínače (uvozeny jednou pomlčkou) se předávají jako hodnota argumentu flags.

Poznámka pro pokročilé

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.region al=LC81920252013215LGN00_B4@landsat v=obce@ruian_praha

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

Shortcuts

PyGRASS umožňuje emulovat způsob volání modulů z příkazové řádky pomocí tzv. „shortcuts“. Příklad volání modulu g.region (řádek 15):

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

g.region(raster=vis)

Vstup

Některé moduly přijímají vstup ze souboru, např. r.recode s parametrem rules (řádky 32-33). 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 45.

from subprocess import PIPE
report = Module('r.stats', flags='pl', input=r_ndvi, separator=':', stdout_=PIPE)
for line in report.outputs.stdout.splitlines():
    trida, popisek, procento = line.split(':')