Próbkowanie

Próbkowanie (sampling) to proces wybierania (selekcji) podzbioru lokalizacji danych przestrzennych z większego obszaru geograficznego w celu analizy cech całego obszaru. Stworzona reprezentatywna próba danych, może zostać wykorzystana do wyciągnięcia wniosków na temat populacji, przy jednoczesnym zminimalizowaniu czasu, kosztów i wysiłku wymaganego do zebrania i przetworzenia tych danych.

Rodzaje próbkowania:

losowe – lokalizacje wybierane są całkowicie przypadkowo.
stratyfikowane (łac. stratum – warstwa) – obszar dzielony jest na podgrupy (stratum) na podstawie pewnych cech (np. użytkowanie terenu, wysokość). Następnie pobierane są losowe próbki w obrębie każdego stratum, co zapewnia zbalansowaną reprezentację różnych podgrup.
regularne / systematyczne – próbki wybierane są w określonych i równomiernych odstępach, np. co 100 m.
celowe – próbki wybierane są na podstawie wiedzy badacza i określonych kryteriów.

Próbkowanie losowe

Istnieje wiele sposobów losowego doboru próbek używając różnych rozkładów prawdopodobieństwa (np. jednostajnego, normalnego, dwumianowego czy Poissona) w zależności od charakterystyki danych i konkretnych wymagań analizy. W przypadku danych przestrzennych wykorzystywane są specjalistyczne rozkłady przestrzenne, które są bardziej złożone od tradycyjnych ze względu na uwzględnienie relacji przestrzennych (patrz pakiet spatstat.random w R).

W najprostszym przypadku należy wykorzystać losowanie ze zwracaniem, gdzie jako obiekt wejściowy definiujemy listę, która reprezentuje minimalny oraz maksymalny zasięg przestrzenny danej osi, np. x_coords = list(range(xMinimum, xMaximum + 1)). Z racji, iż dane rastrowe są reprezentacją dyskretną obszarów badań, trzeba zdefiniować interwał próbkowania, np. co 10 m, aby uniknąć sytuacji, w której próbki pochodzą z identycznej lokalizacji (komórki). Zakres przestrzenny możemy zdefiniować na dwa sposoby – dla całego rastra lub dla określonego obszaru. Jeśli korzystamy z QGIS, to dodatkowo możemy, pobrać aktualnie wyświetlany zasięg warstwy w aplikacji.

Uwaga! Generowanie liczb losowych jest procesem losowym, to znaczy że za każdym razem, gdy kod zostanie wykonany, to otrzymamy inną sekwencję liczb. Aby uzyskać tę samą sekwencję liczb, należy ustawić ziarno losowości (random.seed()) najlepiej na początku skryptu, dzięki czemu wyniki będą powtarzalne.

import random
random.seed(999) # ziarno losowości

x_coords = list(range(0, 100 + 1, 10))
# losowanie ze zwracaniem
x_coords = random.choices(x_coords, k = 5)
print(x_coords)

[80, 0, 90, 60, 50]

Do stworzenia punktu wymagana jest para współrzędnych X (długość geograficzna) i Y (szerokość geograficzna) oraz klasa QgsPointXY. Wygenerowane punkty możemy przechowywać w liście (co jest prostszym rozwiązaniem) lub stworzyć nową warstwę wektorową (patrz temat “Przetwarzanie danych wektorowych: Tworzenie nowej warstwy”).

from qgis.core import QgsPointXY

x_coords = list(range(0, 100 + 1, 10))
x_coords = random.choices(x_coords, k = 5)
y_coords = list(range(0, 100 + 1, 10))
y_coords = random.choices(y_coords, k = 5)

points = []
for x, y in zip(x_coords, y_coords):
    points.append(QgsPointXY(x, y))

print(points)

[<QgsPointXY: POINT(10 10)>, <QgsPointXY: POINT(80 100)>, <QgsPointXY: POINT(30 20)>, <QgsPointXY: POINT(80 20)>, <QgsPointXY: POINT(90 60)>]

Próbkowanie stratyfikowane

Procedura wygenerowania losowego punkty w poligonie składa się z dwóch etapów:

Wylosowanie punktu, który znajduje się w zakresie przestrzennym (bounding box) poligonu.
Sprawdzenie czy wylosowany punkt znajduje się w poligonie (lub czy poligon zawiera ten punkt).

from qgis.core import QgsGeometry

