Wstęp

R, tuż obok Pythona, jest jednym z najbardziej popularnych języków programowania używanych do analizy danych. R znalazł szczególne zastosowanie w takich dziedzinach jak ekologia, leśnictwo, hydrologia czy teledetekcja dzięki swojej prostocie, automatyzacji oraz ogromnemu wsparciu przez jego społeczność. Jest również używany w profesjonalnych zastosowaniach przez takie firmy jak Microsoft, Facebook, Twitter czy Google.

Niniejszy poradnik prezentuje w jaki sposób pozyskać najpopularniejsze dane przestrzenne, tj. ortofotomapy i Cyfrowe Modele Wysokościowe z zasobów Głównego Urzędu Geodezji i Kartografii używając języka R i pakietu rgugik. Dodatkowo, zawiera on praktyczne przykłady wykorzystania tych danych.

Jeśli chciałbyś rozszerzyć swoją wiedzę i zdobyć praktyczne umiejętności analizy danych nieprzestrzennych, to zachęcamy do uczestnictwa w bezpłatnym kursie online przygotowanym przez Uniwersytet Warszawski dla początkujących w języku polskim lub przeczytania książki Elementarz programisty: Wstęp do programowania używając R. Więcej na temat analizy danych przestrzennych można natomiast przeczytać po angielsku na stronie Geocomputation with R.

Instalacja

Niezbędnym narzędziem w tym poradniku będzie interpreter R, który można pobrać na system:

Dodatkowo, pomocnym narzędziem jest edytor kodu - RStudio, który można pobrać z tego miejsca dla różnych systemów operacyjnych. Oba te programy są całkowicie darmowe, a ich instalacja przebiega w standardowy sposób.

Pierwsze uruchomienie

Uruchamiamy RStudio i wybieramy File > New Project… > New Directory > New Project. W ostatnim oknie podajemy nazwę nowego folderu (np. dane_gugik) oraz jego lokalizację (np. Moje Dokumenty). Następnie klikamy File > New File > R Script (lub kombinacja klawiszy “CTRL + SHIFT + N”). Otworzy się nowe puste okienko i możemy przystąpić do pisania naszego kodu.

Na początku musimy zainstalować potrzebne pakiety. Do instalacji nowych pakietów służy funkcja install.packages().

pakiety = c("sf", "raster", "remotes")
install.packages(pakiety)

Do zmiennej pakiety przypisaliśmy trzy nazwy pakietów (zbiorów dodatkowych funkcji) w R: sf obsługuje dane wektorowe, raster dane rastrowe, a pakiet remotes umożliwi pobranie nam wesji testowej pakietu rgugik.

remotes::install_github("kadyb/rgugik")

Powyższa linia kodu pobierze najnowszą wersję pakietu rgugik, który umożliwia przeszukiwanie i pobieranie danych z zasobów GUGiK.

Wczytanie danych

Zacznijmy od pobrania i wczytania danych bez używania pakietu rgugik. Ze strony Generalnej Dyrekcji Ochrony Środowiska pobierzmy dane wektorowe zawierające rezerwaty przyrody w formacie .shp. Rozpakujmy je do podfolderu o tej samej nazwie ("Rezerwaty") w głównym folderze dane_gugik.

W kolejnym kroku załadujmy pakiety, które zainstalowaliśmy wcześniej. Do tego celu służy funkcja library().

library(sf)
library(raster)

Teraz wczytajmy nasze dane przy pomocy funkcji read_sf(), która jako argument przyjmuje ścieżkę do pliku. Oprócz tego, musimy ustawić odpowiednie kodowanie polskich znaków (WINDOWS-1250) jako argument tej funkcji.

rezerwaty = read_sf("Rezerwaty/RezerwatyPolygon.shp",
                    options = "ENCODING=WINDOWS-1250")

Sprawdźmy jak wygląda struktura naszych danych:

