Reklasyfikacja rastra w R

Często w przetwarzaniu obrazów chcemy przypisać konkretnemu zakresowi wartości rastra jedną niepowtarzalną liczbę, na przykład przy tworzeniu różnego rodzaju masek. W teledetekcji taka metoda zamiany wartości pikseli na klasy nazywa się klasyfikacją progową. Przykładem takiego opracowania jest maska wód z zobrazowania satelitarnego, gdzie piksele reprezentujące zbiorniki przyjmują wartość 1. Pozostałe mają wartości 0.

Jak to zrobić w R? Nic prostszego. Najpierw wczytujemy zobrazowanie korzystając z funkcji raster w bibliotece raster. Ja skorzystam z zobrazowanie satelity Sentinel pozyskanego w kanale 8 (podczerwonym). Wyświetlamy wczytane zobrazowanie:

plot(b8,col=grey.colors(255))

W programie GISowym zwektoryzowałem wcześniej cztery obszary, gdzie występują zbiorniki wodne. Wykorzystam tą warstwę do określenia wartości progowej dla maski. Wczytujemy utworzoną warstwę używając funkcji readOGR. Korzystając z funkcji extract odczytujemy wartości pikseli dla obszarów ze zbiornikami wodnymi i zapisujemy je do zmiennej ext:

ext = extract(b8, obszary)

Zmienna ext zawiera listę czterech wektorów z wartościami pikseli odczytanych dla poszczególnych obszarów:

ext

[[1]]

  [1] 274 274 276 273 266 267 269 266 275 260 273 270 266 267 266 270 272

 [18] 258 262 266 259 259 253 260 274 267 267 267 268 265 266 264 258 265

…

[[2]]

  [1] 264 262 261 256 258 264 261 256 256 260 262 255 258 256 256 256 256

 [18] 258 258 259 261 257 257 256 257 256 256 256 256 256 256 257 259 256

…

[[3]]

…

[[4]]

…

Funkcja unlist zamieni listę na jeden wektor zawierający wartości pikseli dla wszystkich obszarów:

ext = unlist(ext)

ext

[1] 290 296 290 288 290 286 281 276 274 276 283 285 282 281 274 268

  [17] 260 268 276 275 282 281 277 270 269 264 263 264 273 270 270 275

…

Policzmy średnią z wartości w ext:

avg = mean(ext)

Wartość średnia będzie naszym progiem. Reklasyfikujmy piksele obrazu zgodnie z założeniem, że wszystkie wartości poniżej średniej będą miały wartość 1. Pozostałe przyjmują wartość 0. Skorzystamy z funkcji reclassify. Dla funkcji tej musimy zdefiniować raster, który będziemy reklasyfikować oraz wektor zmiany wartości, który ma postać:

c(próg_lewy, próg_prawy, nowa_wartość, próg_lewy, próg_prawy, nowa_wartość, …)

Progami może być wartość nieskończona definiowana jako -Inf i Inf. Ostatecznie naszą maskę tworzymy używając wyrażenia:

mask = reclassify(b8,c(-Inf,avg,1,avg,Inf,0))

Wyświetlamy wynik:

plot(mask)

Otrzymaliśmy maskę zbiorników wodnych na podstawie zobrazowania satelitarnego wykonanego w bliskiej podczerwieni. Narzędziem możemy tworzyć dowolne maski różnego rodzaju danych rastrowych.