Wczytanie danych

Procedura wczytania danych wygląda identycznie tak jak w poprzednim przypadku.

library("terra")
files = list.files("dane/landsat", pattern = "\\.TIF$", full.names = TRUE)
landsat = rast(files)
names(landsat) = paste0("B", 1:7)
poly = vect("dane/powiat_sremski.gpkg")

Jednak tym razem dodatkowo wczytamy raster z klasami pokrycia terenu (Sentinel-2 Global Land Cover).

cat = rast("dane/S2GLC_T33UXT.tif")
cat
## class       : SpatRaster 
## dimensions  : 5559, 8535, 1  (nrow, ncol, nlyr)
## resolution  : 10.0005, 10.00056  (x, y)
## extent      : 602547.5, 687901.7, 5742527, 5798120  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 33N (EPSG:32633) 
## source      : S2GLC_T33UXT.tif 
## name        : S2GLC_T33UXT

Prosta wizualizacja klas pokrycia terenu.

plot(cat, type = "classes")

Klasy przedstawione są jako ID, co nie jest dla nas zbyt przydatne. Oprócz danych rastrowych, musimy wczytać jeszcze tabelę atrybutów, która znajduje się w osobnym pliku S2GLC_T33UXT.csv. Możemy to zrobić za pomocą funkcji read.csv().

leg = read.csv("dane/S2GLC_T33UXT.csv")
head(leg)
##   ID     RGB           Klasa
## 1  0 #FFFFFF          Chmury
## 2 62 #D20000        Zabudowa
## 3 73 #FDD327 Obszary uprawne
## 4 75 #B05B10         Winnice
## 5 82 #239800  Lasy liściaste
## 6 83 #086209    Lasy iglaste

W pierwszej kolumnie znajduje się ID klasy, w drugiej wartość koloru RGB w zapisie szesnastkowym, w trzeciej nazwa klasy.

Operacje na rastrach

Jak można zauważyć, nasze rastry różnią się rozdzielczością przestrzenną (scena Landsat ma 30 m, natomiast klasy pokrycia terenu 10 m). Kiedy rastry posiadają różne rozdzielczości, to nie jest możliwe wykonywanie na nich operacji matematycznych. W takim przypadku musimy sprowadzić je do jednakowej rozdzielczości. Ten proces nazywa się przepróbkowaniem (resampling).

Przepróbkowanie można wykonać:

  1. Z większej rozdzielczości do mniejszej, np. 100 m -> 500 m (downsampling, downscaling).
  2. Z mnieszej rozdzielczości do większej, np. 500 m -> 100 m (upsampling, upscaling).

Z racji iż proste metody przepróbkowania nie zwiększają ilości informacji przy zwiększaniu rozdzielczości przestrzennej, to lepiej wykonać przepróbkowanie do niższej rozdzielczości.

Przepróbkowanie w terra można wykonać przy pomocy funkcji resample(). Dostępne są różne metody przepróbkowania, ale w przypadku danych kategorycznych, koniecznie trzeba wykorzystać algorytm najbliższego sąsiada (method = "near") lub wartość modalną, dominantę (method = "mode"). Jeśli tego nie zrobimy, to ID kategorii zostaną zmienione.

cat = resample(cat, landsat, method = "near")
res(cat) # wyświetl rozdzielczość po przepróbkowaniu
## [1] 30 30

Kolejne operacje (czyli przycinanie, maskowanie, skalowanie i usunięcie wartości odstających) przebiegają tak jak w poprzednim przykładzie.

# przycinanie, maskowanie
landsat = crop(landsat, poly, mask = TRUE)
cat = crop(cat, poly, mask = TRUE)

# skalowanie
landsat = landsat * 2.75e-05 - 0.2

# usunięcie wartości odstających
landsat[landsat < 0] = NA
landsat[landsat > 1] = NA

Przygotowanie danych

Teraz przygotujmy nasze dane do klasyfikacji. Tym razem użyjemy modelu drzewa decyzyjnego, który wymaga danych w postaci ramki danych (data frame). Jedna kolumna to zmienna modelowana/zależna (u nas klasy pokrycia terenu), a pozostałe kolumny to zmienne wyjaśniające/niezależne (kanały spektralne).

