Algorytmy danych geoprzestrzennych
Wprowadzenie do PyQGIS

Krzysztof Dyba

PyQGIS to biblioteka języka Python udostępniająca interfejs programistyczny do aplikacji QGIS, dzięki czemu możliwa jest interakcja pomiędzy kodem napisanym w Pythonie, a środowiskiem i elementami QGIS (np. projekt, warstwy, funkcje, itd.). Biblioteka odgrywa kluczową rolę w automatyzacji zadań oraz tworzeniu niestandardowych narzędzi i wtyczek do rozszerzenia funkcjonalności QGIS. PyQGIS jest częścią projektu QGIS, który jest darmowy i otwarty.

Przykładowe zastosowania:

Automatyzacja (przetwarzanie wsadowe, skrypty do czyszczenia i walidacji danych, generowanie map),
Tworzenie nowych narzędzi (wtyczki, specjalistyczne algorytmy geoprzetwarzania, niestandardowe interfejsy),
Zaawansowane analizy przestrzenne (wykorzystanie zewnętrznych bibliotek, np. numpy, scikit-learn).

Instalacja

Instalacja jest bardzo prosta – należy pobrać program QGIS z oficjalnej strony https://qgis.org/download/ i zainstalować według wyświetlanych instrukcji. Środowisko Python i podstawowe biblioteki są automatycznie instalowane wraz z QGIS.

Pierwsze kroki

Zasadniczo bibliotekę PyQGIS możemy wykorzystać na dwa sposoby:

Zintegrowana konsola Pythona w QGIS.

QGIS ma wbudowaną konsolę Pythona, w której można wykonywać skrypty interaktywnie. Aby ją włączyć kliknij Wtyczki > Konsola Pythona w pasku narzędzi. Można również wykorzystać skrót klawiszowy: CTRL + ALT + P. Oprócz konsoli, wbudowany jest bardzo prosty edytor kodu.

Samodzielne skrypty.

PyQGIS można również wykorzystać poza aplikacją QGIS w samodzielnych skryptach Pythona.

Sprawdź czy wszystko działa prawidłowo wykonując następujące polecenie:

from qgis.core import *

Dane wektorowe

Wczytywanie

W folderze dane znajduje się plik powiaty.gpkg, który zawiera granice powiatów w wielkopolsce z Państwowego Rejestru Granic. Żeby go wczytać należy wykorzystać klasę QgsVectorLayer i zdefiniować trzy argumenty:

ścieżkę do pliku,
nazwę warstwy pod jaką zostanie wyświetlony,
backend do wczytania danych (data provider), np. ogr, spatialite, postgres.

Ścieżkę do pliku można zdefiniować na dwa sposoby. Pierwszy to podanie ścieżki bezwzględnej, tj. wskazanie dokładnej lokalizacji, w której znajduje się plik. Na przykład:

sciezka = "C:\Users\Krzysztof\Desktop\dane.csv"
# lub sciezka = "C:\\Users\\Krzysztof\\Desktop\\dane.csv"

Nie jest to jednak rekomendowana metoda, ponieważ zakłada, że struktura katalogów pomiędzy różnymi komputerami i systemami jest identyczna. Drugi sposób polega na podaniu ścieżki względnej. W tym przypadku podajemy lokalizację pliku względem bieżącego katalogu roboczego lub projektu. Aby dowiedzieć się, gdzie znajduje się katalog roboczy, możemy użyć metody getcwd() z biblioteki os, a do jego zmiany metody chdir(). Uwaga: domyślna ścieżka katalogu roboczego w konsoli QGIS, a systemowej są różne!

import os
print(os.getcwd())

C:\Users\Krzysztof\Desktop

Zauważ również, że w zależności od systemu jest wykorzystywany różny separator katalogów w hierarchii. Aby prawidłowo wskazać ścieżkę do pliku należy wykorzystać metodę os.path.join().

print(os.path.join("projekt", "dane", "pomiary.csv"))

projekt\dane\pomiary.csv

Spróbujmy teraz wczytać nasz plik.

from qgis.core import QgsVectorLayer

# określenie ścieżki do pliku
sciezka_do_pliku = os.path.join("algorytmy-geoprzestrzenne", "dane", "powiaty.gpkg")

# wczytanie warstwy wektorowej
wektor = QgsVectorLayer(sciezka_do_pliku, "powiaty", "ogr")

# sprawdzenie czy warstwa została załadowana prawidłowo
print(wektor.isValid())

True

Warstwa została wczytana prawidłowo. Następnie możemy wyświetlić ją w QGIS (pod warunkiem, że kod wykonujemy w QGIS).

from qgis.core import QgsProject
QgsProject.instance().addMapLayer(wektor)

Wyświetlanie metadanych

W kolejnym kroku możemy wyświetlić podstawowe informacje o warstwie, tj.:

rodzaj geometrii,
zakres przestrzenny,
liczba obiektów oraz liczba atrybutów,
układ współrzędnych.

# rodzaj geometrii
rodzaj_geometrii = wektor.geometryType()
rodzaj_geometrii = rodzaj_geometrii.name # zwraca nazwę
print("Rodzaj geometrii:", rodzaj_geometrii)

# zakres przestrzenny
zakres = wektor.extent() # obiekt QgsRectangle
zakres = zakres.toString() # konwersja na tekst
print("Zakres przestrzenny:", zakres)

# liczba obiektów
liczba_obiekty = wektor.featureCount()
print("Liczba obiektów:", liczba_obiekty)

# liczba atrybutów
liczba_atrybuty = wektor.fields()
liczba_atrybuty = liczba_atrybuty.count()
print("Liczba atrybutów:", liczba_atrybuty)

