Algorytmy danych geoprzestrzennych
Przetwarzanie danych wektorowych

Krzysztof Dyba

Na początku wczytajmy warstwę wektorową powiaty.gpkg, a następnie sprawdźmy podstawowe metadane.

import os
from qgis.core import QgsVectorLayer

# wczytanie warstwy wektorowej
filepath = os.path.join("dane", "powiaty.gpkg")
vector = QgsVectorLayer(filepath, "powiaty", "ogr")
print(vector.isValid())

True

# wyświetlenie metadanych
print(vector.featureCount())
print(vector.fields().count())
print(vector.fields().names())

35
5
['fid', 'JPT_KOD_JE', 'JPT_NAZWA_', 'IIP_IDENTY', 'IIP_WERSJA']

Wczytana warstwa składa się z 35 obiektów oraz 5 atrybutów (inne określenie to kolumny bądź pola).

Dostęp do obiektów

Iterowanie po obiektach warstwy wektorowej umożliwia dostęp lub modyfikowanie geometrii/atrybutów tych obiektów. Do tego celu służy metoda getFeatures(). Geometria obiektu dostępna jest za pomocą metody geometry() (klasa QgsGeometry), natomiast atrybuty za pomocą metody attributes().

W poniższym przykładzie wykorzystamy funkcję islice z biblioteki itertools, która zwróci tylko 5 pierwszych obiektów.

from itertools import islice

# iteracja po obiektach
for feature in islice(vector.getFeatures(), 5):
    # dostęp do atrybutów obiektu
    attrs = feature.attributes()
    # dostęp do geometrii obiektu
    geom = feature.geometry()
    
    print(attrs)
    # print(geom)

[1, '3012', 'powiat krotoszyński', 'dc54ce6a-5272-45fc-a49c-824dbf170f2c', '2024-09-05T13:41:42+02:00']
[2, '3002', 'powiat czarnkowsko-trzcianecki', 'b2e00ebb-412b-4ae8-9afc-fe68f65bb1ac', '2021-05-11T10:49:15+02:00']
[3, '3009', 'powiat kolski', '06a4a975-dc89-4ac3-9d51-cf188532908a', '2012-09-27T08:59:03+02:00']
[4, '3029', 'powiat wolsztyński', '6f988fb3-2b2c-467e-9d07-c512d9215c05', '2024-08-27T13:56:13+02:00']
[5, '3026', 'powiat śremski', '2fc58bb6-d095-4f37-8f10-874797cb97fb', '2024-09-03T13:42:56+02:00']

Obliczanie powierzchni

W celu obliczenia powierzchni poligonów musimy wykonać pętle po wszystkich obiektach znajdujących się w warstwie i uzyskać dostęp do ich geometrii. Następnie, należy wykorzystać metodę area(), która oblicza pole powierzchni w jednostkach układu współrzędnych warstwy (zazwyczaj są to metry kwadratowe) w układzie planarnym (kartezjańskim). Opcjonalnie, możemy dokonać konwersji jednostek, np. na kilometry kwadratowe czy hektary. W niniejszym przykładzie wszystkie obiekty w warstwie posiadają geometrie.

for feature in islice(vector.getFeatures(), 5):
    area = feature.geometry().area()
    area = area / 1000**2 # konwersja na km^2
    print(area)

712.7437283077061
1806.6288683868954
1009.273272845488
679.340479768619
573.5639814826102

Dodawanie atrybutów

Dodawanie atrybutów do warstwy wektorowej obejmuje modyfikację tabeli atrybutów poprzez dodanie nowych atrybutów i opcjonalnie ich wypełnienie. Operacja składa się z kilku kroków:

Na początku należy włączyć tryb edycji warstwy startEditing(), aby ją zmodyfikować.
Kolejny krok to zdefiniowanie nowego atrybutu wykorzystując klasę QgsField. Każdy atrybut posiada nazwę, typ danych i opcjonalne ograniczenia (długość, precyzja). Następnie należy użyć metod addAttribute() oraz updateFields(), aby dodać nowe pole do tabeli atrybutów i ją zaktualizować. Zdefiniowanie kilku atrybutów wymaga przekazanie ich jako listy i dodanie w pętli.
Jeśli chcemy uzupełnić nowo utworzony atrybut, to musimy zastosować pętle po obiektach i zaktualizować wartość wybranego atrybutu używając metody changeAttributeValue(), w której kolejno określimy ID obiektu, indeks atrybutu oraz wartość.
W ostatnim kroku należy zatwierdzić i zapisać wprowadzone zmiany używając commitChanges().

