library("terra")

Wstęp

Wczytajmy dane wektorowe Luksemburgu wykorzystane na pierwszych zajęciach.

sciezka = system.file("ex/lux.shp", package = "terra")
wektor = vect(sciezka)
wektor
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 12, 6  (geometries, attributes)
##  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
##  source      : lux.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA       POP
##  type        : <num>    <chr> <num>    <chr> <num>     <num>
##  values      :     1 Diekirch     1 Clervaux   312 1.808e+04
##                    1 Diekirch     2 Diekirch   218 3.254e+04
##                    1 Diekirch     3  Redange   259 1.866e+04

Warto zauważyć, że funkcja vect() poza standardowym odczytem oferuje więcej zaawansowanych możliwości:

  • wczytanie określonego zakresu przestrzennego
  • wczytanie określonych warstw
  • wczytanie określonego rodzaju danych (geometrie lub atrybuty)
  • tworzenie zapytań bazodanowych (SQL)
  • stworzenie wskaźnika do pliku (proxy)

Zasadniczo, do atrybutów warstwy wektorowej możemy odwołać się na dwa sposoby używając:

  • znaku dolara $ – zwraca wektor bez geometrii
  • nawiasu kwadratowego [] – zwraca geometrię z wybranymi atrybutami
wektor$NAME_2
##  [1] "Clervaux"         "Diekirch"         "Redange"          "Vianden"         
##  [5] "Wiltz"            "Echternach"       "Remich"           "Grevenmacher"    
##  [9] "Capellen"         "Esch-sur-Alzette" "Luxembourg"       "Mersch"
wektor[, "NAME_2"]
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 12, 1  (geometries, attributes)
##  extent      : 5.74414, 6.528252, 49.44781, 50.18162  (xmin, xmax, ymin, ymax)
##  source      : lux.shp
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :   NAME_2
##  type        :    <chr>
##  values      : Clervaux
##                Diekirch
##                 Redange

Oprócz tego, możemy pozyskać wszystkie atrybuty w postaci ramki danych używając funkcji as.data.frame().

Operacje

Obliczanie powierzchni

W celu obliczenia powierzchni poligonów należy zastosować funkcję expanse(). Ta funkcja automatycznie wykonuje obliczenia dla układów wyrażonych w stopniach, więc nie jest wymagana projekcja do planarnego układu współrzędnych. Można również określić jednostki, w których zwrócone będą wyniki przy pomocy argumentu unit.

powierzchnia = expanse(wektor, unit = "km")
data.frame(nazwa = wektor$NAME_2, powierzchnia)
##               nazwa powierzchnia
## 1          Clervaux    312.28321
## 2          Diekirch    218.67403
## 3           Redange    259.45481
## 4           Vianden     76.20041
## 5             Wiltz    263.17426
## 6        Echternach    188.28214
## 7            Remich    128.99150
## 8      Grevenmacher    210.35449
## 9          Capellen    185.63077
## 10 Esch-sur-Alzette    251.32202
## 11       Luxembourg    237.11300
## 12           Mersch    233.32996

Generowanie punktów

Można wygenerować punkty o rozkładzie regularnym (method = "regular") lub losowych (method = "random") na podstawie wejściowej geometrii używając funkcji spatSample(). Istnieje również możliwość próbkowania stratyfikowanego (wtedy dla każdego poligonu zostanie wygenerowanych \(n\) punktów).

proba = spatSample(wektor, size = 100, method = "random")
plot(wektor)
plot(proba, add = TRUE)

Generowanie otoczki wypukłej

Otoczka wypukła jest najmniejszym wielokątem wypukłym ograniczającym dany zbiór punktów. Innymi słowy, jest to wielokąt, którego wierzchołki stanowią najbardziej zewnętrzne punkty zbioru. Do jej wygenerowania służy funkcja convHull().

otoczka = convHull(wektor)
plot(wektor)
plot(otoczka, add = TRUE, border = "red")

Generowanie buforów

Bufory można wygenerować wykorzystując funkcje buffer(). Bufory obliczane są dla każdej geometrii osobno, w związku z czym, jeśli chcemy stworzyć jeden bufor musimy zastosować funkcję agregującą geometrie, tj. aggregate(); łączy ona wiele geometrii w jedną. Funkcja buffer() wymaga wskazania odległości bufora (argument width) i należy odnotować, że:

  • domyślną jednostką długości dla układu geograficznego są metry (nie stopnie)
  • argument width jest zwektoryzowany, zatem można określić różne odległości dla kolejnych geometrii
bufor = buffer(aggregate(wektor), width = 1000)
plot(wektor)
plot(bufor, add = TRUE, border = "red")

Generowanie centroidów

Centroid jest to punkt określający geometryczny środek wielokąta. Dla wielokątów wypukłych (tj. kąty takiej figury są mniejsze niż 180°) wyliczany jest na podstawie średniej arytmetycznej współrzędnych wierzchołków. Jego wyznaczenie jest możliwe używając funkcji centroids(). W przypadku wielokątów wklęsłych centroid może znajdować się poza obiektem, wtedy może zastosować argument inside = TRUE, który wymusi przybliżoną lokalizację wewnątrz obiektu.

