QGIS - nástroje zpracování, modelář, PyQGIS

V této části se zaměříme na možnosti automatizace výpočetních kroků v prostředí QGIS https://qgis.org. Kontrétně:

• grafický modelář
• Python API (PyQGIS)

Příklad výpočetní úlohy

Nejprve si vyzkoušíme řešit následující výpočetní úlohu pomocí nástrojů zpracování (Processing > Toolbox).

Zadání

Vyberte body z datasetu LUCAS nacházející se v jižní Itálii v oblastech, kde může dojít k zamrzání vody (nadmořká výška nad 1000 m n. m.). Výsledek validujte na základě CORINE Land Cover (CLC).

Vstupní data

1. Stáhněte DTM pokrývající jižní Itálii ve formátu GeoTIFF
   • Dlaždice E40N10
   • https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1
2. Stáhněte CLC pro Evropu ve formátu GeoPackage
   • https://land.copernicus.eu/pan-european/corine-land-cover/clc2018?tab=mapview
3. Stáhněte data LUCAS z projektu ST_LUCAS pro danou zájmovou oblast
   1. Nainstalujte zásuvný modul "ST_LUCAS Download Manager"
      • Plugins > Manage and Install Plugins... > LUCAS (Search)
   2. Pomocí zásuvného modulu stáhněte LUCAS data pro dané zájmové území a pro rok 2018 (pro tento rok již máme stažena data CLC)

Postup výpočtu

1. Identifikujte území s nadmořskou výškou nad 1000 m n.m.
   1. Vytvořte masku pro území s nadmořskou výškou nad 1000 m n.m. Tuto úlohu můžeme řešit pomocí rastrového kalkulátoru (Raster > Raster Calculator...)
      • if(dem > 1000, 1, 0/0)
   2. Rastrovou masku zvektorizujeme
      • Raster > Conversion > Polygonize (Raster to Vector)...
2. Omezte body LUCAS na základě vektorové masky
   • Vector > Geoprocessing Tools > Clip...
3. Vyberte body LUCAS s třídou krajinného pokryvu korespondující s vodou (viz atribut LC1)
   • Vector > Select Features using an Expression
   • Podmínka: "lc1" LIKE 'G%'
4. Výběr uložte do nové vrstvy
   • Pravé tlačítko nad vrstvou > Export > Save selected features as...
5. Nikdy nevěřte jednomu datovému zdroji - proveďte validaci na základě vodních ploch z datasetu CLC
   1. Ořežte plochy CLC na základě vektorové masky zájmového území
      • Vector > Geoprocessing Tools > Clip...
   2. Vyberte plochy s třídou korespondující s vodou (viz atribut Code_18)
      • Select Features using an Expression
      • "Code_18" > 500
   3. Výběr uložte do nové vrstvy
      • Právé tlačítko na vrstvou > Export > Save selected features as...
   4. Vybrané CLC plochy rozšiřte pomocí obalové zóny o 100m
      • Vector > Geoprocessing Tools > Buffer...
   5. Vyberte body LUCAS ležicí v rozšířených CLC plochách
      • Select by Location

Výsledek úlohy je následující:

Grafický modelář

Pro výše uvedenou úlohu vytvořte model: Processing > Graphical Modeler.

Postup:

• do modelu nejprve přidejte vstupní data (záložka Inputs)
• posléze přidávejte výpočetní nástroje (záložka Algorithms)

Výsledný model může vypadat následovně (ke stažení zde):

Python API

Processing

Stáhněte data LUCAS pro vaše zájmové území a proveďte následující výpočet pomocí Python modulu processing:

1. vyberte body LUCAS na první úrovni nomenklatury Land Cover (LC) , viz řádek 4. K tomu použijte nástroj native:extractbyexpression, viz řádek 13.
2. ořežte plochy CLC na základě prostorového rozsahu vybraných bodů LUCAS gdal:clipvectorbyextent, viz řádek 21.
3. dále omezte výběr ploch CLC na první úrovni nomenklatury LC podobně jako u bodu 1, viz řádek 6 a 29
4. vytvořte obalovou zónu kolem vybraných CLC ploch native:buffer, viz řádek 7 a 37
5. na závěr vyberte body LUCAS, které leží v rozšířených CLC plochách native:extarctbylocation, viz řádek 45

 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