Uwaga odnośnie typów danych atrybutów w QGIS! Od wersji 3.38 typy danych zdefiniowane są w klasie QMetaType (np. QMetaType.Type.Double, QMetaType.Type.QString). Wcześniej była to klasa QVariant (np. QVariant.Double, QVariant.String).

from qgis.core import QgsField
from qgis.PyQt.QtCore import QMetaType

# (1) rozpoczęcie trybu edycji warstwy
vector.startEditing()

# (2) dodanie nowego atrybutu do tabeli
new_field = QgsField("pole_km2", QMetaType.Type.Double)
vector.addAttribute(new_field)
vector.updateFields()

# indeks nowego atrybutu
area_idx = vector.fields().indexOf("pole_km2")

# (3) obliczenie powierzchni dla każdego obiektu
for feature in vector.getFeatures():
    area = feature.geometry().area()
    area = area / 1000**2
    # wprowadzenie wartości do atrybutu obiektu
    vector.changeAttributeValue(feature.id(), area_idx, area)

# (4) zapisanie zmian
vector.commitChanges()

True

Sprawdźmy teraz czy wykonany kod zadziałał prawidłowo.

print("pole_km2" in vector.fields().names())

for feature in islice(vector.getFeatures(), 5):
    print(feature.attribute("pole_km2"))

True
712.7437283077061
1806.6288683868954
1009.273272845488
679.340479768619
573.5639814826102

Wybieranie obiektów

Obiekty z warstwy wektorowej można wybrać (wyselekcjonować) na podstawie różnych kryteriów, takich jak wartości atrybutów czy położenie przestrzenne.

Atrybut

Zasadniczo, istnieją dwa sposoby wybierania obiektów z uwzględnieniem wartości atrybutów:

selectByExpression() z klasy QgsVectorLayer.
setFilterExpression() z klasy QgsFeatureRequest.

Pierwszy sposób selectByExpression() to metoda służąca do wybierania obiektów w warstwie na podstawie określonego wyrażenia. Co najważniejsze, wybór odbywa się z uwzględnieniem wszystkich obiektów i jest tymczasowy.

expression = '"pole_km2" > 1300'
vector.selectByExpression(expression)
selected_features = vector.selectedFeatures()
print(len(selected_features)) # liczba obiektów

# for feature in selected_features:
    # print(feature.attributes())

Jeśli kod wykonaliśmy w QGIS, to wybrane obiekty możemy podświetlić na mapie w następujący sposób:

from qgis.core import QgsProject
QgsProject.instance().addMapLayer(vector)
vector.selectByIds([feature.id() for feature in selected_features])

Natomiast drugi sposób, setFilterExpression() stosuje filtr przy wczytywaniu danych (getFeatures()). Obiekty, które nie spełniają kryteriów, nie są dodawane do sesji. Zwiększa to wydajność przetwarzania poprzez zmniejszenie ilości wczytywanych danych, co jest szczególnie przydatne w przypadku dużych zbiorów.

from qgis.core import QgsFeatureRequest

expression = '"pole_km2" > 1300'
# stworzenie żądania używając wyrażenia
request = QgsFeatureRequest().setFilterExpression(expression)

# for feature in vector.getFeatures(request):
    # print(feature.attributes())

Lokalizacja

Metoda selectByRect() umożliwia dokonanie prostej selekcji obiektów używając zakresu przestrzennego (bounding box) zdefiniowanego przez prostokąt (QgsRectangle) i cztery współrzędne (xmin, ymin, xmax, ymax). Bardziej zaawansowane zapytania przestrzenne pozwalają na wybór obiektów na podstawie ich relacji przestrzennej z innymi geometriami (np. wybór obiektów w obrębie poligonu, przecięcie linii).

from qgis.core import QgsRectangle

# wybranie obiektów używając prostokąta
rect = QgsRectangle(340000, 480000, 380000, 520000) # xmin, ymin, xmax, ymax
vector.selectByRect(rect)
selected_features = vector.selectedFeatures()
print(len(selected_features)) # liczba obiektów

# for feature in selected_features:
    # print(feature.attributes())

Tworzenie nowej warstwy

Do tworzenia nowej warstwy wektorowej stosuje się tę samą klasę, co do wczytywania, czyli QgsVectorLayer, z tą różnicą iż określa się:

