Wprowadzenie

Celem pracy jest próba analizy i wizualizacji danych przestrzennych z wykorzystaniem podstawowych funkcjioprogramowania R na przykładzie danych rastrowych.

Metody

Do analizy wykorzystano numeryczny model terenu reprezentujący wysokość topograficzną terenu. Arkusz został pobrany ze strony https://mapy.geoportal.gov.pl/. Jako obszar badań wybrano teren Skałek Przegorzalskich położonych w zachodniej części Krakowa. Dane pozyskane w 2017 roku.

Wspomniany arkusz (nmt_1) został wczytany w R za pomocą funkcji raster(). W tym celu najpierw został wczytany pakiet do obsługi danych rastrowych:

library(raster)
## Loading required package: sp

Dalej została wczytana warstwa nmt_1:

nmt_1 = raster("C:\\Users\\daria\\Desktop\\Studia EGP\\AiWDP\\Zadanie\\nmt_1.asc")

Podstawowe właściwości warstwy nmt_1 były sprawdzone za pomocą wywołania obiekta:

nmt_1
## class      : RasterLayer 
## dimensions : 2341, 2263, 5297683  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 560372.5, 562635.5, 241940.5, 244281.5  (xmin, xmax, ymin, ymax)
## crs        : NA 
## source     : C:/Users/daria/Desktop/Studia EGP/AiWDP/Zadanie/nmt_1.asc 
## names      : nmt_1

Zgodnie z otrzymanymi informacjami klasa obiektu to raster o rozdzielczości 1 x 1. Układ współrzędnych niejest zdefiniowany (crs: NA). Układ współrzędnych 1992 został nadany w następujący sposób:

crs(nmt_1) = "+init=EPSG:2180"
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj
## = prefer_proj): Discarded datum European_Terrestrial_Reference_System_1989 in
## Proj4 definition

Następnie ponownie wywołujemy obiekt i sprawdzamy jego właściwości:

nmt_1
## class      : RasterLayer 
## dimensions : 2341, 2263, 5297683  (nrow, ncol, ncell)
## resolution : 1, 1  (x, y)
## extent     : 560372.5, 562635.5, 241940.5, 244281.5  (xmin, xmax, ymin, ymax)
## crs        : +proj=tmerc +lat_0=0 +lon_0=19 +k=0.9993 +x_0=500000 +y_0=-5300000 +ellps=GRS80 +units=m +no_defs 
## source     : C:/Users/daria/Desktop/Studia EGP/AiWDP/Zadanie/nmt_1.asc 
## names      : nmt_1

Wizualizacja danych przestrzennych realizowana za pomocą funkcji plot(). Wykres jest przedstawiony w następnym rozdziale. Dalej wyświetlimy podstawowe statystyki badanych danych, w tym wartość maksimalną, wartość minimalną, wartość średnią, odchylenie standardowe oraz różnicę pomiędzy wartością maksymalną a minimalną.

cellStats(nmt_1, max)
## [1] 382.59
cellStats(nmt_1, min)
## [1] 198.94
cellStats(nmt_1, mean)
## [1] 268.7672
cellStats(nmt_1, sd)
## [1] 50.07202
cellStats(nmt_1, max) - cellStats (nmt_1, min)
## [1] 183.65

Korzystając z funkcji crop (), utworzono raster o mniejszym zasięgu niż oryginał:

nmt_1_fragment = crop(nmt_1, c(560000, 561500, 242500, 243000))

Korzystając z funkcji cellStats(), wyświetlimy wartości maksymalną i minimalną dla powstałego fragmentu NMT.

cellStats(nmt_1_fragment, max)
## [1] 354.97
cellStats(nmt_1_fragment, min)
## [1] 216.75

Raster o mniejszej rozdzielczości (2 x 2) został utworzony przy użyciu funkcji aggregate():

nmt_1_10 = aggregate(nmt_1, 2)

Oprócz tego, oryginalny raster został reklasyfikowany na 3 grupy z wykorzystaniem funkcji reclassify():

recl_nmt_1 = reclassify(nmt_1, c(0, 220, 1, 220, 270, 2, 270, 390, 3))

