Klasyfikacja nienadzorowana

library("terra")
library("rstac")
options(timeout = 600)

Pozyskanie danych

Niniejszą analizę wykonamy na podstawie danych z Sentinela-2. Jako źródło danych wykorzystamy katalog STAC i usługę Earth Search. Przeprowadzimy również proste filtrowanie dostępnych produktów definiując w zapytaniu: kolekcję (collections), zakres przestrzenny w układzie WGS 84 (bbox), interwał czasowy w standardzie RFC 3339 (datetime) oraz atrybut zachmurzenia sceny (eo:cloud_cover).

stac_source = stac("https://earth-search.aws.element84.com/v1")
stac_source |>
  stac_search(
    collections = "sentinel-2-c1-l2a",
    bbox = c(22.5, 51.1, 22.6, 51.2),
    datetime = "2023-03-01T00:00:00Z/2023-10-31T00:00:00Z") |>
  ext_query(`eo:cloud_cover` < 10) |>
  post_request() -> obrazy

Jako wynik zapytania zostały zwrócone metadane scen spełniające zadane warunki. Następnie możemy wybrać przykładową scenę z 21 września 2023 r. oraz ograniczyć liczbę kanałów do czterech podstawowych w rozdzielczości 10 m, tj. niebieskiego, zielonego, czerwonego i bliskiej podczerwieni w celu uproszczenia analizy. Finalnie, używając funkcji assets_url() zostaną zwrócone odpowiednie odnośniki do pobrania rastrów.

kanaly = c("blue", "green", "red", "nir")

obrazy |>
  items_filter(properties$`s2:tile_id` == "S2A_OPER_MSI_L2A_TL_2APS_20230921T151458_A043076_T34UEB_N05.09") |>
  assets_select(asset_names = kanaly) |>
  assets_url() -> sentinel
sentinel

## [1] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B02.tif"
## [2] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B03.tif"
## [3] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B08.tif"
## [4] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/9/S2A_T34UEB_20230921T094329_L2A/B04.tif"

Zauważ, że kanał 8 (bliska podczerwień) jest przed kanałem 4 (czerwony). Należy mieć na uwadze, że nieodpowiednia kolejność może spowodować problemy w dalszej części projektu.

W kolejnym kroku, używając funkcji dir.create() stworzymy nowy katalog na dysku, do którego zostaną pobrane zobrazowania. Oprócz tego, należy zdefiniować ścieżki i nazwy plików. W tym celu będą przydatne dwie funkcje – file.path() do stworzenia ścieżek do plików oraz basename() do wyodrębnienia nazwy pliku wraz z rozszerzeniem z URL.

dir.create("sentinel")
rastry = file.path("sentinel", basename(sentinel))

Teraz możemy pobrać nasze dane używając funkcji download.file() w pętli.

for (i in seq_along(sentinel)) {
  download.file(sentinel[i], rastry[i], mode = "wb")
}

Przygotowanie danych

Po pobraniu danych musimy stworzyć listę plików (rastrów), które zamierzamy wczytać. W tym celu możemy wykorzystać funkcję list.files(), która jako argument przyjmuje ścieżkę do folderu z plikami. Oprócz tego musimy wskazać jaki rodzaj plików chcemy wczytać (pattern = "\\.tif$") oraz zwrócić pełne ścieżki do plików (full.names = TRUE).

rastry = list.files("sentinel", pattern = "\\.tif$", full.names = TRUE)
rastry

## [1] "sentinel/B02.tif" "sentinel/B03.tif" "sentinel/B04.tif" "sentinel/B08.tif"

Kiedy utworzyliśmy już listę plików, możemy je wczytać za pomocą funkcji rast(). Pobrane zobrazowania pokrywają duży obszar (około 12 000 km\(^2\)). Dla uproszczenia analizy zdefiniujmy mniejszy zakres przestrzenny do wczytania używając funkcji ext() oraz przekazując obiekt SpatExtent jako argument win w funkcji rast().