data = cbind(values(cat), values(landsat)) # połączenie kolumn w macierzy
data = as.data.frame(data) # konwersja macierzy do ramki danych
data = na.omit(data) # usunięcie brakujących wartości

Jedna z klas w zbiorze danych to chmury, które zostały sklasyfikowane podczas tworzenia mapy z klasami pokrycia terenu. Wiadomo, że zachmurzenie jest zmienne w czasie, dlatego powinniśmy usunąć tę klasę (ID = 0).

# usuń piksele reprezentujące klasę chmury
# ! to operator negacji (NOT)
data = data[!data$S2GLC_T33UXT == 0, ]

W celu ułatwienia analizy danych, możemy zmienić ID klasy na nazwę. Dopasowanie nazw klas do ID można wykonać za pomocą funkcji merge(). Dodatkowo, zmieńmy jeszcze nazwę kolumny z S2GLC_T33UXT na klasa.

data = merge(data, leg[, -2], by.x = "S2GLC_T33UXT", by.y = "ID")
data = data[, -1] # usuń pierwszą kolumnę z ID klasy
colnames(data)[8] = "klasa" # zmień nazwę kolumny
data$klasa = as.factor(data$klasa) # zmień typ danych na kategoryczny

Dane wejściowe wyglądają teraz następująco.

head(data)
##          B1        B2        B3        B4        B5       B6        B7    klasa
## 1 0.0547875 0.0691150 0.1028300 0.1315400 0.2248750 0.281140 0.2568575 Zabudowa
## 2 0.0206875 0.0266550 0.0520650 0.0432375 0.2762725 0.146555 0.0936725 Zabudowa
## 3 0.0377650 0.0485725 0.0869900 0.0841575 0.3838800 0.217615 0.1450150 Zabudowa
## 4 0.0580325 0.0659800 0.0898775 0.0866875 0.1910500 0.192205 0.1502950 Zabudowa
## 5 0.0558050 0.0724425 0.1128950 0.1346475 0.2250125 0.265960 0.2239675 Zabudowa
## 6 0.0285250 0.0311100 0.0529725 0.0406250 0.3414475 0.151945 0.0838550 Zabudowa

Przed klasyfikacją warto sprawdzić częstość występowania poszczególnych kategorii. Jeśli część kategorii pojawia się bardzo często, a niektóre prawie wcale, to wtedy mamy problem z niezbalansowanym zbiorem danych. W takiej sytuacji, kiedy model posiada zbyt mało przykładów którejś klasy (np. bagna), to nie jest możliwe żeby nauczył się rozpoznawać tę klasę. Oprócz tego, wynik jakości klasyfikatora jest zbyt optymistyczny (tzn. w rzeczywistości działa gorzej niż na zbiorze uczącym).

Częstość występowania klas można sprawdzić za pomocą funkcji table() na kolumnie kategorycznej i następnie zamienić to na postać procentową używając prop.table().

tabela = table(data$klasa)
prop.table(tabela) * 100
## 
##                Bagna         Lasy iglaste       Lasy liściaste 
##            4.0356612           11.9760694           10.0572456 
##      Obszary uprawne Roślinność trawiasta           Torfowiska 
##           51.4040823           15.2939704            0.6978963 
##             Zabudowa      Zbiorniki wodne 
##            5.5545476            0.9805271

Klasyfikacja

# install.packages(c("rpart", "rpart.plot"))
library("rpart") # model klasyfikacyjny
library("rpart.plot") # wizualizacja modelu

Każdy opracowany model powinien zostać poddany walidacji (weryfikacji) na niezależnym zbiorze danych. Oznacza to, że powinniśmy sprawdzić poprawność prognozowania naszego modelu na innym zbiorze danych, który nie został wykorzystany na etapie modelowania (uczenia). Innymi słowy, chcemy wykluczyć sytuację, w której nasz model mógłby nauczyć się jednie klasyfikować obiekty na naszym obszarze analizy, a zwracałby błędne wyniki na innym (lecz podobnym) obszarze.

Istnieje wiele metod walidacji, ale wykorzystamy tutaj najprostszą polegającą na podziale wejściowego zbioru danych na treningowy oraz testowy. Generalnie przyjmuje się, że proporcja między tymi dwoma zestawami powinna wynosić około 70-80% do 30-20%. Tak jak poprzednio musimy wylosować próbę o ustalonej wielkości używając funkcji sample().