centroidy = centroids(wektor, inside = FALSE)
centroidy_wewnatrz = centroids(wektor, inside = TRUE)
plot(wektor)
plot(centroidy, add = TRUE, col = "blue")
plot(centroidy_wewnatrz, add = TRUE, col = "orange")

Obliczanie odległości

W kolejnym kroku możemy obliczyć jak oddalone są centroidy od siebie używając funkcję distance(). Wymieniona funkcja również automatycznie wykonuje obliczenia dla układów geograficznych (jednostki w stopniach) i wynik domyślnie zwracany jest w metrach. Jeśli podamy jeden argument w funkcji, to zostaną wyliczone odległości każdego obiektu z każdym.

odleglosci = distance(centroidy)

Jako wynik otrzymujemy obiekt klasy dist, który reprezentuje jako wektor dolną połowę macierzy odległości. Wykorzystując funkcję as.matrix() możemy przetransformować ten obiekt do pełnej macierzy.

macierz_odleglosci = as.matrix(odleglosci)
macierz_odleglosci[1:5, 1:5]
##          1         2        3         4        5
## 1     0.00 24290.360 31366.47 19330.917 16147.66
## 2 24290.36     0.000 18794.53  7486.065 17280.29
## 3 31366.47 18794.533     0.00 24595.137 15579.42
## 4 19330.92  7486.065 24595.14     0.000 17986.54
## 5 16147.66 17280.292 15579.42 17986.538     0.00

Dla ułatwienia możemy również kolumnom i wierszom nadać nazwy jednostek administracyjnych.

colnames(macierz_odleglosci) = rownames(macierz_odleglosci) = centroidy$NAME_2
# View(macierz_odleglosci) # wyświetl

Relacje przestrzenne

Do określenia relacji przestrzennych między obiektami służy funkcja relate(). Dla przykładu możemy sprawdzić czy poligon zawiera (relation = "contains") wylosowany punkt. Jako wynik funkcji zwracana jest macierz z wartościami logicznymi.

proba = spatSample(wektor, size = 5, method = "random")
plot(wektor)
plot(proba, add = TRUE)

relate(wektor, proba, relation = "contains")
##        [,1]  [,2]  [,3]  [,4]  [,5]
##  [1,] FALSE FALSE FALSE FALSE  TRUE
##  [2,] FALSE FALSE FALSE FALSE FALSE
##  [3,]  TRUE FALSE FALSE  TRUE FALSE
##  [4,] FALSE FALSE FALSE FALSE FALSE
##  [5,] FALSE FALSE FALSE FALSE FALSE
##  [6,] FALSE FALSE FALSE FALSE FALSE
##  [7,] FALSE FALSE FALSE FALSE FALSE
##  [8,] FALSE FALSE FALSE FALSE FALSE
##  [9,] FALSE FALSE FALSE FALSE FALSE
## [10,] FALSE  TRUE  TRUE FALSE FALSE
## [11,] FALSE FALSE FALSE FALSE FALSE
## [12,] FALSE FALSE FALSE FALSE FALSE

W wierszach zawarte są jednostki administracyjne, natomiast w kolumnach wylosowane punkty.

Docinanie

Podobnie jak w przypadku danych rastrowych, możemy zmniejszyć zasięg warstwy wektorowej używając funkcję crop().

zasieg = ext(c(5.9, 6.3, 49.6, 49.9))
wektor_dociety = crop(wektor, zasieg)
plot(wektor)
plot(wektor_dociety, add = TRUE, border = "red")

Alternatywnie docięcie można wykonać używając funkcji intersect().

Reprojekcja

Reprojekcja danych wektorowych wygląda identycznie tak jak w przypadku danych rastrowych.

wektor_3857 = project(wektor, "EPSG:3857")

Oczywiście w tej sekcji zostały zaprezentowane jedynie wybrane funkcje do przetwarzania danych wektorowych.

Zadanie

4. Pobierz granice powiatów z Geoportalu, następnie wylicz ich centroidy i wskaż te powiaty, które są położone najbliżej i najdalej. Wskazówki:

  • Odległość obiektu od samego siebie wynosi 0 m i trzeba to wykluczyć.
  • Do znalezienia indeksu o minimalnej lub maksymalnej wartości można wykorzystać funkcje which() z argumentem arr.ind = TRUE.

Zastanów się również jak na wyniki analizy wpływają:

5. Na stronie https://dane.gov.pl znajdziesz numeryczne modele terenu w siatce 100 m dla wszystkich województw. Wybierz dowolne województwo i stwórz dwie mapy przedstawiające średnią oraz odchylenie standardowe wysokości powiatów. Jako rozwiązanie zadania przygotuj powtarzalny skrypt.

Dane udostępnione są w formacie tekstowym. Do ich wczytania możesz wykorzystać funkcję read.table(). Sprawdź jaki jest separator kolumn, separator dziesiętny, nagłówek kolumn oraz typ kolumn. W przypadku dużych zbiorów danych, można przeprowadzić analizę na losowej próbie.