head(rezerwaty)
gid nazwa kodinspire geometry
2509 Jezioro Łuknajno PL.ZIPOP.1393.RP.1113 MULTIPOLYGON (((674824.2 66…
2510 Stary Czapliniec PL.ZIPOP.1393.RP.28 MULTIPOLYGON (((676241.8 66…
2511 Rezerwat Skalny im. Jana Czarnockiego PL.ZIPOP.1393.RP.1031 MULTIPOLYGON (((611690.8 33…
2512 Jata PL.ZIPOP.1393.RP.1126 MULTIPOLYGON (((721593.2 45…
2513 Dębina PL.ZIPOP.1393.RP.30 MULTIPOLYGON (((661007.4 50…
2514 Jedlina PL.ZIPOP.1393.RP.14 MULTIPOLYGON (((684158.7 47…

Jest to tabela (ramka danych) składająca się z czterech kolumn (gid, nazwa, kodinspire oraz geometry) i 1674 obiektów (rezerwatów przyrody).

Na potrzebny analizy wybierzmy sobie jeden - rezerwat Krajkowo położony w województwie wielkopolskim.

# selekcja poprzez atrybut
krajkowo = rezerwaty[rezerwaty$nazwa == "Krajkowo", ]

Wyświetlić geometrię poligonu możemy na dwa sposoby:

  1. Wykorzystując funkcję plot i bezpośrednio podając kolumnę z geometrią obiektu: plot(krajkowo$geometry)
  2. Wykorzystując funkcje plot oraz st_geometry, która pozyska geometrię z poligonu. W przypadku pierwszym musimy znać nazwę kolumny z geometrią, natomiast w drugim przypadku geometria zostanie przekazana automatycznie (jest to bezpieczniejszy sposób).
# najpierw pobieramy geometrię obiektu,
# a następnie wyświetlamy ją (funkcja zagnieżdżona)
plot(st_geometry(krajkowo))

Ortofotomapa

Pobieranie

Do pobrania danych ortofoto użyjemy pakietu rgugik. Załadujemy go do R identycznie jak poprzednie pakiety używając funkcji library().

library(rgugik)

Następnie wykorzystamy funkcję orto_request(), która pokaże nam jakie są dostępne ortoobrazy na analizowanym obszarze. Jako argument tej funkcji podamy nasz poligon Krajkowo (obiekt krajkowo).

req_df = orto_request(krajkowo)

Tabelę wynikową z dostępnymi ortoobrazami możemy wyświetlić używając poniższego kodu1.

# wyświetlamy 10 pierwszych wyników oraz 6 pierwszych kolumn
req_df[1:10, 1:6]
sheetID year resolution composition sensor CRS
N-33-142-B-d-4-2 2004 0.50 B/W Analog PL-1992
N-33-142-B-d-4-4 2004 0.50 B/W Analog PL-1992
N-33-142-B-d-4-4 2010 0.25 RGB Digital PL-1992
N-33-142-B-d-4-2 2010 0.25 RGB Digital PL-1992
N-33-142-B-d-4-4 2010 0.25 CIR Digital PL-1992
N-33-142-B-d-4-2 2010 0.25 CIR Digital PL-1992
N-33-142-B-d-4-2 2016 0.25 CIR Digital PL-1992
N-33-142-B-d-4-2 2016 0.25 RGB Digital PL-1992
N-33-142-B-d-4-4 2016 0.25 CIR Digital PL-1992
N-33-142-B-d-4-4 2016 0.25 RGB Digital PL-1992

W ramach ćwiczeń pobierzemy kompozycję w podczerwieni obrazującą najbardziej aktualne dane. Mamy już naszą tabelę z dostępnymi ortoobrazami (req_df), więc w kolejnym kroku możemy dokonać selekcji poprzez wybór tych wierszy, dla których kolumna composition przyjmuje wartość "CIR".

# wybieramy obrazy w podczerwieni i nadpisujemy obiekt req_df
req_df = req_df[req_df$composition == "CIR", ]

Następnie posortujmy tabelę według roku wykonania zdjęcia, gdzie na początku znajdą się te najaktualniejsze zobrazowania.

# minus oznacza kolejność malejącą
req_df = req_df[order(-req_df$year), ]

Wyświetlmy tabelę ponownie i wybierzmy najnowsze kompozycje.

req_df[, 1:6]
sheetID year resolution composition sensor CRS
7 N-33-142-B-d-4-2 2016 0.25 CIR Digital PL-1992
9 N-33-142-B-d-4-4 2016 0.25 CIR Digital PL-1992
22 N-33-142-B-d-4-2 2013 0.25 CIR Digital PL-1992
24 N-33-142-B-d-4-4 2013 0.25 CIR Digital PL-1992
5 N-33-142-B-d-4-4 2010 0.25 CIR Digital PL-1992
6 N-33-142-B-d-4-2 2010 0.25 CIR Digital PL-1992
req_df = req_df[req_df$year == 2016, ]

Zwróćmy uwagę, że w wyniku otrzymaliśmy dwa obiekty (obrazy). Oznacza to, że nasz rezerwat Krajkowo jest zobrazowany na dwóch zdjęciach w ramach jednej serii. Atrybut seriesID służy więc do łączenia mniejszych obrazów w większe mozaiki.

# wybór kolumn (atrybutów) poprzez ich indeks
# 1: "sheetID", 2: "year", 3: "resolution", itd.
# zostały wybrane kolumny od 1 do 6 oraz 9
req_df[, c(1:6, 9)]
sheetID year resolution composition sensor CRS seriesID
7 N-33-142-B-d-4-2 2016 0.25 CIR Digital PL-1992 69837
9 N-33-142-B-d-4-4 2016 0.25 CIR Digital PL-1992 69837

Do pobrania ortoobrazow z serwerów GUGiK służy funkcja tile_download(), która jako główny argument przyjmuje naszą wyselekcjonowaną tabelę.

Pobierane dane domyślnie zapisywane są do aktualnego katalogu roboczego2.

tile_download(req_df)

Przetwarzanie

Wczytajmy pobrane ortoobrazy za pomocą funkcji brick() z pakietu raster, która umożliwia ona pracę na rastrach składających się z kilku kanałów - w naszym przypadku trzech: NIR, R, G.

img1 = brick("69837_329609_N-33-142-B-d-4-2.TIF")
img2 = brick("69837_329613_N-33-142-B-d-4-4.TIF")

Następnie przeprowadzimy operacje: połączenia, docięcia do obszaru rezerwatu i zamaskowania pikseli poza zasięgiem poligonu. Ten proces może zająć kilka minut.

# łączenie rastrów
img = merge(img1, img2)
# docięcie rastrów do poligonu
img = crop(img, krajkowo)
# zamaskowanie pikseli poza poligonem
img = mask(img, krajkowo)

Wyświetlmy efekt końcowy używając funkcji plotRGB(). Funkcja plotRGB() tworzy wizualizację poprzez połączenie wartości trzech warstw rastrowych. Poniższa kompozycja przedstawiona jest w podczerwieni, a nie w barwach naturalnych, co może mylnie interpretować z nazwy funkcji.

plotRGB(img)

W ostatnim kroku obliczmy wskaźnik znormalizowany różnicowy wskaźnik wegetacji NDVI wykorzystując kanały podczerwieni (1) i czerwieni (2).

ndvi = (img[[1]] - img[[2]]) / (img[[1]] + img[[2]])
plot(ndvi, main = "NDVI")

Zadziwiającym spostrzeżeniem są względnie niskie wartości NDVI dla obszaru leśnego. Przyczyny są dwie, tj. zdjęcia wykonywane są w okresie jesiennym oraz prawdopodobnie nie zostały skalibrowane. Z tego powodu lepszym źródłem danych do analiz mogą być zdjęcia satelitarne, które podlegają kalibracji spektralnej i stanowią ciągły zasób danych (jeśli nie wystąpi zachmurzenie).

Cyfrowe Modele Wysokościowe

Pobieranie

W tej sekcji przeprowadzimy analogiczną procedurę, ale tym razem pobierzemy dane wysokościowe dla tego samego obszaru. Funkcja, która nam do tego posłuży to DEM_request() działająca identycznie jak orto_request().

req_df = DEM_request(krajkowo)

Do analizy wykorzystamy Numeryczny Model Terenu (DTM) oraz Numeryczny Model Pokrycia Terenu (DSM). Posłużą one do obliczenia wysokości drzew i budynków. Aby ułatwić selekcje danych, stworzymy dwie tabele, a na końcu je złączymy.

req_df_DTM = req_df[req_df$format == "ARC/INFO ASCII GRID" &
                    req_df$product == "DTM" &
                    req_df$year == 2011, ]
req_df_DSM = req_df[req_df$format == "ARC/INFO ASCII GRID" &
                    req_df$product == "DSM" &
                    req_df$year == 2011, ]
# łączenie tabel
req_df = rbind(req_df_DTM, req_df_DSM)
req_df[, 1:6]
sheetID year format resolution avgElevErr CRS
3 N-33-142-B-d-4-2 2011 ARC/INFO ASCII GRID 1.0 m 0.15 PL-1992
4 N-33-142-B-d-4-4 2011 ARC/INFO ASCII GRID 1.0 m 0.15 PL-1992
7 N-33-142-B-d-4-2 2011 ARC/INFO ASCII GRID 1.0 m 0.15 PL-1992
8 N-33-142-B-d-4-4 2011 ARC/INFO ASCII GRID 1.0 m 0.15 PL-1992

Pobierzmy dane ponownie używając tile_download().

tile_download(req_df)

Przetwarzanie

Wczytajmy je za pomocą funkcji raster() (w odróżnieniu od funkcji brick() wczytuje tylko rastry składające się z jednego kanału).

# DTM
DTM1 = raster("4072_8123_N-33-142-B-d-4-2.asc")
DTM2 = raster("4072_8125_N-33-142-B-d-4-4.asc")

# DSM
DSM1 = raster("3182_10015_N-33-142-B-d-4-2.asc")
DSM2 = raster("3182_10017_N-33-142-B-d-4-4.asc")

Połączmy nasze rastry.

DTM = merge(DTM1, DTM2)
DSM = merge(DSM1, DSM2)

Nałóżmy je na siebie w kolejności NMPT (DSM), NMP (DTM).

DEM = brick(DSM, DTM)
# nadajmy nazwy naszym warstwom
names(DEM) = c("DSM", "DTM")

Dotnijmy i zamaskujmy do zasięgu poligonu.

DEM = crop(DEM, krajkowo)
DEM = mask(DEM, krajkowo)

I wyświetlmy.

plot(DEM)

Teraz możemy obliczyć wysokość drzew i budynków.

diffDEM = DEM[[1]] - DEM[[2]]

Ostatnim krokiem jest wizualizacja wyniku różnicy.

plot(diffDEM, main = "Wysokość obiektów [m]")

Informacje o rgugik

Strona projektu: https://kadyb.github.io/rgugik

GitHub: https://github.com/kadyb/rgugik/

Obsługiwane zbiory danych:

  • ortofotomapa,
  • Baza Danych Obiektów Ogólnogeograficznych (BDOO),
  • Państwowy Rejestr Nazw Geograficznych (PRNG),
  • geometria działek katastralnych,
  • modele 3D budynków (LOD1, LOD2),
  • Cyfrowe Modele Wysokościowe (NMT, NMPT, chmura punktów).

  1. Oprócz wymienionego sposobu pomocna jest również funkcja View().↩︎

  2. Jego miejsce można sprawdzić funkcją getwd()↩︎