Typ geometrii warstwy oraz układ odniesienia (CRS). QGIS obsługuje kilka typów geometrii, np. punkty (Point), linie (LineString), poligony (Polygon) i multigeometrie (Multi*). Do zdefiniowania układu odniesienia można wykorzystać kod EPSG, np. EPSG:4326. Format zapisu jest następujący: Point?crs=EPSG:4326.
Nazwę warstwy.
Backend (data provider). Nowe warstwy zazwyczaj są tworzone w pamięci (memory).

# typ geometrii i układ współrzędnych
geometry_type = "Polygon"
crs = "EPSG:2180"

# stworzenie nowej warstwy wektorowej
newlayer = QgsVectorLayer(geometry_type+"?crs="+crs, "Nowa warstwa", "memory")
print(newlayer.isValid())

True

Następnie możemy zdefiniować i dodać do warstwy atrybuty, tak jak zaprezentowano to w sekcji “Dodawanie atrybutów”. Należy pamiętać o włączonym trybie edycji warstwy, a na końcu o zatwierdzeniu zmian.

# rozpoczęcie trybu edycji warstwy
newlayer.startEditing()

# stworzenie atrybutów
fields = [
    QgsField("ID", QMetaType.Type.Int),          # ID obiektu
    QgsField("nazwa", QMetaType.Type.QString),   # nazwa obiektu
    QgsField("wartosc", QMetaType.Type.Double)   # wartość
]

# dodanie atrybutów do warstwy
for field in fields:
    newlayer.addAttribute(field)
newlayer.updateFields()

W kolejnym kroku możemy stworzyć obiekty (klasa QgsFeature), które składają się z geometrii oraz atrybutów. W poniższym przykładzie użyjemy wcześniej zdefiniowanego prostokąta rect i metody fromRect().

from qgis.core import QgsFeature, QgsGeometry

# stworzenie obiektów
feature = QgsFeature()
polygon_geometry = QgsGeometry.fromRect(rect)

Innym sposobem konstruowania geometrii obiektów używając współrzędnych jest pośrednie zastosowanie klasy QgsPointXY i odpowiednich metod, np. fromPointXY() (do stworzenia punktów), fromPolylineXY() (do stworzenia linii) czy fromPolygonXY() (do stworzenia poligonów) z klasy QgsGeometry.

Finalnie, nowemu obiektowi należy nadać geometrię (metoda setGeometry()) oraz przypisać atrybutom wartości (metoda setAttributes()), po czym wymagane jest dodanie obiektów do warstwy, aktualizacja zasięgu warstwy oraz zapisanie zmian.

# nadanie geometrii oraz wartości atrybutom
feature.setGeometry(polygon_geometry)
feature.setAttributes([1, "Prostokąt", 999])

# dodanie obiektów do warstwy
newlayer.addFeature(feature)
# aktualizacja zasięgu warstwy
newlayer.updateExtents()
# zapisanie zmian
newlayer.commitChanges()

True

Dokonajmy jeszcze weryfikacji stworzonej warstwy.

for feature in newlayer.getFeatures():
    print(feature.attributes())
    print(feature.geometry())

[1, 'Prostokąt', 999.0]
<QgsGeometry: Polygon ((340000 480000, 380000 480000, 380000 520000, 340000 520000, 340000 480000))>

# wyświetlenie w QGIS
# QgsProject.instance().addMapLayer(newlayer)

Łączenie geometrii obiektów

Powiaty reprezentowane jako osobne geometrie możemy połączyć w jeden obiekt reprezentujący całe województwo. W tym celu można wykorzystać metodę combine z klasy QgsGeometry pod warunkiem, iż wszystkie granice są styczne. Na początku musimy zainicjować pustą geometrię, a następnie łączymy ją z kolejnymi w pętli. Wynikiem operacji jest nowa geometria będąca połączeniem wszystkich geometrii.

# stworzenie nowej pustej geometrii
geometry = None

for feature in vector.getFeatures():
    geom = feature.geometry()
    if geometry is None:                  # pierwsza geometria jako początkowa
        geometry = QgsGeometry(geom)
    else:                                 # łączy kolejne geometrie
        geometry = geometry.combine(geom)
        
combined = QgsFeature()
combined.setGeometry(geometry)

Uwaga! Jeśli granice poligonów nie przylegają do siebie idealnie, to w wyniku tej operacji powstaną artefakty graficzne w postaci linii lub poligonów szczątkowych (sliver polygons).