Zwróć uwagę, że do zdefiniowania zakresu przestrzennego podczas wyszukiwania danych użyliśmy układu WGS 84, natomiast teraz wymagany jest rzeczywisty układ rastra, tj. EPSG:32634. Można to sprawdzić w metadanych katalogu STAC (obiekt obrazy) lub używając funkcji describe() (wymaga ścieżki do pliku).

bbox = ext(505000, 555000, 5629000, 5658000)
r = rast(rastry, win = bbox)

Możemy również zmienić nazwy kanałów spektralnych. Przed tą operacją należy się upewnić czy kanały zostały wczytane w prawidłowej kolejności.

names(r) = kanaly
r

## class       : SpatRaster 
## dimensions  : 2900, 5000, 4  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## window      : 505000, 555000, 5629000, 5658000  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 34N (EPSG:32634) 
## sources     : B02.tif  
##               B03.tif  
##               B04.tif  
##               B08.tif  
## names       :   blue,  green,    red,    nir 
## min values  : >    1, >  352, >  677, >  910 
## max values  : 18560<, 17648<, 17056<, 16514<

Na pierwszy rzut oka, metadane odnoszące się do wartości minimalnych i maksymalnych sugerują, że odbicie spektralne jest z zakresu od jeden do kilkunastu tysięcy. Wyliczmy jednak statystyki z próby używając funkcji summary().

summary(r)

## Warning: [summary] used a sample

##       blue              green              red               nir        
##  Min.   :-0.00220   Min.   :0.00580   Min.   :0.00340   Min.   :0.0021  
##  1st Qu.: 0.01360   1st Qu.:0.02910   1st Qu.:0.01840   1st Qu.:0.1915  
##  Median : 0.02670   Median :0.05430   Median :0.04200   Median :0.2301  
##  Mean   : 0.03305   Mean   :0.05719   Mean   :0.05623   Mean   :0.2420  
##  3rd Qu.: 0.04880   3rd Qu.:0.07960   3rd Qu.:0.08720   3rd Qu.:0.2829  
##  Max.   : 0.64840   Max.   :0.68400   Max.   :0.64760   Max.   :0.6684

Możemy zauważyć rozbieżność pomiędzy metadanymi a rzeczywistymi statystykami. Aby ograniczyć zajmowane miejsce, pobrane rastry zapisane są jako wartości całkowite (UInt16) z pewnym współczynnikiem skalowania (scale) oraz przesunięciem (offset). Funkcja rast() automatycznie wczytuje te parametry z metadanych pliku i następnie dokonuje przekształcenia wartości według poniższego wzoru:

\[ x = x \cdot scale - offset \]

Parametry te można pozyskać za pomocą funkcji scoff() (lub funkcji describe()). Należy również pamiętać, że każdy produkt/kolekcja ma swoje unikalne parametry i konieczne jest zapoznanie się z dokumentacją.

scoff(r)

##      scale offset
## [1,] 1e-04   -0.1
## [2,] 1e-04   -0.1
## [3,] 1e-04   -0.1
## [4,] 1e-04   -0.1

Zasadniczo, wartości odbicia spektralnego mieszczą się w przedziale od 0 do 1, gdzie 0 oznacza brak odbicia (całe światło zostało pochłonięte), a 1 oznacza całkowite odbicie od powierzchni. W praktyce obiekty nie odbijają bądź pochłaniają stuprocentowo światła, niemniej sensory oraz procesy kalibracyjne nie są idealne, więc mogą pojawić się wartości odstające od tego przedziału, tak jak w tej sytuacji.

Można ten problem rozwiązać na dwa sposoby:

1. Zastąpić te wartości brakiem danych (NA).
2. Dociąć do minimalnej i maksymalnej wartości.

Pierwszy sposób może spowodować, że stracimy dużą część zbioru danych. Natomiast drugi sposób może powodować przekłamania.

# sposób 1
r[r < 0] = NA
r[r > 1] = NA

# sposób 2
r[r < 0] = 0
r[r > 1] = 1

