library("terra")

Generowanie danych

Na poprzednich zajęciach wykorzystaliśmy funkcję rast() do wczytania danych rastrowych, niemniej posiada ona więcej zastosowań. Na przykład można użyć ją do konwersji (macierz -> raster) czy wygenerowania nowego rastra z zadanymi parametrami. W przypadku tworzenia nowego rastra można zdefiniować m. in. liczbę wierszy, kolumn i warstw, zakres przestrzennych, układ współrzędnych czy rozdzielczość. Zdefiniujmy “szablon” bez wartości:

raster = rast(nrows = 10, ncols = 20, nlyrs = 1,
              xmin = -180, xmax = 180, ymin = -90, ymax = 90,
              crs = "EPSG:4326")
raster
## class       : SpatRaster 
## dimensions  : 10, 20, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)

Funkcja rast() przyjmuje argument vals, który zostawiliśmy pusty, w związku z czym stworzonemu rastrowi nie zostały przypisane żadne wartości. Możemy to potwierdzić za pomocą funkcji hasValues().

hasValues(raster)
## [1] FALSE

De facto, oznacza to, że stworzyliśmy jedynie metadane rastra (tj. “szablon”). Jeśli chcemy go wypełnić, to musimy przypisać / wygenerować tyle wartości z ilu komórek składa się raster.

n = ncell(raster) # liczba komórek rastra
wartosci = rnorm(n, mean = 0, sd = 1) # wylosuj wartości z rozkładu normalnego
head(wartosci) # wyświetl 6 pierwszych wylosowanych wartości
## [1] -0.4706121 -0.5787918 -0.5195526  1.0166311  3.1263395 -1.5354948
# hist(wartosci) # histogram

Teraz przypiszmy wygenerowane wartości do raster za pomocą funkcji values().

values(raster) = wartosci
raster
## class       : SpatRaster 
## dimensions  : 10, 20, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     lyr.1 
## min value   : -3.504558 
## max value   :  3.126339

Wyświetlmy obiekt raster. Dodatkowo możemy wyświetlić wartości komórek jako etykiety używając funkcji text() z argumentami określającymi liczbę znaków (digits) oraz rozmiar tekstu (cex).

plot(raster)
text(raster, digits = 1, cex = 0.7)

Algebra rastrów

Na rastrach można wykonywać standardowe działania algebraiczne, operacje logiczne oraz funkcje matematyczne.

r1 = sqrt(raster + 10)
r1
## class       : SpatRaster 
## dimensions  : 10, 20, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :    lyr.1 
## min value   : 2.548616 
## max value   : 3.623029
plot(r1)

r2 = abs(raster - r1)
r2
## class       : SpatRaster 
## dimensions  : 10, 20, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     lyr.1 
## min value   : 0.4966896 
## max value   : 6.0531736
plot(r2)

r3 = raster < 0 # operacja logiczna
r3
## class       : SpatRaster 
## dimensions  : 10, 20, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        : lyr.1 
## min value   : FALSE 
## max value   :  TRUE
plot(r3)

Operacje

Docinanie

Nadmierny zasięg rastra wykraczający poza obszar analizy można dociąć za pomocą funkcji crop() używając innego rastra (SpatRaster), wektora (SpatVector) lub zasięgu zdefiniowanego przy użyciu współrzędnych (SpatExtent). W poniższym przykładzie wyznaczymy zasięg używając funkcji ext().

zasieg = ext(-100, 100, -50, 50) # xmin, xmax, ymin, ymax
r = crop(raster, zasieg)
r
## class       : SpatRaster 
## dimensions  : 6, 12, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -108, 108, -54, 54  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     lyr.1 
## min value   : -2.090382 
## max value   :  2.485482
plot(r, ext = ext(raster)) # porównaj z oryginalnym zasięgiem

Należy zauważyć, że zasięg dociętego rastra może różnić się od oczekiwanego, ponieważ zwracane są pełne komórki w wierszach i kolumnach.

Zmiana rozdzielczości

Rozdzielczość przestrzenną rastra można zmienić za pomocą trzech funkcji:

  • disagg() (zwiększa rozdzielczość, czyli komórki stają się mniejsze)
  • aggregate() (zmniejsza rozdzielczość, czyli komórki stają się większe)
  • resample() (przepróbkowanie do zdefiniowanej siatki)

Pierwsze dwie funkcje wymagają podanie współczynnika agregacji, natomiast ostatnia wymieniona funkcja wymaga wskazania rastra z oczekiwaną geometrią. W przypadku funkcji disagg() dostępne są dwie metody, tj. interpolacji najbliższego sąsiada (near) oraz dwuliniowej (bilinear). Ta pierwsza stosowana jest najczęściej w przypadku danych kategorycznych. Zauważ również, że zastosowanie tej metody powoduje podział komórki na mniejsze części, co nie wpływa na efekt wizualizacji.

r1 = disagg(raster, fact = 2, method = "bilinear")
r2 = aggregate(raster, fact = 2, fun = "mean")
par(mfrow = c(1, 3)) # wyświetl 3 rastry obok siebie
plot(raster, main = "Raster wejściowy")
plot(r1, main = "Upsampling")
plot(r2, main = "Downsampling")

Rozdzielczość rastrów można także sprawdzić za pomocą funkcji res(). Sprawdźmy jeszcze jak zastosować funkcję resample() w praktyce.

szablon = rast(nrows = 20, ncols = 40,
               xmin = -180, xmax = 180, ymin = -90, ymax = 90)
r3 = resample(raster, szablon, method = "bilinear")

Reprojekcja