# układ współrzędnych
uklad = wektor.crs()
uklad = uklad.authid() # zwraca kod EPSG
print("CRS:", uklad)

Rodzaj geometrii: Polygon
Zakres przestrzenny: 281949.3437500000000000,360241.7500000000000000 : 507167.3125000000000000,645513.6250000000000000
Liczba obiektów: 35
Liczba atrybutów: 6
CRS: EPSG:2180

Zapisywanie

Do zapisu danych wektorowych można wykorzystać dwie klasy:

QgsVectorFileWriter – zapis danych wyłącznie jako pliki na dysku.
QgsVectorLayerExporter – możliwe wykorzystanie różnych backendów (data providers), np. zapis do bazy danych.

Pierwszy sposób używając QgsVectorFileWriter.

from qgis.core import QgsVectorFileWriter, QgsCoordinateTransformContext

output_path = "test.geojson"
# opcje zapisu
options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "GeoJSON"
options.fileEncoding = "UTF-8"

writer = QgsVectorFileWriter.writeAsVectorFormatV3(
    layer = wektor,
    fileName = output_path,
    transformContext = QgsCoordinateTransformContext(),
    options = options
)

if writer[0] != 0:
    print("Błąd zapisu")
else: 
    print("OK")

Alternatywne podejście z zastosowaniem QgsVectorLayerExporter.

from qgis.core import QgsVectorLayerExporter

output_path = "test.geojson"
output_format = "GeoJSON"

QgsVectorLayerExporter.exportLayer(
    layer = wektor,
    uri = output_path,
    providerKey = "ogr",
    destCRS = wektor.crs(),
    options = {"driverName": output_format} # jako słownik
)

if writer[0] != 0:
    print("Błąd zapisu")
else: 
    print("OK")

Dane rastrowe

Wczytywanie

Wczytanie danych rastrowych wygląda podobnie jak w przypadku danych wektorowych z tą różnicą, iż służy do tego klasa QgsRasterLayer. Wczytajmy plik DEM.tif, który znajduje się w katalogu dane.

from qgis.core import QgsRasterLayer

sciezka_do_pliku = os.path.join("algorytmy-geoprzestrzenne", "dane", "DEM.tif")
raster = QgsRasterLayer(sciezka_do_pliku, "DEM")
print(raster.isValid())

True

Po prawidłowym wczytaniu warstwy, również możemy ją zwizualizować w QGIS.

from qgis.core import QgsProject
QgsProject.instance().addMapLayer(raster)

Wyświetlanie metadanych

Używając odpowiednich metod na rastrze możemy uzyskać następujące informacje:

liczba wierszy, kolumn i kanałów,
zakres przestrzenny,
układ przestrzenny,
rozdzielczość komórki.

print("Liczba kolumn (Szerokość):", raster.width())
print("Liczba wierszy (Wysokość):", raster.height())
print("Liczba kanałów:", raster.bandCount())
print("Zakres:", raster.extent().toString())
print("CRS:", raster.crs().authid())
print("Rozdzielczość (Rozmiar komórki):", raster.rasterUnitsPerPixelX(),
      raster.rasterUnitsPerPixelY())

Liczba kolumn (Szerokość): 533
Liczba wierszy (Wysokość): 608
Liczba kanałów: 1
Zakres: 253698.3311999999859836,353734.3616000000038184 : 520058.2303000000538304,657570.7552000000141561
CRS: EPSG:2180
Rozdzielczość (Rozmiar komórki): 499.73714652908075 499.7309105263158

Co więcej, każdy z kanałów rastra zawiera dodatkowe informacje, takie jak:

nazwa kanału,
wartość reprezentująca brak danych (NoData),
typ danych,
wartość maksymalna,
wartość minimalna.

# nasz raster składa się tylko z jednego kanału
for band in range(1, raster.bandCount() + 1):
    print("Nazwa kanału:", raster.bandName(band))
    print("Wartość NoData:", raster.dataProvider().sourceNoDataValue(band))
    print("Typ danych:", raster.dataProvider().dataType(band).name)
    print("Wartość minimalna:", raster.dataProvider().bandStatistics(band).minimumValue)
    print("Wartość maksymalna:", raster.dataProvider().bandStatistics(band).maximumValue)

Nazwa kanału: Band 1
Wartość NoData: 9999.0
Typ danych: Float32
Wartość minimalna: -68.71467590332031
Wartość maksymalna: 459.58660888671875

Zapisywanie

Do zapisu danych rastrowych służy klasa QgsRasterFileWriter.

from qgis.core import QgsRasterFileWriter, QgsRasterPipe, QgsCoordinateTransformContext

output_path = "test.tif"

writer = QgsRasterFileWriter(output_path)
writer.setCreateOptions(["COMPRESS=LZW"]) # opcje jako lista
pipe = QgsRasterPipe()
provider = raster.dataProvider()
pipe.set(provider.clone())

status = writer.writeRaster(
    pipe = pipe,
    nCols = provider.xSize(),
    nRows = provider.ySize(),
    outputExtent = provider.extent(),
    crs = provider.crs(),
    transformContext = QgsCoordinateTransformContext()
)

# writeRaster zwraca RasterFileWriterResult
if status != 0:
    print("Błąd zapisu")
else: 
    print("OK")

Zadania

Napisz funkcję, która obliczy wielkość komórki na podstawie zakresu przestrzennego oraz liczby kolumn i wierszy rastra.
Napisz funkcję (wyswietl_metadane), która zaprezentuje metadane warstwy wektorowej oraz rastrowej w czytelny sposób. Do sprawdzenia typu obiektu wykorzystaj funkcję isinstance.