Algorytmy danych geoprzestrzennych
Próbkowanie

Krzysztof Dyba

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:

  1. Wylosowanie punktu, który znajduje się z zakresie (bounding box) poligonu.
  2. Sprawdzenie czy wylosowany punkt znajduje się w poligonie (lub czy poligon zawiera ten punkt).
from qgis.core import QgsGeometry

# przykładowe dane
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

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')

Zadania

  1. Analiza geomorfometryczna. Dokonaj analizy ukształtowania terenu (plik DEM.tif) używając dwóch transektów północ-południe oraz wschód-zachód. Wynik zaprezentuj na rycinie (na osi X powinna znajdować się odległość wyrażona w metrach). Zastanów się również, co zrobić w przypadku, gdy komórki nie posiadają wartości (nan).

  2. Analiza teledetekcyjna. Wygeneruj losową próbę i następnie dla dwóch wybranych kanałów spektralnych z Landsat 8/9:

  • wykonaj wykres rozrzutu uwzględniając przezroczystość punktów oraz zaznacz linię x = y,

  • oblicz statystyki opisowe z wartości różnic między tymi kanałami.

    Dodatkowo można obliczyć współczynnik korelacji Pearsona (numpy.corrcoef()).

  1. Analiza geoinformacyjna. Zaimplementuj funkcję losowania stratyfikowanego sample_strata(wektor, n), 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.