Tworzenie buforów

Stworzenie warstwy z buforami wymaga zainicjowania nowej warstwy poligonowej w pamięci. Następnie iterujemy po wszystkich obiektach i tworzymy bufory używając metody buffer, w której trzeba określić odległość bufora (wartość może być ujemna w przypadku poligonów) oraz liczbę segmentów. Po wygenerowaniu wszystkich obiektów, należy dodać je do warstwy.

# stworzenie pustej warstwy w pamięci
buffer_layer = QgsVectorLayer("Polygon?crs=" + vector.crs().authid(), "Bufor", "memory")

# odległość buforu w jednostkach warstwy (m)
buffer_distance = 10000.0 # 10 km

features = []
for feature in [combined]: # traktujemy jeden obiekt jako listę
    geom = feature.geometry()
    # stworzenie geometrii bufora
    buffer_geom = geom.buffer(buffer_distance, segments = 30)
    buffer_feature = QgsFeature()
    buffer_feature.setGeometry(buffer_geom)
    features.append(buffer_feature)

# dodanie buforów do warstwy wektorowej
buffer_layer.dataProvider().addFeatures(features)
buffer_layer.updateExtents()

# wyświetlenie w QGIS
# QgsProject.instance().addMapLayer(buffer_layer)

Reprojekcja warstwy

Istnieje kilka różnych sposobów reprojekcji warstwy wektorowej:

Iteracja po obiektach i zastosowanie metody transform z klasy QgsCoordinateTransform.
Zapisanie warstwy z nowym układem współrzędnych używając QgsVectorFileWriter.
Wykorzystanie narzędzia native:reprojectlayer.

Pierwsza metoda polega na transformacji wszystkich obiektów w pętli. Należy zwrócić szczególną uwagę na fakt, iż ta operacja modyfikuje bezpośrednio geometrie (można nieświadomie nadpisać oryginalne dane)! Kluczowym elementem jest zainicjowanie klasy QgsCoordinateTransform z trzema argumentami:

źródłowym układem współrzędnych,
docelowym układem współrzędnych,
kontekstem transformacji (związany z parametrami transformacji między układami).

from qgis.core import QgsProject, QgsCoordinateTransform, QgsCoordinateReferenceSystem

source_crs = vector.crs()
target_crs = QgsCoordinateReferenceSystem("EPSG:4326")
transform = QgsCoordinateTransform(source_crs, target_crs, QgsProject.instance())

newlayer = QgsVectorLayer("Polygon?crs=" + target_crs.authid(), "powiaty_4326", "memory")
newlayer.dataProvider().addAttributes(vector.fields())
newlayer.updateFields()

features = []
for feature in vector.getFeatures():
    new_feature = feature
    geom = feature.geometry()
    geom.transform(transform)
    new_feature.setGeometry(geom)
    features.append(new_feature)
    
newlayer.dataProvider().addFeatures(features)
print(newlayer.crs().authid())

# wyświetlenie w QGIS
# QgsProject.instance().addMapLayer(newlayer)

EPSG:4326

Druga metoda jest zdecydowanie łatwiejsza, ponieważ wymaga wyłącznie zdefiniowania wymaganych argumentów przez klasę QgsCoordinateTransform i zapisanie nowej warstwy.

from qgis.core import QgsVectorFileWriter

source_crs = vector.crs()
target_crs = QgsCoordinateReferenceSystem("EPSG:4326")
context = QgsProject.instance().transformContext()

output_path = "powiaty_4326.gpkg"
options = QgsVectorFileWriter.SaveVectorOptions()
options.driverName = "GPKG"
options.ct = QgsCoordinateTransform(source_crs, target_crs, context)

writer = QgsVectorFileWriter.writeAsVectorFormatV3(
    layer = vector,
    fileName = output_path,
    transformContext = context,
    options = options
)

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

Zadania

Oblicz długość granic i dodaj jako nowy atrybut do warstwy.
Napisz funkcję, która obliczy podstawowe statystyki opisowe dla wskazanej kolumny numerycznej. Jeśli kolumna nie istnieje w tabeli, to funkcja zwróci błąd. Jako przykład można wykorzystać atrybut z powierzchnią powiatów.
Napisz funkcję, która zwróci centroidy dla warstwy poligonowej i wynik zapisze na dysku w formacie .gpkg.