Geostatystyczna analiza danych

Próbkowanie preferencyjne

Author

Anna Dmowska 

# Wczytanie niezbędnych bibliotek
library(sf)         # Dane wektorowe
library(tmap)       # Tworzenie map
library(gstat)      # Geostatystyka (kriging, semiwariogramy)
library(terra)      # Dane rastrowe
library(tidyverse)  # Przetwarzanie danych i wykresy (ggplot2)
library(patchwork)  # Łączenie wykresów
library(dismo)

1 Próbkowanie preferencyjne

Próbkowanie preferencyjne (ang. preferential sampling) to sytuacja, w której punkty pomiarowe nie są rozmieszczone w sposób losowy lub regularny na całym badanym terenie, lecz są celowo zagęszczane w określonych obszarach.

W poniższym przykładzie znacznie gęstsze opróbowanie występuje wokół punktów o wysokich wartościach zmiennej Ba. Taki typ próbkowania ma wpływ na wartości statystyk opisowych. Jeśli (jak w poniższym przykładzie) zagęszczono pomiary wokół wysokich wartości zmiennej Ba, średnia z próby będzie znacznie wyższa niż rzeczywista średnia dla całego badanego obszaru.

W celu zniwelowania wpływu nierównomiernego rozmieszczenia punktów stosuje się metody rozgrupowia danych.

# --- KROK 1: Wczytanie i czyszczenie danych ---
# Wczytujemy dane z pliku zewnętrznego
df = read.csv("data/cw2_dane.csv")

# Selekcja kolumn: współrzędne (Xutm, Yutm) oraz zmienna Ba_pref
# na.omit() usuwa rekordy zawierające braki danych (NA)
df = df |> 
  dplyr::select(Xutm, Yutm, Ba_pref) |> 
  na.omit()

# --- KROK 2: Konwersja na obiekt przestrzenny ---
# Przekształcamy ramkę danych (data.frame) na obiekt sf (simple features)
# coords: wskazujemy kolumny ze współrzędnymi geograficznymi
# crs: definiujemy układ współrzędnych (EPSG:32633 to UTM strefa 33N)
punkty_pref = st_as_sf(df, coords = c("Xutm", "Yutm"), crs = "EPSG:32633") 
# Używamy biblioteki tmap do stworzenia mapy lokalizacyjnej
tm_shape(punkty_pref) + 
  # Rysujemy punkty, gdzie kolor wypełnienia (fill) zależy od wartości Ba_pref
  tm_symbols(
    fill = 'Ba_pref', 
    size = 0.6, 
    # Definiujemy skalę kolorów (odwrócona paleta Spectral dla lepszej czytelności)
    fill.scale = tm_scale(values = "-brewer.spectral", n = 7)
  )

2 Rozgrupowanie danych

Rozgrupowanie danych (ang. declustering) to procedura statystyczna stosowana w geostatystyce, której celem jest zniwelowanie wpływu nierównomiernego rozmieszczenia punktów pomiarowych na szacowanie parametrów całej populacji (np. średniej) (Isaaks i Srivastava, 1989; Rozdział 10). Istnieją dwie główne metody rozgrupowania danych:

  • rozgrupowanie poligonalne - waga punktu jest proporcjonalna do powierzchni obszaru, który mu przypada.
  • rozgrupowanie komórkowe - nadaje wagę punktom na podstawie liczby obserwacji w regularnych komórkach siatki.

2.1 Rozgrupowanie poligonalne

Rozgrupowanie poligonalne (ang. polygonal method) polega na wyznaczeniu wokół każdego punktu pomiarowego tzw. wieloboku wpływu, wewnątrz którego każdy fragment terenu leży bliżej danego punktu pomiarowego niż jakiegokolwiek innego. Wagi rozgrupowujące są ustalane na podstawie powierzchni tych wieloboków:

  • punkty w skupiskach otrzymują małe wieloboki, a co za tym idzie - niskie wagi,
  • punkty odizolowane reprezentujące rzadziej opróbowane obszary otrzymują duże wieloboki i wysokie wagi.

Procedura obliczeniowa dla rozgrupowania poligonalnego

Procedura rozgrupowania poligonalnego składa się z kilku kroków:

1. Wyznaczenie poligonu “wpływu” wokół każdego punktu, np. stosując metodę poligonów Woronoja.

