Wprowadzenie
Ogrom informacji meteorologicznych wytwarzanych każdego dnia na świecie wymaga odpowiedniej archiwizacji, możliwości przesyłania i udostępniania ich dla użytkowników. W tym celu muszą być one zapisywane w odpowiednich formatach, umożliwiających uniwersalny dostęp i duży stopień kompresji danych.
Niestety wciąż nie istnieje jeden idealny format danych meteorologicznych, który byłby wykorzystywany przez wszystkie ośrodki na całym świecie. Obecnie możemy natrafić na dziesiątki różnych formatów plików, w zależności od typu danych, oprogramowania używanego do ich wytworzenia czy instytucji odpowiedzialnej za te dane.
Niektóre formaty danych meteorologicznych są proste do odczytania i edycji bez potrzeby instalacji dodatkowego oprogramowania, inne wymagają odpowiednich bibliotek, niekiedy ograniczonych do określonego systemu operacyjnego i architektury komputera.
Dzisiaj poznamy te najbardziej popularne formaty danych meteorologicznych, spróbujemy pobrać je i odczytać za pomocą R.
Wybrane formaty danych meteorologicznych
Na przestrzeni kilkudziesięciu lat wraz z rozwojem cyfryzacji, powstało wiele formatów, w których przechowywane są dane meterologiczne. Niekiedy są to formaty charakterystyczne dla danej instytucji, niekiedy standardy formatu danych określone np. przez WMO (World Meteorological Organization). Najważniejsze z nich to:
- GRIB (GRidded Binary Edition 1 i 2) – opracowane przez WMO
- netCDF (Network Common Data Form Version 3.x i 4.x) – opracowane przez Unidata (UCAR/NCAR)
- HDF i HDF-EOS (Hierarchical Data Format-Earth Observing System) – opracowane przez NCSA/NASA
- GeoTIFF – raster z georeferencją
- BUFR (Binary Universal Form) – opracowane przez WMO
- Inne formaty binarne
- ASCII, ANSI, XML i inne
Standardy WMO
W latach 80. XX wieku, w celu ujednolicenia sposobu kodowania danych meteorologicznych i konieczności kompresji, WMO opracowała dwa binarne formaty danych: GRIB i BUFR.
Na pewno był to postęp w stosunku do wcześniej stosowanych praktyk, jednakże oba te formaty nie uniknęły wad i niedoskonałości. Największą ich słabością jest zależność od tzw. tablic kodujących, które jednoznacznie określają, jakie dane znajdują się w plikach. Jeżli podmiot wytwarzający dane meteorologiczne nie udostępnia odpowiednich tablic, to możemy się tylko domyślać lub zgadywać co na pewno znajduje się w pliku.
Drugą kwestią jest niejednoznaczność i niedopowiedzenia co do tego, jak dokładnie powinny wyglądać pliki w standardach WMO. Same wytyczne niekiedy nie dostarczają odpowiedzi, a struktura pliku dopuszcza pewną dowolność w zapisie danych. Nigdy nie powstało również oprogramowanie jednoznacznie umożliwiające rozkodowanie informacji i ujednolicające wiele tabeli rozkodowujących, jakie można znaleźć w internecie.
Tutaj wybija się pewna przewaga formatu BUFR - umożliwia on zapis takiej tabeli w tym samym pliku w sekcji poprzedzającej dane meteorologiczne.
GRIB
Powstał w roku 1985 i jest standardem danych dla wyników modelowania meteorologicznego w siatce. Od tamtej pory ukazały się trzy wersje tego formatu: GRIB0, GRIB1 i GRIB2. O ile GRIB0 jest bardzo rzadko spotykany, to zarówno GRIB1 jak i GRIB2 cieszą się względnie dużą popularnością.
W formacie GRIB1 każdy rekord zawiera informację o wymiarach w poziomie, jednym poziomie w pionie atmosfery i kroku czasowym.
GRIB2 umożliwia dalszą kompresję danych poprzez umożliwienie zapisywania w jednym rekordzie informacji w wielu poziomach dla każdego kroku czasowego.
Ze względu na dosyć dużą popularność formatu GRIB, powstało kilka narzędzi umożliwiających ich konwersję lub wizualizację. Również wiele bibliotek i pakietów do obsługi danych w siatce umożliwia odczytanie danych w tym formacie. Wybrane narzędzia warte uwagi: GrADS, wgrib/wgrib2, CDO, GDAL, NCL, RNomads, degrib, Pygrib.
BUFR
Utworzony w 1988 r. w celu zapisu obserwacji meteorologicznych, również operacyjnie do prognozy pogody. Używany głównie do danych punktowych, początkowo miał zastąpić depesze meteorologiczne (SYNOP, CLIMAT i inne).
Podobnie jak w przypadku formatu GRIB, format BUFR jest dostosowywany do potrzeb instytucji, które wytwarzają te pliki. W rezultacie pomimo bycia tworzonym zgodnie ze standardem WMO, w takim pliku mogą znajdować się dodatkowe dane lub nieco odmienny sposób kodowania informacji.
Dzięki temu, że tabele rozkodowujące mogą znajdować się wewnątrz pliku BUFR, to czyni ten format samoopisującym (ang. self-describing).
Niektóre narzędzia umożliwiające dekodowanie plików BUFR: ECMWF BUFRDC, NCEP BUFRLIB, BUFRdisplay, GrADS, PyBUFR.
netCDF
To format danych atmosferycznych i oceanicznych utworzony w 1989 roku przez Unidata Program Center (UCAR) w USA. Obecnie jest bardzo szeroko stosowanym, samoopisującym się formatem.
Olbrzymią zaletą tego formatu plików jest fakt, że dowolny fragment może być czytany osobno, bez wcześniejszego wczytywania całego pliku i rozkodowywania informacji.
Jednocześnie, do pliku można dopisywać dane w jednym wymiarze bez edycji istniejącej części danych.
Na ten moment istnieją dwie wersje plików netCDF. Pierwsza z nich, netCDF-3 lub netCDF-classic, była używana przez wiele lat, jednak wraz ze wzrostem stopnia skomplikowania danych zaczęto dostrzegać jej ograniczenia.
Wersja netCDF-4 powstała w odpowiedzi na rosnące oczekiwania użytkowników, a ze względu na to, że jest to format będący niejako hybrydą HDF5 i netCDF-3, to charakteryzuje się wysoką popularnością i uniwersalnością.
Format netCDF jest obsługiwany przez wiele komercyjnych i otwartych programów GIS-owych, a także szereg bibiliotek służących do analiz danych przestrzennych.
HDF
Utworzony przez NCSA (National Center for Supercomputing Applications) w USA na początku lat 90. XX wieku jako format przeznaczony do przechowywania danych naukowych.
Obecnie format HDF jest utrzymywany przez HDF Group, który wspiera rozrastanie się sieci danych modelowych, dostarcza API do przetwarzania danych i pozwala na dodawanie nowych treści.
HDF jest formatem samoopisującym, pozwalającym na interpretację zawartości plików bez dodatkowej informacji z zewnątrz. Składa się z mieszanki powiązanych obiektów, do których można się odwoływać pojedynczo albo kolektywnie. Może być wykorzystywany do przechowywania danych w formie macierzy, tabeli, adnotacji tekstowych, niektórych typów rastrów i metadanych.
Obecnie istnieją dwie wersje formatu HDF: HDF4 i HDF5. HDF4 jest wcześniejszą wersją posiadającą pewne ograniczenia, m. in. limit rozmiaru pliku do 2 GB, co przy dzisiejszych zastosowaniach jest bardzo limitujące.
HDF5 poprawia niektóre problemy z wcześniejszą wersją, zawiera dane które są lepiej uporządkowane i jednorodne, co ułatwia pracę z tymi plikami. Z powodzenia używany jest nie tylko do danych meteorologicznych, ale też do danych giełdowych czy monitoringowych.
Rozróżnia się także format HDF-EOS, który związany jest bardzo ściśle z danymi o geoprzestrzennym znaczeniu (siatki, punkty itp.). Dostęp do tych danych i wybieranie konkretnych treści odbywa się poprzez zadanie odpowiednich współrzędnych i czasu (o ile współrzędne geograficzne zostały nadane). Narzędzia, które pozwalają na odczytanie danych HDF, będą również pracować z HDF-EOS, jednakże standardowe biblioteki HDF będą gorzej radzić sobie z metadanymi charakterystycznymi dla HDF-EOS niż biblioteki HDF-EOS.
Do tej pory powstało wiele narzędzi obsługujących pliki HDF na wielu platformach (Java, MATLAB, Scilab, Python, R), a także program do wizualizacji oparty na Javie (HDFView).
GeoTIFF
Najprostszym sposobem do zapisu danych meteorologicznych jest raster z georeferencją - GeoTIFF. Jest to w zasadzie macierz, która w prosty sposób może być wizualizowania i analizowana w każdym programie i bibliotece GIS. Ten format jest dobry do przechowywania danych klimatologicznych (wieloletnie trendy temperatury, średnie miesięczne wartości itp.). Jednakże, kiedy zależy nam na czwartym wymiarze - czasie, to każda “klatka” musi być wtedy przechowywana w osobnym pliku bądź kanale (band).
Źródła modelowych danych meteorologicznych
Najważniejsze źródła globalnych danych meteorologicznych pochodzących z modeli:
- RDA UCAR - zasoby gromadzące dane z wielu źródeł światowych i krajowych, również z pomiarów meteorologicznych w wielu formatach
- NOAA - dane modelowe z wielu modeli NOAA
- NASA GMAO - produkty modelowania w NASA Global Modeling and Assimilation OffiCe
- ECMWF - Europejskie Centrum Prognoz Średnioterminowych
- CHELSA - dane klimatologiczne opracowane w oparciu o modele globalne
Wybrane pakiety do analizy danych meteorologicznych w R
Pakiet RNCEP
Pakiet RNCEP pozwala na pobieranie, zarządzanie i wizualizację danych meteorologicznych pochodzących z dwóch zestawów reanaliz:
Aby uzyskać dostęp do tych danych, należy uzyskać najpierw odpowiedni token ze strony NOAA. Na podany adres e-mail zostanie wysłany kod, który upoważni nas do pobierania danych. Należy wykorzystać go do ustawienia zmiennej środowiskowej R, najłatwiej wpisując w R Studio poniższy kod:
Sys.setenv(NOAA_KEY = "Twój kod")NCEP.gather
To funkcja pobierająca dane reanaliz z serwera NOAA. Umożliwia pobieranie jednego parametru w zadanym okresie czasu dla określonego obszaru. Aby wiedzieć, jakie są nazwy parametrów, należy sięgnąć do dokumentacji pakietu RNCEP. Poniżej przykładowe użycie funkcji do pobrania parametru air.sig995 (temperatura powietrza na poziomie sigma = 0.995, bliskim powierzchni Ziemi).
# Wczytanie pakietu
library(RNCEP)
# Pobranie danych z serwera NOAA i zapisanie do obiektu temp_2018
temp_2018=NCEP.gather(variable='air.sig995', level='surface',
months.minmax=c(7,8), years.minmax=c(2018),
lat.southnorth=c(46,56), lon.westeast=c(11,27),
reanalysis2 = FALSE, return.units = TRUE)## [1] Units of variable 'air.sig995' are degKNCEP.vis.area
To funkcja do wizualizacji pobranych danych. Wykorzystujemy wcześniej zapisany obiekt do sporządzenia mapki rozkładu temperatury. Przykład wykorzystania dla zadanej daty i godziny:
# Wizualizacja wybranego rozkładu temperatury według daty i godziny
NCEP.vis.area(temp_2018, layer='2018-07-15 12', show.pts=FALSE, draw.contours=TRUE,
cols=heat.colors(64), transparency=.5, axis.args=NULL, map.args=NULL,
grid.args=NULL, title.args=list(main="Temperatura na 2m n.p.t - 15.07.2018 r. 12:00",
xlab = "Długość geograficzna", ylab = "Szerokość geograficzna"),
interp.loess.args=NULL,image.plot.args=list(legend.args=list(text="[K]", cex=1.15, padj=-1, adj=0)),
contour.args=list(labcex=.8), points.args=NULL)
NCEP.aggregate
To funkcja pozwalająca na obliczenie funkcji danego parametru według lat, miesięcy, godzin lub dni. Mogą to być podstawowe statystyki, jak np.: mean, max, min lub np. funkcja. Przykładowo, deklaracja funkcji zwracającej procent obserwacji poniżej danej wartości wygląda następująco: function(x){return(length(which(x < 273.15))/(length(x) - sum(is.na(x)))).
Zagregowane wartości danego parametru możemy również zwizualizować za pomocą wcześniej wspomnianej funkcji NCEP.vis.area.
# Obliczenie średniej temperatury miesięcznej
mean_t2 = NCEP.aggregate(wx.data=temp_2018, YEARS=FALSE, MONTHS=TRUE,
HOURS=FALSE, DAYS=FALSE, fxn='mean')
# Wizualizacja średniej temperatury w lipcu 2018
NCEP.vis.area(wx=mean_t2, layer=1, title.args=list(main='Średnia temperatura na 2m n.p.t. w lipcu 2018',
xlab='Długość geograficzna', ylab='Szerokość geograficzna'),
image.plot.args=list(legend.args=list(text="[K]", cex=1.15, padj=-1, adj=0)),
contour.args=list(labcex=.8), points.args=NULL, show.pts=FALSE)
GRIB w R
Istnieje kilka sposobów odczytywania danych GRIB w R. Jednym z nich jest potraktowanie pliku GRIB jak rastra i wykorzystanie pakietu raster. Inną możliwością jest pakiet stars, który pomaga odczytywać wiele rodzajów plików, w tym netCDF.
Spróbujmy zwizualizować plik GRIB pochodzący z miesięcznych średnich wartości temperatur powietrza reanalizy ECMWF ERA5 (ds630.1 na UCAR RDA).
Pakiet raster
# Wczytanie bibliotek
library(raster)
# Zwrócenie informacji o metadanych pliku
raster("Dane/e5.mnth.mean.an.sfc.128_167_2t.regn320sc.2008010100_2008120100.grb")## class : RasterLayer
## band : 1 (of 12 bands)
## dimensions : 640, 1280, 819200 (nrow, ncol, ncell)
## resolution : 0.2812502, 0.2810172 (x, y)
## extent : -0.1406251, 359.8596, -89.92551, 89.92551 (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +a=6367470 +b=6367470 +no_defs
## data source : C:\Users\aqpro\Desktop\PNoZII\Projekt_3\Dane\e5.mnth.mean.an.sfc.128_167_2t.regn320sc.2008010100_2008120100.grb
## names : e5.mnth.mean.an.sfc.128_167_2t.regn320sc.2008010100_2008120100# Wybranie piątego kanału z rastra
r <- raster("Dane/e5.mnth.mean.an.sfc.128_167_2t.regn320sc.2008010100_2008120100.grb", band=5)
# I wizualizacja
plot(r)
Pakiet stars
# Wczytanie bibliotek
library(stars)
library(cptcity)
# Wczytanie pliku GRIB z dysku twardego
met <- read_stars("Dane/e5.mnth.mean.an.sfc.128_167_2t.regn320sc.2008010100_2008120100.grb")
# Wizualizacja w zadanej palecie kolorów
plot(met, col = cpt())
# Policzenie średniej wartości "piksela" w naszym pliku
met_mean <- st_apply(met, 1:2, mean)
# Wizualizacja
plot(met_mean, col = cpt())
netCDF w R
Pliki netCDF również możemy wczytywać w R na kilka sposobów. Powstało kilka pakietów do tego dedykowanych. Dla przykładu, wykonajmy kilka funkcji w pakiecie ncdf4. To podstawowa biblioteka pozwalająca na edycję plików netCDF - łącznie z dodawaniem nowych zmiennych.
Plik, jaki mamy do dyspozycji, to output z modelu meteorologicznego WRF. Obrazuje on przejście burzy śnieżnej nad wschodnim wybrzeżem Stanów Zjednoczonych w styczniu 2000 roku. Zwizualizujemy zmienną QVAPOR - stosunek zmieszania pary wodnej w pewnej odległości nad powierzchnią Ziemi.
### Wczytywanie pakietu
library(ncdf4)
#Otwieranie pliku netcdf
wrf_2000 <- nc_open('Dane/wrfout_d01_2000-01-24_12')
# Aby uzyskać informacje o wszystkim metadanych w pliku uruchom:
#print(wrf_2000)
# Wybierzmy zmienną opisującą wartość stosunku zmieszania pary wodnej
QVAPOR <- ncvar_get(wrf_2000, varid="QVAPOR")
# Wczytanie pakietów
require(reshape)
require(ggplot2)## Error in get(genname, envir = envir) : nie znaleziono obiektu 'group_map'
## Error in get(genname, envir = envir) :
## nie znaleziono obiektu 'group_split'# Zamiana czterowymiarowego obiektu na ramkę danych z 5 kolumnami (X1, X2, X3, X4 i value)
QVAPOR_df = melt(QVAPOR)
# Wizualizacja piątego kroku czasowego na dziesiątym poziomie zmiennej QVAPOR
ggplot(aes(x = X1, y = X2, fill = value), data = QVAPOR_df[QVAPOR_df$X3==10 & QVAPOR_df$X4==5,]) +
geom_raster()
# Zamknięcie pliku
nc_close(wrf_2000)