Po przeskalowaniu wartości możemy wyświetlić kompozycję RGB. W tym przypadku zamiast funkcji plot() należy użyć funkcji plotRGB() oraz zdefiniować kolejność kanałów czerwonego, zielonego oraz niebieskiego. Często zdarza się, że kompozycje są zbyt ciemne/jasne, wtedy warto zastosować rozciągnięcie kolorów używając argumentu stretch = "lin" lub stretch = "hist".

# plotRGB(r, r = 3, g = 2, b = 1)
plotRGB(r, r = 3, g = 2, b = 1, stretch = "lin")

Klasteryzacja

library("cluster") # klasteryzacja danych

Dane do modelowania muszą zostać przygotowane w odpowiedni sposób. Modele klasyfikacyjne najczęściej na etapie trenowania wymagają macierzy lub ramki danych (data frame). Dane rastrowe można przetworzyć do macierzy przy użyciu funkcji values().

mat = values(r)
nrow(mat) # wyświetla liczbę wierszy

## [1] 14500000

Za pomocą interaktywnej funkcji View() możemy sprawdzić jak wygląda nasza macierz.

View(mat)

W macierzy występują brakujące wartości. Zazwyczaj modele nie obsługują NA, więc musimy je usunąć. Służy do tego dedykowana funkcja na.omit().

mat_omit = na.omit(mat)
nrow(mat_omit)

## [1] 14499310

Teraz przejdziemy do kolejnego etapu analizy, jakim jest wytrenowanie modelu. Istnieje wiele metod i modeli grupowania (patrz CRAN Task View), ale w tym przykładzie użyjemy prostego modelu grupowania metodą k-średnich. Ten model wymaga jedynie, aby podać z góry liczbę grup/klastrów (argument centers). Jest to algorytm stochastyczny, więc za każdym razem zwraca inne wyniki. Żeby analiza była powtarzalna musimy ustawić ziarno losowości – set.seed().

set.seed(123) # ziarno losowości
mdl = kmeans(mat_omit, centers = 4)

W wyniku powyższej operacji otrzymaliśmy m.in.:

1. Obliczone średnie wartości grup dla poszczególnych kanałów (mdl$centers).
2. Wektor ze sklasyfikowanymi wartościami macierzy (mdl$cluster).

Wyświetlmy te obiekty:

mdl$centers

##         blue      green        red       nir
## 1 0.02463308 0.04917406 0.03927713 0.2775827
## 2 0.01521135 0.02944429 0.02342255 0.1818289
## 3 0.02945364 0.06482988 0.04082320 0.3849277
## 4 0.06184890 0.09154669 0.11423999 0.2123605

head(mdl$cluster)

## [1] 2 1 3 4 4 1

Oznacza to, że pierwszy wiersz (reprezentujący pojedyncze oczko siatki) należy do grupy 2, drugi do grupy 1, trzeci do grupy 3, itd.

Walidacja

Nieodłącznym elementem modelowania jest walidacja opracowanych modeli. Wyzwaniem jest wybór właściwej metody grupowania dla konkretnego zbioru danych i określenie odpowiedniej liczby grup. Należy pamiętać, że zwiększenie liczby klastrów zwiększa podobieństwo między obiektami w klastrze, ale przy ich większej liczbie interpretacja staje się trudniejsza.

Najczęstszym sposobem walidacji wyników grupowania jest użycie wewnętrznych metryk, takich jak wskaźnik Dunna, wskaźnik Daviesa-Bouldina lub wskaźnik sylwetki (silhouette index). Na potrzebny niniejszej analizy użyjemy tego ostatniego.

