Každý kdo pracuje s prostorovými daty, se dříve nebo později dostane do situace, kdy pro dosažení požadovaného výsledku používá posloupnost dílčích nástrojů (např. Výběr prvků na základě atributu -> ořezaní vrstvy -> obalová zóna). V praxi se často stává, že stejný postup chceme aplikovat na různá vstupní data, nebo (v horším případě) musíme stejný postup opakovat několikrát na původních vstupních datech, např. z důvodu jejich drobných úprav v průběhu projektu. V takovém případě je vhodné naučit se proces automatizovat – vytvořit si vlastní nástroj, který bude tuto posloupnost operací opakovat za nás. Jelikož celý automatizovaný pracovní postup potom můžeme provést spuštěním jediného nástroje, výhodu takového přístupu pocítíme velmi brzo.

Většina úloh má více možných řešení, z toho důvodu je velice důležité mít všeobecný přehled o dostupných možnostech zpracování dat. Čím větší přehled totiž máme, tím dokážeme sestrojit lepší nástroj (jednodušší, efektivnější, rychlejší, robustnější, …). K automatizaci při zpracování prostorových dat lze využít různých nástrojů, programovacích jazyků a jejich kombinací (Python, R, GRASS, GDAL, Bash, …). Jestliže však potřebujeme vytvořit nástroj, který bude uživatelsky přívětivý, ať už pro naše potřeby nebo pro potřeby ostatních uživatelů, můžeme využít tvorbu nástrojů – uživatelských skriptů – přímo v QGIS.

Framework Zpracování – Processing framework

Jednotlivé nástroje, jak nativní tak od dalších poskytovatelů (GRASS, SAGA, R, …) lze v QGIS spouštět pomocí frameworku Zpracování (Processing framework). Hlavní součástí frameworku je panel Nástroje zpracování (Processing toolbox), který umožňuje spouštět nástroje jednotlivě i hromadně jako dávkový proces. Při hromadném zpracování (Dávkový proces – Batch process) lze spustit najednou několik opakování algoritmu s různými vstupními hodnotami (odlišné vstupní vrstvy, jiná hodnota vzdálenosti obalové zóny atd.). Použití a nastavení frameworku Zpracování jsme podrobněji popsali v materiálech Školení QGIS pro pokročilé.

Pro automatizaci pracovního postupu je nejpodstatnější částí frameworku možnost tvorby vlastních nástrojů – vlastních algoritmů. Zde se nabízí dvě základní možnosti, tvorba nástroje pomocí grafického modelu nebo přímo tvorba nástroje jako skriptu.

Obr.1 Panel Nástroje zpracování

Tvorba grafických modelů

Nejjednodušší možnost jak automatizovat procesy v QGIS je tvorba vlastních nástrojů s využitím grafického modeláře (Graphical modeler). Zde můžeme řetězit či kombinovat funkce z Nástrojů zpracování, včetně vlastních nástrojů – uživatelských skriptů. Výsledkem grafického modelu je pak vlastní spustitelný nástroj (tvorbě modelu je také věnována kapitola v materiálech Školení QGIS pro pokročilé). Nicméně při komplikovanějších požadavcích má tvorba grafických modelů své limity. Ty lze částečně vyřešit vložením malých uživatelských skriptů, ale v případě, že bychom tyto skripty měli psát sami (nebo už máme základní znalosti jazyka Python), je v takové situaci lepší od grafického modelu přejít přímo k tvorbě nástroje jako uživatelského skriptu. Dobrá zpráva je, že grafický modelář umožňuje model exportovat jako Python skript, který lze rovnou použít pro vytvoření nástroje ze skriptu. Pokud tedy zjistíme nějaké omezení v průběhu tvorby modelu, nebo chceme nějaký stávající model upravit, můžeme model exportovat a doladit skript podle našich potřeb.

Obr. 2 Příklad grafického modelu

Tvorba nového nástroje jako Python skriptu

poznámka: Příspěvek se vztahuje k verzi QGIS 2.18 a starším, ve verzi QGIS 3 došlo k významným změnám – přechod na Python 3 a úpravy PyQGIS API. Proto nebude uvedený postup od verze QGIS 3 fungovat.

