Statystyki strefowe (Zonal statistics) w R

Statystyki strefowe to pomocne narzędzie pozwalające na wyznaczenie dowolnej statystyki z pikseli rastra pokrywających się np. z poligonem z warstwy wektorowej. Najczęściej w oprogramowaniu GIS mamy domyślnie do wyboru kilkanaście statystyk np. maksymalna, średnia, minimalna wartość itp. A co w przypadku gdy chcemy obliczyć statystykę wymyśloną przez nas. W oprogramowaniu GIS nie jest to niemożliwe, ale wymaga sporo pracy, np. z wykorzystaniem Modelera. Na szczęście mamy R-a, gdzie można to zrobić w jednej funkcji:) Zacznijmy jednak od prostych rzeczy. Inicjujemy bibliotekę raster w której znajduje się funkcja extract, służąca do ekstrakcji wartości pikseli rastra leżących w danym wektorze:

install.packages(‘raster’)
library(raster)

Inicjujemy bibliotekę rgdal, którą wykorzystamy do wczytania danych wektorowych:

library(rgdal)

Wczytujemy plik rastrowy oraz plik z warstwą wektorową np. z darmowych warstw udostępnianych przez CODGiK (w moim przypadku NMT i warstwa powiatów województwa mazowieckiego):

r = raster(“C:/grid_mazowieckie.tif”)

g = readOGR(“C:/”,”powiaty_maz”)

Wczytane warstwy wyglądają następująco:

Używając funkcji extract z biblioteki raster możemy wyciągnąć wartości pikseli NMT znajdujących się wewnątrz np. dwóch pierwszych powiatów:

h = extract(r,powiaty[1:2,])

Utworzona zmienna h jest listą, zawierającą w poszczególnych elementach wartości wysokości dla konkretnego powiatu. Wartości te wywołujemy używając:

h[[1]]

[1] 98.47 98.97 99.38 99.85 99.68 99.31

W przypadku, gdy dla konkretnego powiatu chcemy policzyć np. średnią wartość wysokości musimy dodać do funkcji extract atrybut fun = mean:

h_mean = extract(r,powiaty[1:2,],fun=mean)

Wywołując w konsoli zmienną h_mean widzimy, że jest to już macierz wartości:

h_mean
         [,1]
[1,] 119.7398
[2,] 119.0318

Jako atrybut fun do funkcji extract możemy wstawić dowolną funkcję utworzoną przez nas. Pamiętamy, że dla każdego z powiatów otrzymujemy wektor wartości i będziemy go musieli przetwarzać. Zdefiniujmy funkcję tworzącą statystykę, która będzie nam zwracała procent pikseli o wartościach większych od średniej wysokości w powiecie. Najpierw utworzymy funkcję nasza_funkcja zawierającą tylko jedną zmienną lokalną x będąco naszym zbiorem wartości pikseli oraz argument logiczny na.rm odpowiedzialny za usuwanie wartości NA:

nasza_funkcja function(x,na.rm=TRUE)

{

Dodajmy do funkcji zmienną zawierającą wartość średnią

avg = mean(x)

Policzmy ile jest wszystkich pikseli:

all = length(x)

i policzmy ile jest tych większych od avg:

more_than_avg = sum(x › avg)

policzmy procent:

prc = (more_than_avg/all)*100

Funkcja będzie zwracała procent:

return(prc)

Zamykamy definicję funkcji:

}

Nasza funkcja ostatecznie powinna wyglądać następująco:

nasza_funkcja function(x,na.rm=TRUE)

{

avg = mean(x)

all = length(x)

more_than_avg = sum(x›avg)

prc = (more_than_avg/all)*100

return(prc)

}

W celu wyznaczenia wartości procentów wstawiamy ją do funkcji extract do atrybutu fun:

pix_prc = extract(r,powiaty[1:2,],fun=nasza_funkcja)

W wyniku otrzymujemy macierz zawierającą procenty ze zdefiniowanej przez nas funkcji:

prc_pixeli

         [,1]
[1,] 46.52973
[2,] 38.13020

Funkcję extract możemy modyfikować w zależności jakie dane chcemy uzyskać. Polecam zajrzeć do dokumentacji technicznej tej funkcji, lub do pomocy:

help(“extract”)