Indeks sylwetki ocenia zbieżność i separację klastrów na podstawie odległości między obiektami w tym samym klastrze i między obiektami w różnych klastrach. Wartości tego wskaźnika mieszczą się w zakresie od -1 do 1. Wartość bliska 1 wskazuje, że obiekt jest dobrze zgrupowany i znajduje się daleko od sąsiednich klastrów. Wartość bliska -1 sugeruje, że obiekt mógł zostać przypisany do niewłaściwego klastra. Wartość bliska 0 wskazuje, że obiekt znajduje się bardzo blisko granicy pomiędzy różnymi klastrami. Ogólnie rzecz biorąc, wyższa wartość tego wskaźnika wskazuje na lepsze wyniki grupowania. Więcej szczegółów można znaleźć w dokumentacji cluster::silhouette().

Spróbujmy teraz obliczyć wartości tego wskaźnika. Zasadniczo wymaga to obliczenie podobieństwa każdego obiektu do każdego obiektu, co w naszym przypadku jest zadaniem niemożliwym (nasz zbiór danych składa się z ponad 14 mln obiektów). Aby to wykonać, musimy wykorzystać mniejszą próbkę (załóżmy \(n=10000\)). W funkcji musimy określić dwa obiekty –- wektor z klastrami oraz macierz niepodobieństwa, którą można wcześniej obliczyć za pomocą funkcji dist().

set.seed(123)
# losowanie indeksów
idx = sample(1:nrow(mat_omit), size = 10000)
head(idx)

## [1] 10406351 11723278 12716970  3269750  2782437  7705241

# obliczenie wskaźnika sylwetki
sil = silhouette(mdl$cluster[idx], dist(mat_omit[idx, ]))
summary(sil)

## Silhouette of 10000 units in 4 clusters from silhouette.default(x = mdl$cluster[idx], dist = dist(mat_omit[idx, ])) :
##  Cluster sizes and average silhouette widths:
##      2972      3007      1144      2877 
## 0.3752375 0.5016613 0.3752901 0.4159651 
## Individual silhouette widths:
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -0.07105  0.28904  0.47232  0.42498  0.58097  0.68149

Średnia zbieżność klastrów wynosi 0,43. Nie jest to najlepszy wynik (powinniśmy spróbować zwiększyć liczbę klastrów lub użyć innej metody grupowania). Możemy również zaprezentować wyniki na wykresie.

kolory = rainbow(4) # wybierz 4 kolory z wbudowanej palety `rainbow`
plot(sil, border = NA, col = kolory, main = "Silhouette Index")

Interpretacja

Istotą grupowania jest utworzenie grupy podobnych obiektów, natomiast naszym zadaniem jest zinterpretowanie tego, co reprezentują utworzone grupy i nadanie im nazwy. Interpretacja jest trudnym zadaniem i często wyniki są niejasne. W tym celu konieczne jest przeanalizowanie statystyk opisowych klastrów oraz wykorzystanie różnych wykresów i kompozycji map. Bardzo przydatna jest także znajomość właściwości spektralnych obiektów.

Spróbujmy zatem zinterpretować uzyskane skupienia, korzystając z wykresu pudełkowego. Największe możliwości wizualizacji danych dostarcza pakiet ggplot2. Tutaj można znaleźć darmowy podręcznik oraz gotowe “przepisy”. Wymieniony pakiet wymaga przygotowania zbioru danych do odpowiedniej postaci, tj. dane muszą być przedstawione jako ramka danych w tzw. formie długiej (wiele wierszy), podczas gdy standardowe funkcje do wizualizacji wymagają formy szerokiej (wiele kolumn). Takiej konwersji można dokonać w prosty sposób używając pakietu tidyr.

library("tidyr") # transformacja danych
library("ggplot2") # wizualizacja danych

Jak zauważyliśmy wcześniej, nasz zbiór danych jest dość duży i nie ma potrzeby prezentowania wszystkich danych. Możemy to zrobić efektywniej, używając wcześniej wylosowanej próby. Połączmy zatem wylosowane wiersze z macierzy z odpowiadającymi im klastrami (cbind()). Następnie macierz zamienimy na ramkę danych (as.data.frame()).

stats = cbind(mat_omit[idx, ], klaster = mdl$cluster[idx])
stats = as.data.frame(stats)
head(stats)