# podział na zbiór treningowy i testowy
set.seed(1) # ziarno losowości
n = round(0.7 * nrow(data)) # wielkość próby 70%
trainIndex = sample(nrow(data), size = n) # wylosuj indeksy
train = data[trainIndex, ] # wybierz próbki treningowe
test = data[-trainIndex, ] # wybierz próbki testowe (nietreningowe)

Po tej czynności dochodzimy do najważniejszego etapu, czyli trenowania modelu klasyfikacyjnego. W tym celu wykorzystamy funkcję rpart(), która wymaga zdefiniowania:

  1. Zmiennej zależnej i zmiennych zależnych za pomocą odpowiedniej formuły.
  2. Zbioru danych treningowych (data = train).
  3. Metody (method = "class").

Odnośnie punkty pierwszego, formułę można zdefiniować na dwa sposoby:

  1. Używając nazw poszczególnych zmiennych: klasa ~ B1 + B2 + B3 + B4 + B5 + B6 + B7.
  2. Używając kropki: klasa ~ .. Kropka zastępuje wszystkie nazwy zmiennych wyjaśniających z ramki danych.

Znak ~ (tylda) oznacza “jest zależne od”, czyli klasa pokrycia terenu jest zależna od kanałów B1 do B7.

mdl = rpart(klasa ~ ., data = train, method = "class")

Po zakończeniu tej operacji możemy sprawdzić jakich reguł klasyfikacyjnych nauczył się model. Drzewo decyzyjne można zwizualizować za pomocą funkcji prp().

prp(mdl)

Kluczowym etapem modelowania jest walidacja modelu. Sprawdźmy zatem jaka jest jego skuteczność. W tym celu musimy wykonać predykcję dla zbioru testowego (funkcja predict()). Zbiór danych testowych musi posiadać dokładnie te same zmienne wyjaśniające co zbiór treningowy. Kolumnę z prawdziwymi (rzeczywistymi) klasami należy usunąć.

# walidacja dla zbioru testowego
pred = predict(mdl, test[, -8], type = "class")
unname(head(pred)) # `unname()` usuwa numer porządkowy wiersza/piksela
## [1] Obszary uprawne Obszary uprawne Obszary uprawne Obszary uprawne
## [5] Zabudowa        Obszary uprawne
## 8 Levels: Bagna Lasy iglaste Lasy liściaste ... Zbiorniki wodne

Wykonaliśmy predykcje dla zbioru testowego. Teraz musimy obliczyć wybraną miarę skuteczności. Jako przykład wybierzemy dokładność (accuracy) definiowaną jako iloraz poprawnych klasyfikacji do wszystkich (poprawnych i niepoprawnych) klasyfikacji.

pop_klas = test$klasa == pred # zwraca wartość logiczną czy klasa jest prawidłowa
sum(pop_klas) / length(pop_klas)
## [1] 0.7138344

Skuteczność naszego modelu wynosi około 71%. Oprócz jednej ogólnej statystyki możemy sprawdzić również błędy klasyfikacji dla poszczególnych klas. Takie zestawienie nazywane jest tabelą pomyłek (confusion matrix).

table(pred = pred, true = test$klasa)
##                       true
## pred                   Bagna Lasy iglaste Lasy liściaste Obszary uprawne
##   Bagna                 1339          150             85              26
##   Lasy iglaste           251        17542           2298             106
##   Lasy liściaste        2133         3159          12183             709
##   Obszary uprawne       2034         1375           3563           93946
##   Roślinność trawiasta   934          229            405            2855
##   Torfowiska               0            0              0               0
##   Zabudowa               900          431            701            1017
##   Zbiorniki wodne        133            5              0               3
##                       true
## pred                   Roślinność trawiasta Torfowiska Zabudowa Zbiorniki wodne
##   Bagna                                  59         23       74             383
##   Lasy iglaste                           69        195       99              12
##   Lasy liściaste                       1058        526      441              10
##   Obszary uprawne                     20077        482     6554               3
##   Roślinność trawiasta                 7190         49      211               0
##   Torfowiska                              0          0        0               0
##   Zabudowa                              933         94     3286               8
##   Zbiorniki wodne                        11          0       15            1431

