0. Wprowadzenie

Pliki NetCDF początkowo mogą sprawiać bardzo dużo problemów w obsłudze ze względu na bardzo skomplikowaną strukturę danych jaka może być zapisana w tym formacie. Jest to obecnie jednak najbardziej podstawowy format udostępniania danych meteorologicznych, więc warto się stopniowo do niego przyzwyczajać… Szczegółowa dokumentacja do formatu NetCDF dostępna jest na stronie: http://www.unidata.ucar.edu/software/netcdf/docs/index.html

1. Instalacja rozszerzeń

Teoretycznie do pracy z plikami zapisanymi w formatach NetCDF-3 oraz NetCDF-4 wystarczy biblioteka ‘ncdf4’. W praktyce (zwłaszcza na początku pracy z plikami .nc) warto na bieżąco wizualizować wyniki pracy w celu sprawdzenia poprawności poszczególnych procedur.

Rozpocznijmy do zainstalowania pakietów:

ncdf4 - pakiet dający dostęp do bibliotek NetCDF
magic - rozszerzenie pozwalające na łatwe obracanie macierzy
fields - pakiet rozszerzający funkcję image dodając skalę (funkcja: image.plot)
maps - pakiet zawierący mapę konturową świata

Następnie aktywizujmy zainstalowane pakiety:

library(fields)
library(maps)
library(ncdf4)
library(magic)

2. Pobierz plik NetCDF

Ze strony http://openmeteo.pl/html/tg_0.25deg_reg_v11.0.nc pobierz plik zawierający średnie dobowe temperatury powietrza z europejskiej wersji reanalizy ECA&D od roku 1950. W celu zmniejszenia objętości pliku dane zostały przycięte do współrzędnych zawierających obszar Polski. Rozdzielczosć tej bazy danych to 0.25x0.25 i zawiera ona informacje jedynie dla punktów lądowych. Zapisz plik w znanej lokalizacji, którą użyjesz w kolejnym punkcie.

3: Utwórz połączenie do pliku

3.1: Wskaż katalog roboczy

Katalog roboczy powinien wskazywać na lokalizację z pobranym plikiem ## 3.2. Utwórz połączenie do pliku nc za pomocą funkcji nc_open. Zapoznaj się z dokumentacją tej funkcji: (?nc_open)

nc <- nc_open("tg_0.25deg_reg_v11.0.nc")