##     blue  green    red    nir klaster
## 1 0.0196 0.0518 0.0288 0.3321       1
## 2 0.0254 0.0469 0.0484 0.1880       2
## 3 0.0411 0.0876 0.0628 0.3515       3
## 4 0.0429 0.0927 0.0600 0.3505       3
## 5 0.0104 0.0270 0.0168 0.2179       2
## 6 0.0241 0.0469 0.0424 0.2396       1

Wyświetlone dane mają formę szeroką (każdy kanał spektralny zapisany jest w osobnej kolumnie). Teraz musimy zmienić formę, w której otrzymamy dwie kolumn – kanał oraz wartość. W tym celu wykorzystamy funkcję pivot_longer().

stats = pivot_longer(stats, cols = 1:4, names_to = "kanal", values_to = "wartosc")

Dla formalności możemy jeszcze zmienić typ danych (klastrów i kanałów) na kategoryczny (factor). W praktyce związane jest to z uproszczeniem struktury danych (przejście ze skali ilorazowej do nominalnej).

stats$klaster = factor(stats$klaster)
stats$kanal = factor(stats$kanal)
head(stats)

## # A tibble: 6 × 3
##   klaster kanal wartosc
##   <fct>   <fct>   <dbl>
## 1 1       blue   0.0196
## 2 1       green  0.0518
## 3 1       red    0.0288
## 4 1       nir    0.332 
## 5 2       blue   0.0254
## 6 2       green  0.0469

Ramka danych jest już przygotowana. Teraz stwórzmy prosty wykres pudełkowy.

ggplot(stats, aes(x = kanal, y = wartosc, fill = klaster)) +
  geom_boxplot()

Zmieńmy kilka domyślnych parametrów żeby poprawić odbiór ryciny.

etykiety = c("Niebieski", "Zielony", "Czerwony", "Bliska\npodczerwień")

ggplot(stats, aes(x = kanal, y = wartosc, fill = klaster)) +
  geom_boxplot(show.legend = FALSE) +
  scale_x_discrete(limits = kanaly, labels = etykiety) +
  scale_fill_manual(values = kolory) +
  facet_wrap(vars(klaster)) +
  xlab("Kanał") +
  ylab("Odbicie") +
  theme_light()

Na podstawie powyższego wykresu możemy przeanalizować właściwości spektralne klastrów, a tym samym zinterpretować, jakie obiekty reprezentują. Uproszczone wnioski mogą być następujące:

• Klaster 1 – znacząca różnica pomiędzy odbiciem w kanale czerwonym a bliskiej podczerwieni wskazuje na obecność roślinności. W porównaniu do klastra 3, iloraz między tymi kanałami jest mniejszy, co może świadczyć o niższym poziomie biomasy. Ten klaster może reprezentować obszary trawiaste i/lub pastwiska.
• Klaster 2 – ten klaster wygląda podobnie do klastra 2, z tą różnicą, że we wszystkich pasmach występują wartości odstające bliskie zera. Wartości te reprezentują obiekty wodne pochłaniające promieniowanie. Oczywiście nie jest to poprawna klasyfikacja – obiekty wodne i lasy powinny być przypisane do oddzielnych klastrów.
• Klaster 3 – w porównaniu do klastra 1, wyższe wartości odbicia w kanale bliskiej podczerwieni mogą świadczyć o wyższym poziomie biomasy, w związku z czym klaster ten może reprezentować obszary o gęstej roślinności.
• Klaster 4 – wyższe wartości odbicia kanałów spektralnych w porównaniu do pozostałych klastrów wskazują na jasne obiekty. W tym przypadku są to odkryte gleby bez pokrywy roślinnej, które zazwyczaj mają kolor jasnobrązowy lub szary.

Finalna mapa

Ostatnim etapem jest stworzenie mapy klasyfikacyjnej pokrycia terenu na podstawie otrzymanego wektora z klastrami (mdl$cluster). Na początku musimy przygotować pusty wektor składający się z całkowitej liczby komórek rastra. Można to sprawdzić za pomocą funkcji ncell(). W naszym przypadku jest to 14 500 000 komórek.