Poprawnie sklasyfikowane obiekty znajdują się na przekątnej. Obiekty sklasyfikowane jako fałszywe pozytywne (false positive) znajdują się w prawej górnej części, natomiast obiekty sklasyfikowane jako fałszywe negatywne (false negative) w lewej dolnej części.

Jeśli opracowany model spełnia nasze oczekiwania, to możemy wykorzystać go do predykcji na całym obszarze. Ponownie wykorzystamy funkcję predict(), ale tym razem jako dane wejściowe użyjemy nasz raster landsat. Dodatkowo, powinniśmy ustawić argument na.rm = TRUE, aby uniknąć predykcji poza obszarem analizy.

pred_map = predict(landsat, mdl, type = "class", na.rm = TRUE)

Jako wynik powyższej operacji otrzymaliśmy mapę z predykcją klas pokrycia terenu. Teraz możemy dokonać wizualizacji. Jednak tym razem musimy użyć odpowiedniego schematu kolorów, który znajduje się w obiekcie leg. Tę operację należy wykonać w dwóch krokach:

  1. Sprawdzić, które klasy rzeczywiście występują na naszym obszarze.
  2. Dopasować kolory do odpowiednich klas.

Odnośnie pierwszego punktu, to najprościej wykorzystać funkcję droplevels(), która usunie puste kategorie z rastra (np. winnica, wrzosowisko). Następnie należy wyświetlić kategorie za pomocą funkcji levels(). Do każdej warstwy rastra mogą być przypisane różne kategorie, więc musimy pobrać je tylko dla pierwszej warstwy oraz wskazać atrybut class.

Dopasowania kolorów można dokonać za pomocą funkcji match(). Jako pierwszy wskazujemy obiekt z naszymi klasami, jako drugi obiekt ramkę danych ze schematem kolorów (leg) i w wyniku otrzymamy dopasowane indeksy kolorów.

lv = droplevels(pred_map) # usuń puste kategorie z rastra
lv = levels(lv) # zwróć kategorie jako ramkę danych
lv = lv[[1]][["class"]] # wybierz pierwszą warstwę i kolumnę z nazwami klas
col_idx = match(lv, leg$Klasa) # dopasuj klasy do odpowiednich kolorów
plot(pred_map, main = "Predykcja klas", col = leg$RGB[col_idx])

W przypadku, gdy raster posiada oryginalne ID klas (tj. 0, 62, 73, 75, itd.), to możemy po prostu przypisać ramkę danych do rastra.

levels(cat) = leg[, c(1, 3)]
lv = levels(droplevels(cat))[[1]][["Klasa"]]
col_idx = match(lv, leg$Klasa)
plot(cat, col = leg$RGB[col_idx], main = "Rzeczywiste klasy")

Podsumowanie

Niniejsze warsztaty stanowią jedynie zarys wykorzystania metod uczenia maszynowego do analizy danych przestrzennych. Oprócz problemu klasyfikacji możemy również rozwiązywać problemy z zakresu regresji, wykrywania anomalii, redukcji wymiarowości, itd., używając bardziej zaawansowanych modeli (np. lasy losowe, modele wzmacnianie, sieci neuronowe, uczenie głębokie).

Co więcej, cała analiza może być bardziej zaawansowana. Przedstawione treści można rozwinąć o kolejne zagadnienia dotyczące:

  • automatycznego wyszukiwania i pobierania danych satelitarnych,
  • integracji różnych sensorów (satelity + drony),
  • przetwarzania danych w chmurze,
  • optymalizacji parametrów modeli,
  • innych metod walidacji wyników (w tym w ujęciu przestrzennym),
  • oceny istotności zmiennych wyjaśniających.

Należy podkreślić, że zdjęcia satelitarne to nie jest jedyne źródło danych przestrzennych i do analizy można wykorzystać zupełnie inne dane, np. cyfrowe modele wysokościowe czy zasięgi występowania gatunków.

Jeśli zainteresowałeś się tym tematem, to jego rozwinięcie znajdziesz w poniższych podręcznikach (w j. angielskim):

Warto również sprawdzić inne pakiety do analizy danych przestrzennych oraz uczenia maszynowego.

W razie problemów pomoc najlepiej szukać na StackOverflow, lub RStudio Community. Często ciekawe dyskusje pojawiają się na Twitterze (#rspatial) oraz GitHubie (https://github.com/rspatial; https://github.com/r-spatial).