rgugik: pobieranie danych z GUGiK w R
Krzysztof Dyba, Jakub Nowosad
2020-08-23
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:
- Wykorzystując funkcję
ploti bezpośrednio podając kolumnę z geometrią obiektu:plot(krajkowo$geometry) - Wykorzystując funkcje
plotorazst_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).