# przygotuj pusty wektor
wek = rep(NA, ncell(r))

Następnie musimy przypisać nasze grupy w wektorze w odpowiednie miejsca, tj. tym, które nie są zamaskowane (NA). Do niezamaskowanych wartości można odwołać się przez funkcję complete.cases().

# zastąp tylko te wartości, które nie są NA
wek[complete.cases(mat)] = mdl$cluster 

W ostatnim kroku należy skopiować metadane obiektu r, ale tylko z jedną warstwą, i przypisać mu wartości wektora wek.

# stwórz nowy raster
clustering = rast(r, nlyrs = 1, vals = wek)

Zaprezentujmy wynik grupowania na mapie, używając odpowiednich kolorów i nazw klastrów.

kolory = c("#9aed2d", "#086209", "#2fbd2f", "#d9d9d9")
kategorie = c("trawa", "lasy/woda", "gęsta roślinność", "odkryta gleba")
plot(clustering, col = kolory, type = "classes", levels = kategorie,
     mar = c(3, 3, 3, 7))

Jeśli wynik jest zadowalający, możemy zapisać go za pomocą funkcji writeRaster(). Taki plik można później wczytać w R lub innym programie obsługującym dane przestrzenne (np. QGIS). Dodatkowo, w przypadku danych kategorycznych, podczas zapisu warto ustawić typ danych jako Byte / INT1U (pod warunkiem, że liczba kategorii nie przekracza 255).

writeRaster(clustering, "clustering.tif", datatype = "INT1U")

Dalsze kroki

Przedstawiona analiza stanowi jedynie zarys wykorzystania klasyfikacji nienadzorowanej do danych satelitarnych i następujące aspekty mogą zostać rozszerzone:

• Post-processing wyników – zastosowana metoda grupowania opiera się na podejściu pikselowym, co oznacza, że relacje między sąsiadującymi komórkami nie są uwzględniane w klasyfikacji. Powoduje to efekt pieprzu i soli. Można to nieco poprawić, stosując filtr modalny (terra::focal()) lub filtr przesiewowy (terra::sieve()). Innym podejściem jest użycie zmiennych, które są obliczane w ruchomym oknie lub w mniejszych skalach przestrzennych.
• Zwiększenie skuteczności klasteryzacji – w naszym grupowaniu lasy i obiekty wodne zostały błędnie zaliczone do jednego klastra oraz strefa zabudowana nie została w ogóle wyodrębniona. W tym przypadku pojawia się kilka potencjalnych rozwiązań, np. wykorzystanie większej liczby klastrów lub innego algorytmu grupującego, użycie innych zmiennych (tj. pozostałych kanałów spektralnych, obliczenie wybranych wskaźników spektralnych czy wysokości obiektów z modelu wysokościowego). Istotna jest także kwestia odpowiedniego doboru tych zmiennych oraz redukcji wymiarowości w przypadku skorelowanych zmiennych.
• Klasyfikacja na większych lub nowych obszarach – nasza analiza obejmowała niewielki obszar, dzięki czemu mogliśmy załadować wszystkie dane do pamięci. Zwykle jednak model trenowany jest na reprezentatywnej próbie, a następnie stosowany do całego obszaru. Niestety, funkcja kmeans() nie posiada metody predykcji, więc musimy napisać ją sami. Nie jest to skomplikowane; wymaga jedynie obliczenia odległości pikseli dla każdego klastra i wybrania najbliższego. Dodatkowo, aby przyspieszyć predykcję, warto byłoby zrównoleglić ten proces. Finalnie, należy pamiętać, że opracowany model może działać niepoprawnie na innych obszarach.
• Porównanie wyników z klasyfikacją nadzorowaną – istnieje wiele map pokrycia terenu (np. CORINE Land Cover, ESA WorldCover, ESRI Land Cover). Warto porównać te dwa podejścia.