import os

lucas_input = os.path.join(os.environ["HOME"], "Downloads", "lucas.gpkg")
lucas_class_l1 = "G"
clc_input = os.path.join(os.environ["HOME"], "Downloads", "02", "U2018_CLC2018_V2020_20u1.gpkg")
clc_class_l1 = "5"
clc_buffer = 100

output_dir = os.path.join(os.environ["HOME"], "Downloads")

print("selecting LUCAS points...")
lucas_selected = os.path.join(output_dir, "lucas_selected.gpkg")
processing.run('native:extractbyexpression', {
            'INPUT': lucas_input,
            'EXPRESSION': f'"lc1" LIKE \'{lucas_class_l1}%\'',
            'OUTPUT': lucas_selected
})

print("clipping CLC by selected LUCAS points...")
clc_clipped = os.path.join(output_dir, "clc_clipped.gpkg")
processing.run('gdal:clipvectorbyextent', {
            'INPUT': clc_input,
            'EXTENT': lucas_selected,
            'OUTPUT': clc_clipped,
})

print("selecting CLC areas...")
clc_selected = os.path.join(output_dir, "clc_selected.gpkg")
processing.run('native:extractbyexpression', {
            'INPUT': clc_clipped,
            'EXPRESSION': f'"Code_18" > {clc_class_l1}00',
            'OUTPUT': clc_selected
})    

print("creating CLC buffers...")
clc_selected_buffer = os.path.join(output_dir, "clc_selected_buffer.gpkg")
processing.run('native:buffer', {
            'INPUT': clc_selected,
            'DISTANCE': clc_buffer,
            'OUTPUT': clc_selected_buffer
        })

print("selecting LUCAS points based on CLC buffers...")
result = os.path.join(output_dir, "luca_clc_selected_buffer.gpkg")
processing.run('native:extractbylocation', {
            'INPUT': lucas_selected,
            'INTERSECT': clc_selected_buffer,
            'PREDICATE': [0],  # intersect
            'OUTPUT': result
        })

print("done")

PyQGIS

Dokumentace:

• PyQGIS Developer Cookbook
• QGIS Python API

Projděte všechny vrstvy a vypište jejich souřadnicový systém.

• QgsMapLayer
• QgsCoordinateReferenceSystem

from qgis.core import QgsProject

# mapLayers() vraci slovnik
for layer in QgsProject.instance().mapLayers().values():
    print(f'{layer.name()}: {layer.crs().authid()}')

Projděte všechny body ve vybrané vrstvě a vypište jejich souřadnice a hodnotu atributu lc1.

• QgsGeometry

from qgis.core import QgsProject

# mapLayersByName() vraci seznam
for feature in QgsProject.instance().mapLayersByName('lucas')[0].getFeatures():
    # souradnice lze vypsat i pomoci feature.geometry(), ale pokud chceme
    # jejich hodnoty, musime z nich vytvorit objekt typu Point
    # print(point.geometry())

    point = feature.geometry().asPoint()
    print(f'[{point.x()}, {point.y()}]: {feature["lc1"]}')

Projděte všechny body ve vybrané vrstvě, vytvořte kolem nich obalovou zónu s poloměrem určeným atributem obs_dist, a vypočtěte její obsah.

from qgis.core import QgsProject

# mapLayersByName() vraci seznam
for feature in QgsProject.instance().mapLayersByName('lucas')[0].getFeatures():
    if feature['obs_dist']:  # nutno pokud jsou nektere prvky NULL
        # druhym parametrem je pocet lomovych bodu pouzitych k vytvoreni kruznice
        # (jedna se o body navic ke ctyrem potrebnym pro ctverec)
        buffer = feature.geometry().buffer(feature["obs_dist"], 5)

        point = feature.geometry().asPoint()
        print(f'[{point.x()}, {point.y()}]: {buffer.area()} m2')