# przykładowe dane (trójkąt)
wkt = "POLYGON((0 0, 20 0, 5 10, 0 0))"
polygon = QgsGeometry.fromWkt(wkt)
bbox = polygon.boundingBox()
print(bbox) # xmin, ymin, xmax, ymax

<QgsRectangle: 0 0, 20 10>

# (1) losowanie punktu z zakresu
x = random.uniform(bbox.xMinimum(), bbox.xMaximum())
y = random.uniform(bbox.yMinimum(), bbox.yMaximum())
pt = QgsPointXY(x, y)
print(pt)

<QgsPointXY: POINT(14.04204171589216621 2.53705020234383216)>

# (2) sprawdzenie czy poligon zawiera punkt
print(polygon.contains(pt))

True

Naturalnie, jeśli chcemy wygenerować punkty dla obiektów, to musimy zastosować pętlę. W przypadku “dziwnych” geometrii, które trudno opróbkować, należy pamiętać żeby pętla była skończona.

Próbkowanie regularne

W poniższym przykładzie stworzymy linię na podstawie dwóch punktów, a następnie wygenerujemy punkty w stałych odstępach wzdłuż tej linii. Kluczowe jest zastosowanie metody interpolate(), która oblicza dokładne położenie punktu na linii w określonej odległości od punktu początkowego.

from qgis.core import QgsLineString

pt1 = QgsPointXY(10, 20)
pt2 = QgsPointXY(50, 70)
line = QgsGeometry(QgsLineString([pt1, pt2])) # QgsGeometry.fromPolylineXY()
print(line)

line_length = line.length()
print("Długość linii:", line_length)

<QgsGeometry: LineString (10 20, 50 70)>
Długość linii: 64.03124237432849

interval = 10
points = []
distance = 0
while distance <= line_length:
    point = line.interpolate(distance)
    points.append(point.asPoint())
    distance = distance + interval
print("Liczba punktów:", len(points))

Liczba punktów: 7

W przypadku próbkowania wzdłuż linii istotne jest sprawdzenie kierunku przebiegu linii, np. czy jest to kierunek wschód-zachód czy na odwrót. Jeśli błędnie przyjmiemy kierunek, to wartości na krzywej będą odwrócone. Do weryfikacji można napisać prostą funkcję, która porówna współrzędne punktu początkowego i końcowego linii.

def line_direction(line):
    # pobierz wierzchołki
    vertices = line.asPolyline()
    
    # wybierz pierwszy i ostatni wierzchołek
    start_point = vertices[0]
    end_point = vertices[-1]
    
    # sprawdzenie osi X
    if start_point.x() < end_point.x():
        x_direction = "ZACH-WSCH"
    elif start_point.x() > end_point.x():
        x_direction = "WSCH-ZACH"
    else:
        x_direction = "" # kierunek stały
        
    # sprawdzenie osi Y
    if start_point.y() < end_point.y():
        y_direction = "PŁD-PŁN"
    elif start_point.y() > end_point.y():
        y_direction = "PŁN-PŁD"
    else:
        y_direction = ""
        
    return x_direction, y_direction

print("Kierunek linii:", line_direction(line))

Kierunek linii: ('ZACH-WSCH', 'PŁD-PŁN')

Alternatywnie również można obliczyć azymut pomiędzy punktem początkowym i końcowym.

Zadania

Napisz funkcję sample_random(raster, n), która posłuży do losowego próbkowania rastra wielokanałowego. Jako wynik zostanie zwrócona ramka danych (biblioteka pandas) składającą się ze współrzędnych X i Y punktów oraz wartości kanałów rastra. Nadaj kolumnom odpowiednie nazwy. Jeśli komórka rastra nie posiada wartości, to zwróć nan w tabeli.
Napisz funkcję sample_regular(layer, interval), która stworzy regularną siatkę punktów dla zakresu przestrzennego podanej warstwy wektorowej bądź rastrowej i określonego interwału (np. co 100 m). Wykorzystaj pętlę while (zamiast for), aby zapewnić wsparcie dla liczb zmiennoprzecinkowych. Jako wynik zostanie zwrócona warstwa punktowa w pamięci (warstwa tymczasowa).
Napisz funkcję sample_strata(vector, n, save_path) służącą do losowania stratyfikowanego, która zagwarantuje, że dla każdego poligonu zostanie zwrócona dokładna liczba punktów wskazana przez użytkownika, a wynik zostanie zapisany jako warstwa wektorowa. Do analizy wykorzystaj dane z pliku powiaty.gpkg.