library("terra")
library("rgugik")

Wstęp

Główny Urząd Geodezji i Kartografii jest istotnym źródłem danych przestrzennych dla Polski. Dane można przeglądać i pobrać z Geoportalu lub wykorzystując różne usługi. W otwartych zbiorach danych znajdziemy m. in.:

  • ortofotomapy
  • cyfrowe modele wysokościowe (CMW):
    • numeryczny model terenu (NMT)
    • numeryczny model pokrycia terenu (NMPT)
    • chmury punktów
  • modele 3D budynków
  • Państwowy Rejestr Granic (PRG)
  • Baza Danych Obiektów Topograficznych (BDOT)
  • i inne

Wyszukiwanie i pobieranie wymienionych zbiorów danych umożliwia pakiet rgugik.

Ortofotomapa

Ortofotomapa to rastrowe, ortogonalne i kartometryczne przedstawienie powierzchni terenu powstałe w wyniku cyfrowego przetwarzania zdjęć lotniczych lub satelitarnych. Podczas ortorektyfikacji usuwane są zniekształcenia geometryczne wynikające z rzeźby terenu przy użyciu cyfrowych modeli wysokości. Ortofotomapa posiada georeferencje, co pozwala na określenie współrzędnych geograficznych dla każdej komórki obrazu.

Cechy ortofotomapy:

  • Rozdzielczość przestrzenna – związana z rozmiarem najmniejszego obiektu, który może zostać wykryty przez czujnik i jest określana przez rozmiar komórki obrazu (piksel). Im mniejsza komórka, tym więcej szczegółów reprezentuje. Zbyt duży rozmiar oznacza, że poszczególne obiekty na zdjęciu przestają być rozpoznawalne.
  • Kompozycja – obrazy analogowe są przedstawione w odcieniach szarości, natomiast obrazy cyfrowe mogą składać się z naturalnych kolorów (RGB) lub bliskiej podczerwieni (NIR)

Wyszukiwanie

W katalogu dane znajduje się warstwa poligonowa Lublin.gpkg, która reprezentuje granicę Lublina w układzie EPSG:2180.

lublin = vect("../dane/Lublin.gpkg")
lublin
##  class       : SpatVector 
##  geometry    : polygons 
##  dimensions  : 1, 0  (geometries, attributes)
##  extent      : 741124, 756250.7, 369483.5, 387356.7  (xmin, xmax, ymin, ymax)
##  source      : Lublin.gpkg
##  coord. ref. : ETRF2000-PL / CS92 (EPSG:2180)

Funkcja perim() pozwala obliczyć długość granicy, natomiast funkcja expanse() wskazuje powierzchnię obszaru.

perim(lublin) / 1000 # wynik w km
## [1] 81.4903
expanse(lublin, unit = "km")
## [1] 147.4631

Kolejny etap dotyczy wyszukania dostępnych danych dla naszego obszaru zainteresowania. Dla uproszczenie wyznaczmy centroid. Jego współrzędne można odczytać używając funkcji crds().

punkt = centroids(lublin)
crds(punkt)
##             x      y
## [1,] 748243.3 379235
plot(lublin, main = "Lublin")
plot(punkt, col = "blue", add = TRUE)

Teraz wyszukajmy dostępne dane za pomocą funkcji ortho_request().

dane = ortho_request(punkt)

Możemy wyświetlić część otrzymanej ramki danych lub alternatywnie przeglądać całość używając funkcji View().

# wyświetl 10 pierwszych wierszy i 6 pierwszych kolumn
dane[1:10, 1:6]
##            sheetID year resolution composition  sensor        CRS
## 1    M-34-34-A-c-1 2003       0.25         B/W  Analog    PL-1992
## 2  M-34-34-A-c-1-4 2022       0.25         CIR Digital    PL-1992
## 3  M-34-34-A-c-1-4 2020       0.25         RGB Digital    PL-1992
## 4  M-34-34-A-c-1-4 2019       0.10         RGB Digital    PL-1992
## 5  M-34-34-A-c-1-4 2018       0.10         RGB Digital    PL-1992
## 6      8.151.08.12 2015       0.25         RGB Digital PL-2000:S8
## 7  M-34-34-A-c-1-4 2015       0.25         RGB Digital    PL-1992
## 8  M-34-34-A-c-1-4 2015       0.25         CIR Digital    PL-1992
## 9  M-34-34-A-c-1-4 2013       0.10         RGB Digital    PL-1992
## 10     8.151.08.12 2013       0.10         RGB Digital PL-2000:S8

Standardowo dane możemy filtrować z uwzględnieniem zadanych parametrów.

dane[dane$year > 2016, ]
dane[dane$composition == "CIR", ]

I sortować, np. według aktualności.

# kolejność malejąca (najnowsze dane)
dane[order(-dane$year), ]

Pobieranie

Jako przykład pobierzmy dwie kompozycje tego samego obszaru wykonane w naturalnych barwach i z kanałem bliskiej podczerwieni z 2022 r. (ID: 76746_1143325_M-34-34-A-c-1-4 i 76745_1143326_M-34-34-A-c-1-4).

id = c("76746_1143325_M-34-34-A-c-1-4", "76745_1143326_M-34-34-A-c-1-4")
dane_sel = dane[dane$filename %in% id, ]

Po selekcji potrzebnych danych, można je pobrać wykorzystując funkcję tile_download(). Możliwe jest również wskazanie katalogu, do którego powinny zostać pobrane obrazy (argument outdir).

Zazwyczaj warto zwiększyć domyślną wartość przekroczenia czasu połączenia (timeout) z domyślnych 60 sekund w przypadku dużych plików lub wolnego połączenia.

options(timeout = 600)
tile_download(dane_sel, outdir = "dane")

Do wylistowania pobranych plików służy funkcja list.files(). Należy wskazać jakie pliki chcemy wczytać (pattern = "\\.tif$") i zapobiegawczo zwrócić pełne ścieżki do plików (full.names = TRUE).

pliki = list.files("dane", pattern = "\\.tif$", full.names = TRUE)
pliki
## [1] "dane/76745_1143326_M-34-34-A-c-1-4.tif"
## [2] "dane/76746_1143325_M-34-34-A-c-1-4.tif"

W ostatnim kroku możemy kolejno wczytać rastry i je wyświetlić.

# kompozycja z bliską podczerwienią
r1 = rast(pliki[1])
plot(r1)

# kompozycja w naturalnych barwach
r2 = rast(pliki[2])
plot(r2)

Zadanie

6. Pobierz minimum dwa sąsiadujące ze sobą kafelki ortofotomapy z tej samej serii i połącz je:

  1. do jednego pliku .tiff używając funkcji merge()
  2. do jednego wirtualnego pliku .vrt używając funkcji vrt()

Sprawdź zajmowaną ilość miejsca przez te pliki na dysku wykorzystując funkcję file.size() (wynik zwracany jest w bajtach). Sprawdź również zawartość pliku .vrt (czym on jest w rzeczywistości?). Następnie, zmniejsz rozdzielczość mozaiki do 10 m i zapisz wynik. Jak zmieniła się jakość w stosunku do oryginalnego zdjęcia?