Do stworzenia poligonów wykorzystujemy funkcję voronoi() z pakietu dismo. Musimy najpierw dodać unikalny identyfikator, aby móc później połączyć wyniki z punktami.

library(dismo)
# Dodanie unikalnego ID dla każdego punktu (potrzebne do późniejszego złączenia danych)
punkty_pref$id = 1:nrow(punkty_pref)

# Tworzenie poligonów Voronoi:
# 1. as(..., "Spatial") - konwersja sf do starszego formatu sp (wymagane przez dismo)
# 2. voronoi(...) - generowanie geometrii poligonów
# 3. st_as_sf(...) - powrót do formatu sf
v = st_as_sf(voronoi(as(punkty_pref, "Spatial")))

2. Wyliczenie powierzchni poligonu wpływu.

# Obliczamy powierzchnię każdego wieloboku (st_area zwraca wynik z jednostką, np. m2)
# as.numeric() usuwa jednostkę, pozostawiając samą liczbę
v$pow = as.numeric(st_area(v))

3. Obliczenie wagi każdego punktu

Waga dla punktu obliczana jest poprzez podzielenie powierzchni poligonu przez powierzchnię całego obszaru, a następnie pomnożenie przez liczbę punktów.

\[w'_j=\frac{area_j}{\sum_{j=1}^{n}area_j} \cdot n\]

, gdzie \(area_j\) to powierzchnia poligonu dla lokalizacji \(j\), \(n\) to łączna liczba obserwacji.

Waga o wartości 1 oznacza, że punkt reprezentuje dokładnie średnią powierzchnię przypadającą na jedną próbkę. Waga > 1 to punkt odizolowany, a < 1 to punkt w klastrze.

# Obliczamy wagę: (powierzchnia_indywidualna / powierzchnia_całkowita) * liczba_punktów
v$waga_rp = v$pow / sum(v$pow) * nrow(punkty_pref)

# Przygotowujemy tabelę pomocniczą (usuwamy geometrię, zostawiamy tylko ID i wagę)
v_df = st_drop_geometry(v[c("id", "waga_rp")])

# Łączymy wyliczone wagi z naszym pierwotnym zbiorem punktów na podstawie ID
punkty_pref = merge(punkty_pref, v_df, by = "id")
head(punkty_pref, n = 3)
Simple feature collection with 3 features and 3 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 617722.6 ymin: 5821319 xmax: 620601.7 ymax: 5821420
Projected CRS: WGS 84 / UTM zone 33N
  id Ba_pref  waga_rp                 geometry
1  1      24 4.205149 POINT (617722.6 5821319)
2  2      35 1.236913 POINT (619546.5 5821335)
3  3      23 1.297242 POINT (620601.7 5821420)

2.2 Rozgrupowanie komórkowe

W rozgrupowaniu komórkowym (ang. cell declustering method) cały obszar dzielony jest na regularną siatkę kwadratowych (lub prostokątnych) komórek. Waga przypisana do próbki jest odwrotnie proporcjonalna do liczby innych próbek znajdujących się w tej samej komórce. Jeśli w komórce jest wiele punktów, każdy z nich otrzymuje mniejszą wagę.

Uwaga! Wynik tej metody zależy od wybranej wielkości komórki siatki. Zbyt małe komórki (zawierające tylko po 1 punkcie) nie zmienią wag, natomiast zbyt duża komórka (obejmująca cały obszar) nada wszystkim punktom taką samą wagę.

Procedura obliczeniowa dla rozgrupowania komórkowego

Procedura rozgrupowania komórkowego składa się z kilku kroków:

1. Stworzenie regularnej siatki dla badanego obszaru.

Rozmiar komórki (cellsize) powinien być dobrany tak, aby w obszarach rzadziej opróbowanych w komórce znajdował się zazwyczaj jeden punkt.

# Definicja siatki na podstawie zasięgu danych. 
# Tworzymy siatkę o oczku 3000 x 3000 jednostek mapy
siatka_n = st_make_grid(st_bbox(punkty_pref), cellsize = 3000)
siatka_n = st_as_sf(siatka_n)

# Dodajemy kolumnę pomocniczą do zliczania punktów
punkty_pref$liczebnosc = 0
st_crs(siatka_n) = st_crs(punkty_pref)

2. Policzenie liczby obserwacji w każdym oczku siatki

Używamy funkcji aggregate, która sprawdza, ile punktów zawiera się w każdym poligonie siatki.

# Agregujemy punkty do siatki, stosując funkcję length (zliczanie)
siatka_nr = aggregate(punkty_pref["liczebnosc"], by = siatka_n, FUN = length)
tm_shape(siatka_nr) +
  tm_polygons(
    fill = "liczebnosc", 
    fill.scale = tm_scale_intervals(values = "-brewer.spectral", n = 5),
    fill.legend = tm_legend(title = "Liczba punktów")
  ) +
  tm_shape(punkty_pref) +
  tm_dots() +
  tm_title("Rozkład zagęszczenia próbek") +
  tm_layout(legend.outside = TRUE)

3. Obliczenie wagi każdego punktu

\[w'_j=\frac{\frac{1}{n_i}}{\text{liczba komorek z danymi}} \cdot n\]

# Usuwamy kolumnę pomocniczą, aby nie dublować danych przy złączaniu
punkty_pref$liczebnosc = NULL

# Przypisujemy informację o liczebności komórki z powrotem do punktów (spatial join)
punkty_pref = st_join(punkty_pref, siatka_nr)

# Obliczamy wagi rk (rozgrupowanie komórkowe)
n_punktow = nrow(punkty_pref)
n_komorek_z_danymi = sum(!is.na(siatka_nr$liczebnosc))

punkty_pref$waga_rk = ((1 / punkty_pref$liczebnosc) / n_komorek_z_danymi) * n_punktow
head(punkty_pref, 3)
Simple feature collection with 3 features and 5 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 617722.6 ymin: 5821319 xmax: 620601.7 ymax: 5821420
Projected CRS: WGS 84 / UTM zone 33N
  id Ba_pref  waga_rp liczebnosc                 geometry   waga_rk
1  1      24 4.205149          4 POINT (617722.6 5821319) 0.7091837
2  2      35 1.236913          4 POINT (619546.5 5821335) 0.7091837
3  3      23 1.297242          6 POINT (620601.7 5821420) 0.4727891

Warto zauważyć, że w samej procedurze krigingu, algorytm obliczający wagi estymacji (\(\lambda\)) automatycznie uwzględnia konfigurację przestrzenną danych, w tym ich potencjalne skupienia, co pozwala na poprawne szacowanie wartości nawet przy nieregularnym próbkowaniu.

Niemniej jednak, rozgrupowanie (poligonalne lub komórkowe) pozostaje kluczowe dla:

  • Statystyki opisowej: Poprawne obliczenie średniej i wariancji dla całego obszaru.

  • Weryfikacji modeli: Porównywanie wyników estymacji z “prawdziwą” średnią z próby.

2.3 Porównanie wartości średnich

Celem rozgrupowania jest zbliżenie statystyk obliczonych z próby obciążonej do wartości rzeczywistych. Jeśli próbkowanie było preferencyjne (zagęszczone w rejonach wysokich wartości), spodziewamy się, że:

  • Średnia arytmetyczna (srednia_pref) będzie znacznie zawyżona.

  • Średnie po rozgrupowaniu (srednia_rp, srednia_rk) będą niższe i bardziej zbliżone do wartości referencyjnej (\(41,267\)).

# Tworzymy wektor porównawczy z trzema różnymi statystykami
porownanie_srednich = c(
  # 1. Zwykła średnia arytmetyczna (obarczona błędem próbkowania preferencyjnego)
  srednia_pref = mean(punkty_pref$Ba_pref),
  
  # 2. Średnia ważona metodą poligonową (wagi wynikające z powierzchni poligonów Voronoi)
  srednia_rp = mean(punkty_pref$Ba_pref * punkty_pref$waga_rp, na.rm = TRUE),
  
  # 3. Średnia ważona metodą komórkową (wagi wynikające z zagęszczenia punktów w oczkach siatki)
  srednia_rk = mean(punkty_pref$Ba_pref * punkty_pref$waga_rk, na.rm = TRUE)
)

# Wyświetlenie wyników
print(porownanie_srednich)
srednia_pref   srednia_rp   srednia_rk 
    51.53957     42.43672     43.26675