Skript pro nový nástroj můžeme vytvořit nebo importovat ze souboru .py přímo v panelu Nástroje zpracování. Pokud chceme vytvořit nový skript přímo v QGIS, poklikáním na položku Skripty ->Tools -> Create new script v panelu Nástroje zpracování se nám otevře nativní QGIS textový editor. Skripty můžeme editovat přímo v tomto okně, nebo lze použít jiný oblíbený textový editor.

Obr. 3 Okno QGIS textového editoru

Jelikož vytvořený skript budeme spouštět jako běžný nástroj, je v první řadě nutné provázat jednotlivé vstupní proměnné s GUI nástrojů zpracování. Tyto informace spolu s názvem skriptu a zařazením skriptu do definované skupiny se uvádí na prvních řádcích. Řádky s těmito informacemi se vyznačují dvojitým křížkem (##) na začátku řádku, a uvádí se v následující struktuře:

##[nazev_parametru]=[typ_parametru] [volitelne_hodnoty]

příklad:

##Nazev noveho nastroje=name

##ciselny parametr=number
##rozsah=extent
##vstupni_vrstva=vector
##vystupni_vrstva=output raster

Obr. 4 Vzhled okna vytvořeného nástroje

Samotné algoritmy potom zapisujeme do skriptu podobným způsobem jako se používají v QGIS Python konzole, tedy pomocí příkazu processing.runalg(). Jako první parametr se uvádí název funkce, dále potom jednotlivé parametry daného algoritmu. Názvy algoritmů můžeme zjistit v Python konzole pomocí metody .alglist(), nebo pokud jsme chtěný algoritmus v minulosti spouštěli z nástrojů zpracování, lze Python příkaz pro spuštění najít v historii (Menu ->Zpracování -> Historie…). Nápovědu k jednotlivým algoritmům lze v konzole zobrazit pomocí metody .alghelp() např. processing.alghelp(“qgis:creategrid”). Výsledek algoritmu, můžeme uložit do proměnné, v takovém případě jako parametr výstupu použijeme None, takto:

grid = processing.runalg("qgis:creategrid",1,Extent,Resolution,Resolution,epsg,None)

pokud bychom chtěli zavolat výstupní data z předchozího příkladu použijeme:

grid['OUTPUT']

pokud bychom chtěli použít algoritmus pro generování výsledné vrstvy skriptu příkaz by vypadal takto (za předpokladu, že máme nadefinováno ##vystupni_vrstva = output vector v hlavičce):

grid = processing.runalg("qgis:creategrid",1,Extent,Resolution,Resolution,epsg,vystupni_vrstva)

Příkazy, typy parametrů, jejich hodnoty a další potřebné informace najdeme v dokumentaci QGIS zde.

Příklad – Tvorba rastru na základě počtů bodů a hodnot atributů bodové vrstvy

Pro ilustraci si uvedeme příklad, který vychází z reálného dotazu jednoho z účastníků školení. Řekněme, že pro náš projekt potřebujeme nástroj, který vytvoří rastrovou vrstvu na základě bodové vrstvy a hodnot jejích atributů. Chceme docílit toho, aby pixely výsledného rastru nesly jednu z volitelných hodnot např. součet, minimum, maximum či průměr daného atributu nebo součet bodů spadajících do jednoho pixelu.

Jelikož žádná takováto funkce v QGIS není, musíme navrhnout sled jednotlivých úkonů, které povedou k výsledku. Protože máme základní přehled o funkcích QGIS, víme, že pomocí základních funkcí lze zjistit počet bodů spadajících do polygonu jiné vrstvy (funkce Spočítat body v polygonu – Counts point in polygon) nebo připojit atribut vektorové vrstvy na základě prostorového dotazu a přitom spočítat základní statistiku např. průměrnou hodnotu atributu (Připojit atributy podle umístění – Join atributes by location). Dále víme, že lze snadno tvořit vektorové mřížky – polygonovou síť (Vytvořit mřížku – Create grid) a převádět vektorové vrstvy na rastr (Převézt na rastr (vektor na rastr)… – Rasterize (vector to raster)) se zachováním jednoho číselného atributu jako hodnoty pixelů. Správnou posloupností těchto funkcí máme řešení na světě.

Počty bodů

Začneme jednodušší částí, a to vytvořením rastru s počtem bodů v oblasti pixelu. Základem bude vytvoření vektorové mřížky v dané oblasti (rozsah/extent) o požadovaných rozměrech čtverců (velikost pixelu budoucího rastru).

##Pocet bodu do rastru=name

##Rozliseni=number
##Rozsah=extent
##Bodova_vrstva=vector
##Rastrovy_vystup=output raster

#vytvoreni mrizky s atributem poctu bodu
mrizka = processing.runalg("qgis:creategrid",1,Rozsah,Rozliseni,Rozliseni,"EPSG:5514",None)
pocet = processing.runalg('qgis:countpointsinpolygon',mrizka['OUTPUT'],Bodova_vrstva,'NUMPOINTS',None)

#prevod na rastr
processing.runalg("gdalogr:rasterize",pocet['OUTPUT'],'NUMPOINTS',1,Rozliseni,Rozliseni,Rozsah,False,5,"",4,75.0,6.0,1.0,False,0,"",Rastrovy_vystup)

Takový nástroj bychom byli schopni vytvořit i pomocí modeleru, ale už i u takto jednoduchého modelu můžeme spatřit jedno omezení, a to v nutnosti pevného nastavení souřadnicového systému vytvářené mřížky (v příkladu – mrizka = processing.runalg…,”EPSG:5514″,…). Takové řešení bude komplikovat práci, pokud budeme chtít vytvářet výstupy s různým souř. systémy, proto už je nutné zde zapojit skriptování.

V uživatelském skriptu bychom takový problém snadno vyřešili vytvořením textové proměnné (např. ##EPSG_kod = string, v příkladu by potom bylo mrizka = processing.runalg…,EPSG_kod,…), kdy uživatel při spuštění nástroje zadá ručně EPSG kód (např. EPSG:5514). U našeho příkladu však vytváříme data na základě existující vrstvy, takže výslednou mřížku chceme ve stejném souř. systému jako vstupní bodová vrstva. Proto by bylo nejjednodušší, aby se tato informace vzala přímo ze vstupní bodové vrstvy a my se o tento parametr dále nemuseli starat. Toto v modeleru také není možné, v uživatelském skriptu si s tím však jednoduše poradíme přidáním pár řádků:

from qgis.core import QgsVectorLayer

PointLayer = QgsVectorLayer(Point_layer, "layer_name", "ogr")
epsg = PointLayer.crs().authid()

...
mrizka = processing.runalg("qgis:creategrid",1,Rozsah,Rozliseni,Rozliseni,epsg,None)
...

V této fázi máme hotový skript na vytvoření rastru s počtem odpovídajích bodů ležící v pixelu.

Základní statistika

Druhá, komplikovanější část úlohy spočívá ve vytváření dalšího rastru, jehož pixely budou mít hodnotu vypočítanou z atributu bodové vrstvy. Princip zůstává stejný, jen místo funkce pro počet bodů použijeme funkci Připojit atributy podle umístění (Join atributes by location) pro připojení atributů bodové vrstvy k odpovídajícím polygonům mřížky. Tato funkce vypočítá z atributu zvolenou statistiku (např. suma, minimum, maximum či průměr). Pokud bychom tento postup chtěli využít v modeleru, objeví se několik dalších komplikací, při jejichž řešení se bez skriptování neobejdeme (např. výběr statistiky z rolovací nabídky a její přidělení hodnoty vytvořeného atributu pixelům rastru).

Výhodou také je, že ve skriptech můžeme definovat různé podmínky posloupností funkcí, např. když jedna funkce z nějakého důvodu nefunguje, využije se funkce jiná atd. Ve výsledném skriptu našeho příkladu se vyzkouší funkce z knihovny GDAL a teprve pokud by výpočet selhal, použije se funkce z GRASS. Také je zde ukázka toho, jak lze získat rozsah (extent) z vrstvy mřížky vytvořené během procesu a použít jej jako parametr pro vytvoření rastru (na toto je i v modeleru příslušná funkce). Tento postup definování rozsahu volíme v tomto případě z důvodu různého přístupu při tvoření rastru a vektorové mřížky při zadávání rozsahu. Rastr vytvoří pouze ty pixely které by spadaly větší polovinou do zvoleného rozsahu, kdežto mřížka se vytvoří tak, aby zvolený rozsah byl vždy “zaplněn” polygony, v takovém případě má vrstva zpravidla větší rozsah než námi zvolený (Obr. 5). Pokud tedy chceme, aby rastr odpovídal mřížce, tedy aby nebyl ve skutečnosti menší než námi zvolený rozsah, můžeme použít zmíněný postup.

Obr. 5a Ruční výběr rozsahu

Obr. 5b Výsledek použití vybraného rozsahu. Výsledná rastrová vrstva má rozsah menší než definovaný, výsledná vektorová vrstva (mřížka) má rozsah větší než definovaný

Points to raster statistics=name
##Raster=group
##Resolution=number
##Extent=extent
##Point_layer=vector
##Attribute_values=field Point_layer
##Statistic=selection sum;mean;min;max;median
##Raster_output=output raster
##Raster_output_count=output raster

from qgis.core import QgsVectorLayer

if Statistic == 0:
    stat = "sum"
elif Statistic == 1:
    stat = "mean"
elif Statistic == 2:
    stat = "min"
elif Statistic == 3:
    stat = "max"
elif Statistic == 4:
    stat = "median"

#get epsg from point layer 
PointLayer = QgsVectorLayer(Point_layer, "layer_name", "ogr")
epsg = PointLayer.crs().authid()

#generate vector grid with point stats
grid = processing.runalg("qgis:creategrid",1,Extent,Resolution,Resolution,epsg,None)
joined = processing.runalg("qgis:joinattributesbylocation", grid['OUTPUT'],PointLayer,['intersects'],0.0,1,stat,0,None)

#get extent from vector grid
joinedLayer = QgsVectorLayer(joined['OUTPUT'], "layer_name", "ogr")

extent2 = joinedLayer.extent()
xmin = extent2.xMinimum()
xmax = extent2.xMaximum()
ymin = extent2.yMinimum()
ymax = extent2.yMaximum()

#get created attribute (joined i.e."sumid")
attr = Attribute_values#chosen attribute
attr = str(stat)+str(attr)#"sum"+"id"

#generate stats raster
try:
    processing.runalg("gdalogr:rasterize", joined['OUTPUT'],attr,1,Resolution,Resolution,"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),False,5,"",4,75.0,6.0,1.0,False,0,"",Raster_output)
except:
    pass
    progress.setInfo("gdalogr:rasterize for point statistics fialed, running GRASS v.to.rast")
    processing.runalg("grass7:v.to.rast.attribute",joined['OUTPUT'],0,attr,"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),0,-1,0.0001,Raster_output)

#generate countpoints raster
countpoints = processing.runalg("qgis:countpointsinpolygon",joined['OUTPUT'],Point_layer,"NUMPOINTS",None)

try:
    processing.runalg("gdalogr:rasterize", countpoints['OUTPUT'],"NUMPOINTS",1,Resolution,Resolution,"%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),False,5,"",4,75.0,6.0,1.0,False,0,"",Raster_output_count)
except:
    pass
    progress.setInfo("gdalogr:rasterize for counted points fialed, running GRASS v.to.rast")
    processing.runalg("grass7:v.to.rast.attribute",countpoints['OUTPUT'],0,"NUMPOINTS","%f,%f,%f,%f"% (xmin, xmax, ymin, ymax),0,-1,0.0001,Raster_output_count)

Obr. 6 Vzhled výsledného nástroje

Obr. 8 Rastr součtu číselného atributu bodové vrstvy, vygenerovaný pomocí vytvořeného nástroje

Post navigation