Stworzony obiekt (ze względu na stopień skomplikowania) jest listą (obiektem klasy ‘list’). Na obiektach tej klasy pracowaliśmy dotychczas w bardzo niewielkim stopniu. Warto zatem zapoznać się ze specyfiką tego typu obiektów, np. korzystając z podręcznika P. Biecka (http://www.biecek.pl/R/PrzewodnikPoPakiecieRWydanieIIIinternet.pdf), strona 75.

Teoretycznie obiekt (jak każdy inny) można podejrzeć za pomocą funkcji str. Bardziej adekwatne w tym przypadku będzie użycie nazw list zagnieżdżonych wewnątrz obiektu ‘nc’

names(nc)
##  [1] "filename"    "writable"    "id"          "safemode"    "format"     
##  [6] "is_GMT"      "groups"      "fqgn2Rindex" "ndims"       "natts"      
## [11] "dim"         "unlimdimid"  "nvars"       "var"

W praktyce, spośród wszystkich informacji, które są zapisane w pliku niezbędnych jest tylko kilka spośród setek informacji i kilka metadanych:

  • ndims - liczba wymiarów przechowywanych w pliku
  • natts - liczba atrybutów
  • dim - metadane i wartości dotyczące wymiarów (najczęściej: długość, szerokość geograficzna i czas)
    • Item 2a
  • nvars - liczba zmiennych (parametrów) przechowywanych w pliku
  • var - metadane i wartości dla wszystkich zmiennych w pliku

4. Rozpoznanie struktury dla wymiaru czasu:

Nazwy wymiarów możemy sprawdzić wpisując komendę, która wydobędzie nam nazwy z podlisty dim

names(nc$dim)
## [1] "latitude"  "longitude" "time"

Większość informacji dla zmiennej czasu można pobrać z podlisty: nc$dim$time. Zawartość podlist zawsze warto przejrzeć w pierwszej kolejności za pomocą funkcji ‘str’
Spójrzmy na pierwszych kilka wartości (head, n=ilość rzędów), które są zawarte w obrębie zmiennej vars :

head(nc$dim$time$vals, 50)
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22
## [24] 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
## [47] 46 47 48 49

Są to kolejne liczby naturalne… Co zatem mogą one oznaczać? Tą informację można sprawdzić w zmiennej units przechowującej jednostki dla danego wymiaru

nc$dim$time$units
## [1] "days since 1950-01-01 00:00"

Kolejne wartości wektora oznaczają liczbę dni od 1. stycznia 1950 roku.
Jako, że czas jest często wykorzystywaną zmienną warto zapisać go do nowego obiektu o krótkiej nazwie, np. czas.
W pierwszej kolejności musimy utworzyć obiekt, który będzie naszym punktem początkowym, czyli datę 01-01-1950, która będzie rozpoznawalna dla R. Najprostsza klasa takiego obiektu to “Date”, zatem:

daty <- as.Date("1950-01-01")

W kolejnym kroku dodamy do naszego obiektu kolejne dni, które będą odpowiadały zawartości obiektu nc$dim$time$vals i nadpiszemy obiekt daty

daty <- as.Date("1950-01-01")+nc$dim$time$vals
# i sprawdźmy profilaktycznie pierwszych kilka dat:
head(daty)
## [1] "1950-01-01" "1950-01-02" "1950-01-03" "1950-01-04" "1950-01-05"
## [6] "1950-01-06"
# oraz zakres wartości:
range(daty)
## [1] "1950-01-01" "2014-12-31"

5: Zadanie (1):

Utwórz zmienną lon i lat przechowującą odpowiednio wartości dla długości i szerokości geograficznej. ## 5.1. Podaj ile jest gridów dla długości i szerokości geograficznej ## 5.2. Jakie podano jednostki dla obu wymiarów? ## 5.3. Podaj jaki jest zakres wartości dla obu utworzonych obiektów

…

6: Pobieranie wartości zmiennych

Pobieranie wartości zmiennych (wartości poszczególnych parametrów) wygląda nieco inaczej niż pobranie wymiarów # 6.1. Nazwy zmiennych W najszybszy sposób:

names(nc$var)
## [1] "tg"

W naszym pliku jest tylko jedna zmienna nazwana ‘tg’. Sprawdź zatem zawartość listy przechowującej wartości dla analizowanej zmiennej:

nc$var$tg

W praktyce zdecydowanie szybciej można pobrać wartości zmiennych (i wymiarów) za pomocą komendy ncvar_get. Sprawdź jednocześnie wymiar utworzonej (wielowymiarowej) macierzy:

tg <- ncvar_get(nc,varid = "tg")
dim(tg)
## [1]    44    26 23741

Otrzymane powyżej wartości oznaczają kolejno wymiary dla współrzędnej X (lon), Y (lat) oraz T (czas)

7: Wizualizacja (diagnostyka)

W celu sprawdzenia czy nasze dane wczytały się poprawnie wykorzystaj funkcję image do wizualizacji danych 2-wymiarowych. Aby wyświetlić pierwszą warstwę czasową wystarczy, że pobierzemy wszystkie wymiary dla osi X, Y, oraz zdefiniujemy dla osi T pierwszą warstwą czasową, czyli

image(tg[,,1])

Domyślnie nie jest konieczne definiowane w funkcji image argumentów dla osi X i Y (tworzone są domyślnie wartości od 0 do 1). Jeśli chcemy poprawić nasz wykres o odpowiednie zakresy na osiach wpiszmy jako 2 pierwsze argumenty wektory wartości lon i lat

image(lon,lat,tg)