Do transformacji przestrzennego układu współrzędnych służy funkcja project(). Tak jak w przypadku poprzednich funkcji można wykorzystać różne metody interpolacji. Przetransformujmy nasz aktualny układ EPSG:4326 wyrażony w stopniach do odwzorowania Mollweidego (+proj=moll) wyrażonego w metrach.

r_moll = project(raster, "+proj=moll", method = "bilinear")
r_moll
## class       : SpatRaster 
## dimensions  : 22, 45, 1  (nrow, ncol, nlyr)
## resolution  : 788849.2, 806777.6  (x, y)
## extent      : -18040096, 17458119, -8729059, 9020048  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=moll +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs 
## source(s)   : memory
## name        :     lyr.1 
## min value   : -2.761979 
## max value   :  1.889298

Porównajmy otrzymany raster z oryginalnym.

raster
## class       : SpatRaster 
## dimensions  : 10, 20, 1  (nrow, ncol, nlyr)
## resolution  : 18, 18  (x, y)
## extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326) 
## source(s)   : memory
## name        :     lyr.1 
## min value   : -3.504558 
## max value   :  3.126339
par(mfrow = c(1, 2))
plot(raster, main = "World Geodetic System 84")
plot(r_moll, main = "Odwzorowanie Mollweidego")

Jeśli posiadasz dane rastrowe i wektorowe w różnych układach współrzędnych, to zalecana jest reprojekcja danych wektorowych z dwóch powodów:

  • jest szybsza
  • jest odwracalna (transformacja odbywa się bez utraty precyzji)

Statystyki globalne

Wyliczenie statystyk komórek warstwy rastrowej można wykonać za pomocą funkcji global(). Dla przykładu:

# data.frame() łączy poniższe obiekty do jednego
data.frame(
  global(raster, "sum"),
  global(raster, "mean"),
  global(raster, "sd")
)
##            sum       mean       sd
## lyr.1 18.61371 0.09306856 0.997373

Statystyki strefowe

Jeśli posiadamy dwa rastry, tj. pierwszy numeryczny, a drugi kategoryczny (określający strefy), to możemy wyliczyć statystyki strefowe. Najpierw wykorzystajmy funkcje sample() do wygenerowania kategorii oznaczonych numerami od 1 do 3 (w rzeczywistej analizie mogą one reprezentować, np. zbiornik wodny, las i strefę zabudowaną). Musimy ustawić argument replace = TRUE, aby wykonać losowanie ze zwracaniem.

kategorie = sample(1:3, size = ncell(raster), replace = TRUE)
head(kategorie)
## [1] 3 2 2 3 3 2

Następnie kopiujemy raster wejściowy raster wykorzystując funkcję rast zastępując przy tym jego wartości vals = kategorie i nowy raster nazywamy strefy.

strefy = rast(raster, vals = kategorie)
plot(strefy, main = "Strefy", col = c("blue", "green", "red"))

Teraz możemy obliczyć statystyki strefowe.

stat_strefy = zonal(raster, strefy, fun = "mean")
stat_strefy
##   lyr.1     lyr.1.1
## 1     1  0.10082205
## 2     2 -0.03337859
## 3     3  0.21434049

Reklasyfikacja

Reklasyfikacja wartości rastra możliwa jest za pomocą funkcji classify(). Załóżmy, że chcemy wszystkie wartości poniżej 0 zamienić na brakujące wartości (NA), a wszystkie wartości powyżej 0 zamienić na 10. W tym celu musimy stworzyć tabelę klasyfikacyjną składającą się z trzech kolumn:

  1. dolny przedział wartości
  2. górny przedział wartości
  3. nowa wartość

Tabelę najprościej stworzyć definiując uprzednio wektor i następnie transformując go do macierzy.

# wektor
tabela = c(
  -Inf, 0, NA,
  0, Inf, 10
)
# zamiana wektora na macierz
tabela = matrix(tabela, ncol = 3, byrow = TRUE)
tabela
##      [,1] [,2] [,3]
## [1,] -Inf    0   NA
## [2,]    0  Inf   10

Teraz przeprowadźmy reklasyfikację.

reklasyfikacja = classify(raster, tabela)
plot(reklasyfikacja)

Istnieje również prostsza alternatywa w postaci subst() ukierunkowana na zastępowanie wartości.

Okno ruchome

Ostatnią omawianą funkcją jest focal() służącą do obliczania wartości na podstawie sąsiadujących komórek rastra w ruchomym oknie. Można zdefiniować dowolny kształt sąsiedztwa (np. prostokąt, koło) oraz nadać komórkom odpowiednie wagi. Zazwyczaj niniejszą funkcję wykorzystuje się do wygładzania danych numerycznych za pomocą średniej czy mediany oraz danych kategorycznych za pomocą mody.

r = focal(raster, w = 3, fun = "mean")
par(mfrow = c(1, 2))
plot(raster, main = "Raster wejściowy")
plot(r, main = "Raster wygładzony")

Omówiliśmy tylko podstawowe i najczęściej stosowane funkcje do przetwarzania danych rastrowych. Pakiet terra oferuje ich zdecydowanie więcej (sprawdź dokumentację).

Zadanie

3. Z portalu OpenTopography pobierz numeryczny model pokrycia terenu (np. NASADEM lub SRTM) dla małego fragmentu Polski i wykonaj następujące czynności:

  • sprawdź metadane pobranego rastra i wyświetl go
  • transformuj układ współrzędnych do EPSG:2180
  • przytnij zasięg do mniejszego obszaru wyznaczając go przy pomocy funkcji draw(x = "extent") (wymaga wskazania dwóch wierzchołków)
  • wygładź wartości używając średniej arytmetycznej w oknie ruchomym
  • oblicz statystyki rastra, tj. wartość minimalna i maksymalna, średnia oraz odchylenie standardowe
  • zapisz raster na dysku z kompresją LZW