library("terra")

Wprowadzenie

Dane satelitarne

Dane satelitarne odgrywają kluczową rolę w systemach informacji geograficznej i wykorzystywane są w szerokim zakresie do kartowania, analizy pokrycia terenu, monitorowania zmian w środowisku, zarządzania kryzysowego czy rolnictwa precyzyjnego. Źródłem danych mogą być:

sensory optyczne (pozyskują dane w różnych spektrach fali elektromagnetycznej, np. światło widzialne, bliska podczerwień, podczerwień krótkofalowa czy podczerwień termalna)
sensory radarowe (emitują impulsy radarowe, a następnie mierzą odbitą energię. W przeciwieństwie do sensorów optycznych mogą być wykorzystywane do obrazowania powierzchni przez chmury oraz w nocy)

W kontekście ogromnych zbiorów danych satelitarnych nie sposób pominąć zasadniczych technologii ułatwiających ich przetwarzanie:

SpatioTemporal Asset Catalogs – ustandaryzowany sposób organizowania i opisywania zbiorów danych przestrzennych ułatwiający ich wyszukiwanie oraz dostęp (https://stacspec.org).
Byte serving (range requests) – technika używana w protokole HTTP do przesyłania tylko określonej części pliku z serwera do klienta. Jest to szczególnie przydatne w przypadku dużych plików, gdzie pobieranie całego pliku może być nieefektywne lub niepotrzebne. Ta technika umożliwia wydajny dostęp do wymaganych zakresów danych bez przesyłania niepotrzebnych części.
- Cloud Optimized Geotiff – specjalny format przechowywania danych geoprzestrzennych (takich jak zdjęcia satelitarne) zaprojektowany z myślą o środowiskach chmurowych. Opiera się na standardowym formacie Geotiff, ale jest optymalizowany pod kątem wydajnego dostępu i przetwarzania w chmurze. Oprócz byte serving, stosowane są podglądy (overviews), czyli warstwy rastra w niższej rozdzielczości, umożliwiając jego renderowanie bez konieczności dostępu do danych w pełnej rozdzielczości (https://www.cogeo.org/).
- Virtual File Systems – umożliwia jednolity dostęp do różnych typów przechowywania danych (https://gdal.org/user/virtual_file_systems.html). Jednym z modułów jest vsicurl, który zapewnia bezpośredni dostęp do zdalnych plików przez protokoły HTTP(S) i FTP bez konieczności pełnego pobierania, tym samym minimalizując transfer danych.

Przetwarzanie potokowe

Przetwarzanie potokowe w kontekście analizy danych odnosi się organizacji przepływu pracy w sposób liniowy (sekwencyjny). Oznacza to, że dane przepływają przez serię etapów (funkcji), w których wynik jednego etapu służy jako wejście do kolejnego. W praktyce sprawia to, że kod jest bardziej przejrzysty i łatwiejszy do zrozumienia.

W R wbudowanym operatorem przypływu jest |>. Do jego zapisu można wykorzystać skrót klawiszowy CTRL + SHIFT + M, jednak wymaga to zaznaczenia opcji Use native pipe operator (zakłada Code > Editing) w RStudio.

Przykładowo, możemy wylosować 10 liczb z rozkładu normalnego, następnie obliczyć wartość bezwzględną, posortować w kolejności rosnącej i finalnie nadpisać obiekt dane.

dane = rnorm(10)
dane |> 
  abs() |>
  sort() -> dane
dane

##  [1] 0.1542144 0.1747443 0.1949715 0.2234891 0.2378379 0.7143219 1.0003207
##  [8] 1.0516586 1.3784942 2.0626081

Pozyskiwanie danych

Do obsługi katalogów STAC w R służy pakiet rstac, który umożliwia wyszukiwanie i pobieranie danych zgodnie z tym standardem.

library("rstac")

W pierwszym kroku należy zdefiniować źródło danych używając funkcji stac(). Wykaz źródeł można znaleźć na STAC Index. W naszym przykładzie wykorzystamy usługę Earth Search, która dostarcza m. in. zdjęcia satelitarne z Landsata oraz Sentinela.

stac_source = stac("https://earth-search.aws.element84.com/v1")
stac_source

## ###rstac_query
## - url: https://earth-search.aws.element84.com/v1/
## - params:
## - field(s): version, base_url, endpoint, params, verb, encode

W wyniku tej operacji otrzymaliśmy obiekt rstac_query, który zawiera informacje o zapytaniu HTTP, które zostanie wysłane do serwera. Możemy wyróżnić dwie metody wysyłania żądań:

GET (get_request()) – służy głównie do pobierania danych z serwera. Dane przesyłane są poprzez adres URL jako parametry zapytania.
POST (post_request()) – służy głównie do przesyłania danych na serwer. Dane są wysyłane w treści (body) żądania HTTP.

Wybór metody zależy od serwera.

Spróbujmy teraz odpytać serwer jakie zbiory danych znajdują się na nim.

kolekcje <- stac_source |>
  collections() |>
  get_request()
kolekcje

## ###Collections
## - collections (9 item(s)):
##   - sentinel-2-pre-c1-l2a
##   - cop-dem-glo-30
##   - naip
##   - cop-dem-glo-90
##   - landsat-c2-l2
##   - sentinel-2-l2a
##   - sentinel-2-l1c
##   - sentinel-2-c1-l2a
##   - sentinel-1-grd
## - field(s): collections, links, context

W odpowiedzi otrzymaliśmy obiekt, który jest wielopoziomową listą list. Można samodzielnie sprawdzić jego strukturę używając funkcji View(). Eksploracja takich obiektów jest bardziej skomplikowana w porównaniu do obiektów jednowymiarowych (np. wektorów). Przykładowo, jeżeli chcemy sprawdzić pozostałe dostępne zbiory danych, to należy wykorzystać funkcje lapply(), która wykonuje iterację po każdym elemencie danych wejściowych (w tym przypadku listy), stosując przy tym podaną funkcję. Jako wynika zwracana jest lista, którą następnie można uprościć do wektora używając funkcji unlist().

Dodatkowo, dla uproszczenia możemy zastosować funkcję lambda (inaczej funkcję anonimową). Pozwala to zamienić słowo kluczowe function() na znak \(), np. zamiast lapply(lista, function(x) x^2) będzie lapply(lista, \(x) x^2).

kolekcje_nazwy = unlist(lapply(kolekcje$collections, \(x) x$id))
kolekcje_nazwy

## [1] "sentinel-2-pre-c1-l2a" "cop-dem-glo-30"        "naip"                 
## [4] "cop-dem-glo-90"        "landsat-c2-l2"         "sentinel-2-l2a"       
## [7] "sentinel-2-l1c"        "sentinel-2-c1-l2a"     "sentinel-1-grd"

Uwaga! Dane z Sentinela-2 znajdują się w dwóch kolekcjach, tj. sentinel-2-l2a (starsza) i sentinel-2-c1-l2a (nowsza). Aktualnie trwa aktualizacja produktów do nowszej kolekcji sentinel-2-c1-l2a i ta powinna być preferowana, natomiast starsza kolekcja w przyszłości nie będzie obsługiwana.

Teraz wyszukajmy dostępne produkty w ramach kolekcji sentinel-2-c1-l2a. Do wyszukiwania danych służy funkcja stac_search() i umożliwia ona zdefiniowanie takich parametrów jak:

kolekcja (collections)
zakres przestrzenny w układzie WGS 84 (bbox)
interwał czasowy w standardzie RFC 3339 (datetime)
maksymalna liczba produktów (limit)

W tym przypadku żądanie musimy wysłać za pomocą metody POST.

stac_source |>
  stac_search(
    collections = "sentinel-2-c1-l2a",
    bbox = c(22.5, 51.1, 22.6, 51.2), # xmin, ymin, xmax, ymax (WGS84)
    datetime = "2023-01-01T00:00:00Z/2023-12-31T00:00:00Z", # RFC 3339
    limit = 5) |>
  post_request() -> obrazy
obrazy

## ###Items
## - matched feature(s): 288
## - features (5 item(s) / 283 not fetched):
##   - S2A_T34UEB_20231230T094411_L2A
##   - S2A_T34UFB_20231230T094411_L2A
##   - S2A_T34UEB_20231227T093411_L2A
##   - S2A_T34UFB_20231227T093411_L2A
##   - S2B_T34UEB_20231225T094326_L2A
## - assets: 
## aot, blue, cloud, coastal, granule_metadata, green, nir, nir08, nir09, preview, product_metadata, red, rededge1, rededge2, rededge3, scl, snow, swir16, swir22, thumbnail, tileinfo_metadata, visual, wvp
## - item's fields: 
## assets, bbox, collection, geometry, id, links, properties, stac_extensions, stac_version, type

Zostało znalezionych 288 produktów, które spełniają zadane warunki. Możemy zobaczyć jakie warstwy (assets) i atrybuty (fields) są powiązane. Sprawdźmy przykładowe cechy produktu (properties).

# pierwsze 15 cech
names(obrazy$features[[1]]$properties)[1:15]

##  [1] "created"              "platform"             "constellation"       
##  [4] "instruments"          "eo:cloud_cover"       "proj:epsg"           
##  [7] "proj:centroid"        "mgrs:utm_zone"        "mgrs:latitude_band"  
## [10] "mgrs:grid_square"     "grid:code"            "view:azimuth"        
## [13] "view:incidence_angle" "view:sun_azimuth"     "view:sun_elevation"

Mamy dostęp do 44 cech, które opisują produkt, m. in. data pozyskania, nazwa platformy satelitarnej, procentowe zachmurzenie sceny, układ współrzędnych, ID sceny, itd.

Sprawdźmy zachmurzenie wyszukanych scen. Określone jest one w atrybucie eo:cloud_cover. W tym celu ponownie musimy wykorzystać funkcję lapply().

unlist(lapply(obrazy$features, \(x) x$properties$"eo:cloud_cover"))

## [1] 71.93396 68.93441 83.81743 94.01084 99.99993

Wyszukane sceny charakteryzują się wysokim zachmurzeniem w związku z czym ich przydatność jest ograniczona. Oczywiście możemy dokonać prostej filtracji modyfikując nasze poprzednie zapytanie. Do filtrowania przeznaczona jest funkcja ext_query().

stac_source |>
  stac_search(
    collections = "sentinel-2-c1-l2a",
    bbox = c(22.5, 51.1, 22.6, 51.2),
    datetime = "2023-01-01T00:00:00Z/2023-12-31T00:00:00Z",
    limit = 5) |>
  ext_query(`eo:cloud_cover` < 10) |>
  post_request() -> obrazy

unlist(lapply(obrazy$features, \(x) x$properties$"eo:cloud_cover"))

## [1] 0.067071 0.043523 1.478151 0.678711 0.000757

Alternatywnie, zamiast używać wielopoziomowych list, możemy dokonać konwersji obiektu obrazy do “przestrzennej ramki danych” (simple features data frame) za pomocą funkcji items_as_sf(). Jednakże, na ten moment jest ona eksperymentalna i mogą pojawić się niespodziewane problemy.

df = items_as_sf(obrazy)
plot(sf::st_geometry(df)[1], main = "Zasięg sceny", axes = TRUE)

Po dokonaniu selekcji scen, możemy przystąpić do ich pobrania. Zasadniczo, można to wykonać na kilka różnych sposobów używając:

funkcji assets_download() z pakietu rstac
funkcji download.file() wbudowanej w R
modułu /vsicurl/ dostępnego w pakietach obsługujących GDAL, np. terra
narzędzi GDAL (np. gdalwarp) dostępnych m. in. w powłoce OSGeo4W na Windows

assets_download()

W funkcji assets_download() możemy zdefiniować, które warstwy chcemy pobrać (argument asset_names). Jednak ma ona pewien mankament związany z dość skomplikowanym zapisem ścieżek do pobranych plików wynikający ze struktury katalogów STAC.

obrazy |>
  items_select(1) |>
  assets_download(asset_names = c("blue", "green", "red"))

Sceny satelitarne zazwyczaj zawierają niskorozdzielcze miniaturki w kompozycji RGB. W przypadku danych z Sentinela-2 zapisane są one jako thumbnail.jpg. Jeśli chcemy szybko podejrzeć dane, to możemy odwołać się do tej warstwy i pozyskać do niej odnośnik przy pomocy funkcji assets_url(). Do wczytania pliku .jpg wymagany jest dodatkowy pakiet, np. imager, który posiada funkcję load.image().

library("imager")
par(mar = c(0, 0, 0, 0)) # usuń marginesy

obrazy |>
  items_select(1) |>
  assets_select(asset_names = "thumbnail") |>
  assets_url() |>
  load.image() |>
  plot(axes = FALSE)

download.file()

Tak jak w powyższym przypadku, możliwe jest pozyskanie odnośników do pozostałych warstw i pobranie danych przy użyciu funkcji download.file(), wymagającej wskazania odnośników, z których zostaną pobrane dane oraz ścieżek zapisu.

obrazy |>
   items_select(1) |>
   assets_select(asset_names = c("blue", "green", "red")) |>
   assets_url() -> urls
urls

## [1] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/11/S2A_T34UEB_20231107T093339_L2A/B02.tif"
## [2] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/11/S2A_T34UEB_20231107T093339_L2A/B03.tif"
## [3] "https://e84-earth-search-sentinel-data.s3.us-west-2.amazonaws.com/sentinel-2-c1-l2a/34/U/EB/2023/11/S2A_T34UEB_20231107T093339_L2A/B04.tif"

Utwórzmy nowy folder o nazwie sentinel używając funkcji dir.create(). Do stworzenia ścieżek zapisu należy podać nazwę pliku wraz z rozszerzeniem. W tym celu można zastosować funkcję basename() do wyodrębnienia tych nazw z odnośników. Następnie należy połączyć ścieżkę do katalogu ze ścieżką do pliku wykorzystując funkcję file.path().

dir.create("sentinel")
rastry = file.path("sentinel", basename(sentinel))

Finalnie, zastosujemy pętlę, która pobierze wszystkie pliki do lokalizacji określonej przez nas.

for (i in seq_along(sentinel)) {
  download.file(sentinel[i], rastry[i], mode = "wb")
}

/vsicurl/

Ostatnim sposobem jest wykorzystanie modułu /vsicurl/ pochodzącego z GDAL bez potrzeby pobierania całego zbioru danych. Działanie jest proste – wystarczy dodać przedrostek /vsicurl/ do odnośników z rastrami używając funkcji paste0().

urls = paste0("/vsicurl/", urls)
r = rast(urls)
r

## class       : SpatRaster 
## dimensions  : 10980, 10980, 3  (nrow, ncol, nlyr)
## resolution  : 10, 10  (x, y)
## extent      : 499980, 609780, 5590200, 5700000  (xmin, xmax, ymin, ymax)
## coord. ref. : WGS 84 / UTM zone 34N (EPSG:32634) 
## sources     : B02.tif  
##               B03.tif  
##               B04.tif  
## names       :   B02,   B03,   B04 
## min values  :     1,     1,   473 
## max values  : 18640, 17728, 17104

Więcej informacji na temat wyszukiwania i pozyskiwania danych za pomocą standardu STAC znajdziesz w oficjalnym tutorialu Download data from a STAC API using R, rstac, and GDAL napisanym przez Michael Mahoney.

Zadanie

9. Pobierz dowolną scenę satelitarną z Sentinela-2 o niskim zachmurzeniu (maksymalnie 20%) i następnie:

sprawdź metadane rastrów
przygotuj wizualizację RGB
pobierz losową próbę 10 tys. punktów dla kanału niebieskiego, zielonego, czerwonego oraz bliskiej podczerwieni i zaprezentuj statystyki opisowe oraz porównaj histogramy
dodatkowo dla kanału czerwonego oraz bliskiej podczerwieni wykonaj wykres rozrzutu oraz oblicz współczynnik korelacji Pearsona
oblicz znormalizowany różnicowy wskaźnik wegetacji (NDVI) i przygotuj wizualizację

Dostęp do danych satelitarnych

Krzysztof Dyba

Wprowadzenie

Dane satelitarne

Przetwarzanie potokowe

Pozyskiwanie danych

Zadanie