W dalszej części pracy wykorzystujemy funkcję terrain() do obliczenia nachylenia, ekspozycji i szorstkościterenu na badanym fragmencie NMT o mniejszym zasięgu.:

nachylenie = terrain(nmt_1_fragment, "slope")
ekspozycja = terrain(nmt_1_fragment, "aspect")
szorstkosc = terrain(nmt_1_fragment, "roughness")

Korzystając z rastrów nachylenia i ekspozycji, utworzymy mapę cieniowania za pomocą funkcji hillShade():

cieniowanie = hillShade(nachylenie, ekspozycja, 55, 270)

Pakiet RColorBrewer został wczytany w R dla rozszerzenia możliwości wizualizacji danych w zakresie kolorów:

library(RColorBrewer)

Wyniki, wizualizacje i wnioski

W toku pracy badanej warstwie nadano układ współrzędnych 1992, dane pierwotne oraz fragment danych o mniejszym zasięgu zwizualizowano na poniższych wykresach.

plot(nmt_1, main = "NMT wybranego obszaru Krakowa")

Raster o mniejszym zasięgu niż oryginał pokazano na poniższym wykresie, który został wykonany za pomocą palety barwną terrain.colors() (liczba kolorów - 50).

plot(nmt_1_fragment, main = "Fragment NMT", col = terrain.colors(50))

Statystyki obliczone za pomocą oprogramowania pozwoliły na zidentifikowanie najwyższego punktu (382,59 m npm) oraz najniższego punktu (198,94 m npm) badanego obszaru. Obliczono również inne statystyki.

Z przedstawionych wykresów wynika, że najwyższe fragmenty badanego obszaru znajdują się w północno-zachodniej części mapy, natomiast wschodnia i południowa części obszaru charakteryzują się mniejszą wysokością.

W trakcie pracy uzyskano również raster o niższej rozdzielczości (2 x 2 zamiast oryginalnego 1 x 1), który jest przedstawiony na poniższym wykresie.

plot(nmt_1_10, main = "NMT o mniejszej rozdzielczości", col = terrain.colors(50))

Dodatkowo oryginalny raster został reklasyfikowany na trzy grupy według wartości wysokości: od 0 do 220, od 220 do 270 i od 270 do 390. Otrzymana mapa jest pokazana na kolejnych dwóch wykresach (z domyślną skalą barwną oraz ze skalą barwną “Accent” pakietu RColorBrewer, która poprawiła wizualizację reklasyfikowanego rastra).

plot(recl_nmt_1, main = "Reklasyfikowany NMT") 

plot(recl_nmt_1, main = "Reklasyfikowany NMT", col = brewer.pal(3, "Accent"))

Poniższe wykresy przedstawiają różne charakterystyki rzeźby terenu: nachylenie, ekspozycję, szorstkość, a także cieniowanie, które zostały obliczone dla fragmentu NMT o mniejszym zasięgu (560000, 561500, 242500, 243000 zamiast oryginalnego 560372,5, 562635,5, 241940,5, 244281,5).

plot(nachylenie, main = "Nachylenie", col = brewer.pal(9, "YlGn"))

plot(ekspozycja, main = "Ekspozycja")

plot(cieniowanie, main = "Cieniowanie", col = grey(1:100/100))

plot(szorstkosc, main = "Szorstkość", col = brewer.pal(9, "Purples"))

Generalnie warto zauważyć, że analiza danych rastrowych przeprowadzona w pracy pozwoliła na wyznaczenie na badanym obszarze fragmentów terenu o większej i mniejszej wysokości, o większym i mniejszym nachyleniu, a także obliczenie innych charakterystyk rzeźby terenu oraz wskaźników statystycznych. Przeprowadzono też transformację oryginalnych danych rastrowych do celów analizy poprzez reklasyfikację, zmianę rozdzielczości i zasięgu. Korzystanie z różnych palet kolorów pozwoliło na zastosowanie rozmaitych opcji wizualizacji badanych danych. W rezultacie postawiony cel pracy można uznać za pomyślnie zrealizowany.