Podstawy nauk o Ziemi II

Wprowadzenie

Wielokrotnie w trakcie zajęć lub też w przyszłej pracy będą Państwo pracować z nowymi typami danych pochodzącymi z różnych źródeł. Zapewne w większości przypadków będą one reprezentowane za pomocą najbardziej znanych i uniwersalnych formatów zapisu danych, opatrzone dokładnym opisem wszystkich zmiennych. Niekiedy jednak zdarzy się, że pliki w bazie danych będą przechowywane w formacie pliku, który ma na celu maksymalną kompresję informacji – tylko za pomocą znaków alfanumerycznych, bez znaków podziału kolumn, bez nagłówków – surowe dane. I z takimi danymi będziemy mieć dzisiaj do czynienia.

The first rule of data processing is to look at your data! - Dennis Shea (NCAR)

Praca z tego typu plikami, które są niejako “zaszyfrowane” wymaga dodatkowych operacji – przestudiowania dołączonej dokumentacji (oby była!), w której znajdziemy dokładny opis struktury plików, a także poszukiwania potencjalnych narzędzi, które automatycznie za nas przetworzą je do formatu łatwiejszego do interpretacji i dalszej analizy. Kiedy takich nie znajdziemy, a spodziewamy się w przyszłości wracać do tych plików, to oczywiście można samemu utworzyć skrypt przyspieszający ich przetwarzanie.

W międzyczasie mamy na uwadze to, do czego wykorzystamy dane i jakie analizy nas interesują – aby wziąć to pod uwagę podczas porządkowania plików.

Z powyższych względów poniższy skrypt został podzielony na kilka części:

Zapoznanie się z formatem danych w bazie danych meteorologicznych ISD NOAA,
Poznanie narzędzi przyspieszających ich przetworzenie,
Wstępna analiza statystyczna danych w celu odnalezienia trendów i zależności między parametrami,
Opis projektu do samodzielnego wykonania.

Większość pracy przeprowadzona zostanie z wykorzystaniem środowiska R.

1. Opis bazy ISD NOAA

Baza ISD NOAA (Integrated Surface Database of the National Oceanic and Atmosferic Administration) jest obecnie jednym z największych zbiorów bezpłatnie udostępnianych meteorologicznych danych pomiarowych. Zawiera ona obserwacje uśredniane w czasie jednej godziny dla 12922 naziemnych punktów pomiarowych z 237 krajów świata.

W roku 1998 amerykańskie ośrodki naukowe (NCDC, FCC) przy współpracy z siłami powietrznymi i marynarką wojenną Stanów Zjednoczonych dostrzegły potrzebę utworzenia jednolitej bazy pomiarów klimatologicznych przede wszystkim do zastosowań badawczych¹. Korzystając z wielu zewnętrznych źródeł zaczęto gromadzić ogromne ilości obserwacji meteorologicznych w różnych formatach danych i integrować je do jednej globalnej bazy w zunifikowanym formacie plików. Pierwsze dane pochodzą z 1901 roku², a najstarsze pomiary z terenu Polski z roku 1928. Bieżące dane pomiarowe są uzupełniane co kilka dni po przeprowadzeniu testów jakości i transformacji do formatu ISH.

Pomiary meteorologiczne w formacie ASCII udostępniane są bezpłatnie przez serwer FTP, jednak pobieranie ich w sposób bezpośredni jest żmudne i pracochłonne. Format w jakim są przechowywane dane wymaga dodatkowego przetworzenia, aby możliwa była ich swobodna analiza.

W kolejnych podrozdziałach zamieszczono strukturę plików oraz najważniejsze informacje na temat zawartych w nich parametrach, a także sposób ich przechowywania w plikach.

I choć przetworzenie danych w klasyczny sposób jest możliwe, to warto skorzystać z gotowych rozwiązań, aby móc przejść od razu do ich analizy.

Struktura plików

Dane pomiarowe z bazy ISD NOAA przechowywane są za pomocą znaków ASCII w plikach formatu ISH (Integrated Surface Hourly - Compressed Archive). Informacje dotyczące tego, jakie zmienne i w jakich jednostkach są udostępniane znajdują się w szczegółowej dokumentacji formatu. Wartość każdej z dostępnych zmiennych zajmuje określone pozycje w rekordzie pliku, tak więc bez znajomości jego dokładnej struktury każda linijka jest tylko zakodowanym rzędem cyfr i liter. Poza 34 zmiennymi, które są obligatoryjne, czasem znajdziemy również dodatkowe informacje w pliku. Wśród najważniejszych zmiennych możemy wyróżnić wartości wybranych pomiarów meteorologicznych, informacje dotyczące stacji i źródła obserwacji, a także wyniki testów jakości dla poszczególnych serii pomiarowych.

W poniższej tabeli zebrano najważniejsze parametry dotyczące stacji pomiarowej oraz typy pomiarów meteorologicznych przechowywane w plikach z bazy ISD NOAA, a także techniczne informacje dotyczące struktury plików.

Parametr	Jednostka	Mnożnik	Brak danych	Pozycja w rekordzie
Numer USAF	-	-	-	5-10
Numer WBAN	-	-	-	11-15
Data pomiaru (YYYYMMDD)	-	-	-	16-23
Czas pomiaru (GMT-UTC)	-	-	-	24-27
Długość/szerokość geogr.	[°]	1000	99999	29-34 i 35-41
Wysokość n.p.m.	[m]	1	9999	47-51
Kierunek wiatru	[°]	1	999	61-63
Prędkość wiatru	[m/s]	10	9999	66-69
Temperatura powietrza	[°C]	10	9999	88-92
Temperatura punktu rosy	[°C]	10	9999	94-98
Ciśnienie (na poziomie morza)	[hPa]	10	99999	100-104

Pobieranie i przetwarzanie plików

Pobieranie plików z bazy

Jeżeli chcielibyśmy pozyskać oraz przetworzyć manualnie pliki z bazy ISD NOAA, to możemy to zrobić na kilka sposobów, które oferuje serwis internetowy NOAA:

Baza FTP – domyślnie i najprościej, jednak wymaga ręcznego wyszukania numeru stacji, która nas interesuje (zawartego w pliku isd-history.csv). Jest to sposób dobry, jeżeli faktycznie potrzebujemy tylko danych dla jednego punktu z wybranego roku lub używamy oprogramowania pośredniego do pobierania danych (np. w R funkcja download.file).
Wyszukiwarka oraz mapka – pozwalają na pobieranie niewielkiego zbioru wybranych parametrów w oparciu o zadane kryteria.
API (JSON) – za pomocą prostego kodu możemy pobierać ograniczony zestaw danych. Przydatne również do innych danych udostępnianych przez NOAA.

Rozpakowanie plików z bazy FTP

Zaczynamy od rozpakowania skompresowanego pliku GZIP. W systemie Windows możemy użyć odpowiedniego programu do dekompresji bądź konsoli. Korzystając z R możemy to zrobić tak:
system("gunzip -r /ścieżka i nazwa pliku/"), intern = FALSE, ignore.stderr = TRUE)

Edycja formatu ISH

Do automatycznego rozkodowania informacji w pliku ISH możemy skorzystać z dostarczonego przez NOAA programu ishJava. Opis krok po kroku jak go uruchomić znajduje się w dokumentacji. Jednak ten sposób jest bardzo niewskazany, ponieważ plik wynikowy z programu ishJava pomimo, że pozwala na wydobycie z pliku ISH wartości 34 głównych parametrów, to wymaga wielu dalszych (i trudnych) edycji, ze względu na znak spacji (czasem wielokrotnej) będący separatorem kolumn, a także na znak “*” oznaczający brak danych, którego ciężko się pozbyć.

O wiele lepszym i efektywnym pomysłem jest skorzystanie z R. Szczególnie dlatego, że do tej pory powstały gotowe opisy oraz pakiety pozwalające w szybki sposób na edycję plików w formacie ISH i przetworzenie ich do uporządkowanej ramki danych. Jeżeli chcielibyśmy sami przetworzyć plik ISH to musimy pamiętać o kolejnych krokach edycji:

odczytanie formatu o stałej szerokości kolumn – tutaj pomocna jest funkcja read.fwf:
- definiujemy szerokość kolumn (na podstawie dokumentacji formatu):
  szerokosci <- c(4, 6, 5, 4, 2, 2, 2, 2, 1, 6, 7, 5, 5, 5,
  4, 3, 1, 1, 4, 1, 5, 1, 1, 1, 6, 1, 1, 1, 5, 1, 5, 1, 5, 1)
- i wczytujemy plik:
  data <- read.fwf("/sciezka i nazwa pliku/"), szerokosci)
wybranie kolumn z danymi, które nas interesują,
nadanie im nagłówków,
przemnożenie wartości przez odpowiednie mnożniki,,
zastąpienie liczb oznaczających braki danych oznaczeniem NA
przekształcenie daty i godziny na odpowiedni format i zebranie tej informacji w jedną kolumnę.

Zainteresowanym polecam opis autorstwa P.L. Delamater, A.O. Finley oraz C. Babcock – Downloading and Processing NOAA hourly weather station data, gdzie cały proces wczytywania i przetwarzania plików ISH w R jest bardzo dobrze opisany z wykorzystaniem pętli dla edycji wielu plików jednocześnie.

2. Pakiet R worldmet

W ostatnich latach powstało co najmniej kilka pakietów R pomocnych w obsłudze plików pochodzących z bazy ISD NOAA, takich jak isdparser czy rnoaa. Jednak na potrzeby dzisiejszych zajęć będziemy pracować z wykorzystaniem pakietu worldmet, który zawiera dosłownie kilka funkcji, jednak bardzo przydatnych. Jego twórcą jest David Carslaw, znany przede wszystkim jako autor pakietu do analiz danych dotyczących jakości powietrza o nazwie openair ³^,⁴ (do którego jeszcze wrócimy).

Poniżej przytoczono dwie najważniejsze funkcje pakietu worldmet, z których będziemy korzystać w trakcie ćwiczeń.

getMeta

Funkcji getMeta będziemy używać do uzyskania metadanych dotyczących stacji. Zapytania możemy tworzyć w oparciu o znaną nazwę stacji, lokalizację geograficzną czy po prostu nazwę kraju. Możemy również na podstawie uzyskanych danych wygenerować mapę z użyciem pakietu leaflet. Domyślna konstrukcja zapytania wygląda w następujący sposób:

getMeta(site = "heathrow", lat = NA, lon = NA, country = NA,
state = NA, n = 10, end.year = "current", plot = TRUE, fresh = TRUE,
returnMap = FALSE)

Przyjrzyjmy się tabeli wynikowej (nie wszystkie z 15 kolumn pokazano) dla 10 stacji znajdujących się najbliżej punktu o współrzędnych geograficznych (50°N, 0°E):

isd <- getMeta(lat=50, lon=0, n=10, plot = FALSE, returnMap = FALSE)

Domyślnie tabela zawiera kolumny informujące nas o kodach (USAF oraz WBAN) i nazwie stacji, jej przestrzennej lokalizacji a także informacje o przedziale czasowym, dla którego w bazie ISD NOAA są dostępne dane.

Ponieważ nasze zapytanie opiera się o konkretną lokalizację geograficzną, to odpowiedź posiada dodatkowe kolumny informujące nas o przesunięciu stacji względem zadanego punktu (longR, latR) oraz o odległości w kilometrach pomiędzy stacją a zadanym punktem (dist). Najważniejsze informacje zawiera kolumna code, ze względu na to, że jest to podstawowy identyfikator stacji wymagany przez funkcję importNOAA (następny rozdział). Warto zauważyć, że kod ten powstaje przez połączenie numerów USAF i WBAN przedzielonych myślnikiem.

importNOAA

To podstawowa funkcja pozwalająca na import danych pomiarowych z bazy ISD NOAA. Jak wspomniano wcześniej, aby skonstruować takie zapytanie, musimy znać kod stacji lub wektor kodów stacji.

Oprócz tego podajemy także informacje dotyczące roku (bądź lat), dla których chcemy pobrać dane. Dodatkowe opcje pozwalają na pobranie 6- lub 12-godzinnego pomiaru opadu (precip) oraz aktualnych warunków pogodowych (PWC). Trzeba równocześnie mieć pewne wyobrażenie na temat spodziewanego rozmiaru danych (w bazie ISD NOAA znajduje się ponad 600GB skompresowanych plików). Konstrukcja funkcji wygląda następująco:

importNOAA(code = "037720-99999", year = 2014, hourly = TRUE, precip = FALSE, PWC = FALSE, parallel = TRUE, quiet = FALSE)

isd_codes <- isd$code  
isd_data <- importNOAA(code=isd_codes, year=2015)

W przypadku braku części danych dostaniemy szczegółowe informacje na ten temat (sterowane przez słowo kluczowe quiet), co ma duże znaczenie przy pobieraniu dużej ilości obserwacji z różnych punktów naraz.

Pobrane dane zawierają dodatkowe parametry, o których jeszcze nie wspominaliśmy, a które są dokładniej opisane w dokumentacji formatu ISH. Chodzi tutaj głównie o widoczność (visibility) i informacje dotyczące zachmurzenia (wysokość i stopień zachmurzenia kolejnych warstw chmur, np. cl_1_height i cl_1).

Oprócz parametrów, które możemy znaleźć domyślnie w bazie ISD NOAA, funkcja importNOAA oblicza również wilgotność względną powietrza (RH) na podstawie danej temperatury powietrza oraz temperatury punktu rosy. Należy zatem pamiętać, że jest to wartość oszacowana a nie zmierzona (rzeczywista).

3. Dane dostępne dla Polski

Sprawdźmy teraz ile danych dla naszego kraju jest dostępnych w bazie ISD NOAA. Historycznie możemy pobrać pomiary meteorologiczne dla 82 różnych lokalizacji w Polsce (głównie stacjonarne stacje klimatologiczne). Skonstruujmy proste zapytanie zwracające położenie tych stacji na interaktywnej leaflet’owej mapie.

getMeta(country="PL", end.year="current", returnMap=TRUE)

Po kliknięciu na dany punkt możemy od razu zobaczyć nazwę stacji, jej kod oraz okres, dla jakiego dostępne są dane pomiarowe. Gdybyśmy ponownie wykorzystali funkcję getMeta używając punktu odniesienia (lat, lon), to mielibyśmy podaną także informację o odległości między tym punktem a stacją w linijce “Distance (km)”.

Pomimo, że czas operacyjny dla niektórych stacji jest imponujący (np. stacja Okęcie realizuje pomiary od 1928 roku), to nie oznacza, że w całym tym okresie zachowana była ciągłość pomiarów. Jest pewien pośredni sposób na sprawdzenie kompletności danych bez konieczności ich pobierania. W tym celu należy pobrać plik isd-inventory z bazy ISD NOAA, gdzie można sprawdzić liczbę rekordów dla danego roku dla danego kodu stacji. Ta metoda nie jest idealna, daje tylko pogląd na temat ilości danych i czasu operacyjnego danej stacji w roku. Problemem tej inwentaryzacji jest fakt, że w bazie ISD NOAA są przechowywane nie tylko dane zakodowane w depeszy SYNOP (standardowe kodowanie powierzchniowych obserwacji meteorologicznych), ale także inne, np. METAR (używany w meteorologii lotniczej i niekiedy wysyłany co 30 minut). To powoduje, że ilość rekordów w danym miesiącu może przekroczyć rzeczywistą liczbę godzin pomiarowych poprzez dublowanie pomiarów różnego typu.

Niemniej jednak, sprawdźmy ile mamy danych dla każdej stacji. Zacznijmy od pobrania i wczytania do R pliku isd-inventory.csv.

download.file("ftp://ftp.ncdc.noaa.gov/pub/data/noaa/isd-inventory.csv", destfile="isd-inventory.csv")
isd_inventory <- read.csv("isd-inventory.csv", header=TRUE, sep=",", dec=".", na.strings="")

Trochę magii z ggplot2 i wykres wynikowy przedstawia się w sposób następujący:

# Uzyskanie metadanych dla stacji z Polski oraz wybranie z utworzonego obiektu kodów USAF
PL_stations=getMeta(country="PL", end.year="current", returnMap=FALSE)
PL_codes <- PL_stations$USAF
# Wybranie z pliku isd-inventory.csv tych stacji, które znajdują się w Polsce
PL_inventory <- isd_inventory[isd_inventory$USAF %in% PL_codes,]

# Obliczenie sumy rekordów dla danego roku
PL_inventory$ALL <- rowSums(PL_inventory[4:15])

# Wybranie odpowiednich kolumn z obiektu i połączenie z wynikami 
# na temat ilości pomiarów dla stacji w Polsce
PL_names <- subset(PL_stations, select = c(USAF, STATION, latitude, longitude, BEGIN, END))
PL_inventory <- merge(PL_inventory, PL_names)

# Obliczenie sumy pomiarów i ilości stacji, dla których są pomiary w danym roku
PL_inventory %>% 
  group_by(YEAR) %>% 
  summarise(suma_pomiarow = sum(ALL),
            suma_stacji = length(unique(STATION))) -> PL_inv_YEAR

# Dodanie informacji koniecznych do wyświetlania wykresu
year <- seq(1928, 2019, by = 1)
year <- as.data.frame(year)

PL_inv_YEAR <- merge(year, PL_inv_YEAR, by.x = "year", by.y = "YEAR", all.x=TRUE)

czy_pomiary <- data.frame(ymin = -Inf,
                          ymax = Inf,
                          xmin = c(1920,   1943.5,  1963.5),
                          xmax = c(1927.5, 1951.5,  1972.5))

# Wyłączenie notacji naukowej
options("scipen"=100, "digits"=2)

y_labs <- format(seq(0, 1200, by = 100), big.mark = ' ')

# Utworzenie wykresu
ggplot(PL_inv_YEAR, aes(y=suma_pomiarow, x=year)) +
  geom_bar(stat='identity', aes(fill=year), show.legend = FALSE) +
  scale_y_continuous(name = 'Liczba rekordów pomiarowych [w tys.]',
                     limits = c(0, 1200000),
                     expand = c(0, 0),
                     breaks = seq(0, 1200000, by = 100000),
                     labels = y_labs,  
                     sec.axis=sec_axis(~./10000, 
                                       name='Liczba stacji',
                                       breaks = seq(0, 120, by = 10))) +
  scale_x_continuous(name = 'Rok',
                     limits = c(1920, 2020),
                     expand = c(0, 0),
                     breaks = seq(1920, 2020, by = 10)) +
  geom_line(aes(y=suma_stacji*10000), color = 'darkgreen', size = 1) +
  geom_point(aes(y=suma_stacji*10000), color= 'darkgreen', size = 1.5) +
  geom_rect(data = czy_pomiary, 
            aes(ymin = ymin, ymax = ymax,
                xmin = xmin, xmax = xmax), fill = 'red',
            alpha = 0.15, inherit.aes = FALSE) +
  theme_minimal() +
  labs(title = 'Liczba pomiarów i stacji dla Polski w bazie ISD NOAA', 
        subtitle = paste0('Stan na ',format(Sys.time(),'%d %B, %Y')))

Sumaryczna ilość pomiarów dla wszystkich stacji ukazana jest za pomocą niebieskich słupków, czerwone obszary oznaczają kompletny brak danych w pewnym okresie, a zielona linia pokazuje liczbę stacji, dla których są dostępne pomiary w danym roku.

Można na tej podstawie wywnioskować, że ilość stacji w Polsce raportujących pomiary do bazy ISD NOAA wciąż rośnie, a tym samym – dostępnych pomiarów jest coraz więcej.

Kiedy spojrzymy jak to rozkłada się według stacji (tabela poniżej), to możemy zauważyć, że najwięcej pomiarów dostępnych jest dla stacji meteorologicznych zlokalizowanych w bezpośredniej okolicy lotnisk. Warto również zauważyć, że większość stacji raportujących pomiary wciąż jest aktywna.

Kiedy już wiemy, czego możemy spodziewać się w bazie ISD NOAA oraz jak wygląda struktura plików, to nie pozostaje nic innego jak tylko przejść do analizy danych.

4. Przykładowa analiza statystyczna

Wyniki analiz danych klimatologicznych wykorzystywane są do wielu celów, m. in.:

do sporządzania prognoz krótkoterminowych,
do oceny zmian klimatu na podstawie wieloletnich danych (ostatnio modne),
w innych dziedzinach (np. w ochronie jakości powietrza czy gospodarce wodnej).

Dzisiaj zajmiemy się najbardziej ogólną analizą statystyczną dla jednej wybranej stacji pomiarowej. Aby zapewnić jak najdłuższy okres kompletnych danych wybierzmy stację na lotnisku Okęcie w Warszawie.

okecie <- importNOAA(code="123750-99999", year=seq(1928,2019, by=1))

Wiemy już, że dla wszystkich stacji w Polsce występują okresowe braki danych. Przyjrzyjmy się jak to wygląda dla stacji Okęcie:

# Dodanie informacji o roku do ramki danych
okecie$years <- format(okecie$date, "%Y")

# Policzenie ilości pomiarów dla danego roku dla stacji Okęcie
okecie %>%
  group_by(years) %>%
  summarise(liczba_pomiarow = n()) -> okecie_summary

# Uzupełnienie lat, dla których nie ma pomiarów dla stacji Okęcie
okecie_summary <- merge(year, okecie_summary, by.x = "year", by.y = "years", all.x=TRUE)

# Dodanie informacji do wykreślenia wykresu
czy_pomiary_okecie <- data.frame(ymin = -Inf,
                          ymax = Inf,
                          xmin = c(1920,   1939.5, 1941.5,  1964.5),
                          xmax = c(1927.5, 1940.5, 1951.5,  1972.5))

# Utworzenie wykresu
ggplot(okecie_summary, aes(x=year, y=liczba_pomiarow, color=liczba_pomiarow)) +
  geom_point(size=3, show.legend = FALSE) + 
  scale_y_continuous(name = 'Liczba rekordów pomiarowych',
                     limits = c(0, 10000),
                     expand = c(0, 0)) +
  scale_x_continuous(name = 'Rok',
                     limits = c(1920, 2020),
                     expand = c(0, 0),
                     breaks = seq(1920, 2020, by = 10)) +
  geom_rect(data = czy_pomiary_okecie, 
            aes(ymin = ymin, ymax = ymax,
                xmin = xmin, xmax = xmax), fill = 'red',
            alpha = 0.15, inherit.aes = FALSE) +
  theme_minimal() +
  theme(plot.margin = unit(c(.5,.5,.5,.5), "cm")) +
  labs(title = 'Liczba pomiarów dla stacji pomiarowej Okęcie (PL 123750-99999) w bazie ISD NOAA', 
       subtitle = paste0('Stan na ',format(Sys.time(),'%d %B, %Y')))

Poza niektórymi latami stacja Okęcie odznacza się wysoką kompletnością danych (ponad 8000 rekordów pomiarowych dla większości lat). Aby zachować ciągłość obserwacji skupmy się na latach 1973-2019. Wybierzmy również tylko te parametry, na których skoncentruje się analiza statystyczna – data oraz najważniejsze pomiary meteorologiczne.

okecie_47 <- filter(okecie, years > 1972 & years < 2020)

okecie_47s <- subset(okecie_47, select = c(1, 9, 10, 11, 12, 13, 14, 15, 16))

W kolejnych analizach, ze względu na ilość rekordów danych, będziemy posługiwać się raczej wynikami w formie graficznej niż tabelarycznej – łatwiejsze do interpretacji.

Statystyki opisowe

Pierwszym krokiem w statystycznej analizie danych jest przyjrzenie się podstawowym statystykom opisowym w celu podsumowania zbioru danych i wyciągnięcia ogólnych wniosków na tej podstawie⁵. Tak też to wygląda w analizie danych meteorologicznych. Na potrzeby dzisiejszego ćwiczenia ograniczmy się jedynie do:

ilości obserwacji – n(x)
średniej arytmetycznej – mean(x)
odchylenia standardowego – sd(x)
minimum oraz maksimum – min(x) oraz max(x)
amplitudy – max(x) - min(x)
kwantyli 1/4 oraz 3/4 – quantile(x, c(.25, .75))

Na szybko możemy posłużyć się funkcją summary(x) lub stat.desc(x) z pakietu pastecs, które zwracają więcej statystyk. Wykorzystamy ostatnią możliwość, jako że w odpowiedzi dostajemy elegancko uporządkowaną ramkę danych.

# Wczytanie biblioteki i utworzenie obiektu z podsumowaniem statystyk opisowych dla ramki danych
library(pastecs)
desc_okecie <- stat.desc(okecie_47s)

# Wybranie interesujących nas statystyk i parametrów
desc_okecie <- desc_okecie[c(1,9,13,4, 5, 6), c(2:9)] 

# Ręczne policzenie kwantyli
q_wd <- quantile(okecie_47s$wd, c(1/4, 3/4), na.rm=TRUE)
q_ws <- quantile(okecie_47s$ws, c(1/4, 3/4), na.rm=TRUE)
q_hgt <-quantile(okecie_47s$ceil_hgt, c(1/4, 3/4), na.rm=TRUE)
q_vis <-quantile(okecie_47s$visibility, c(1/4, 3/4), na.rm=TRUE)
q_temp <-quantile(okecie_47s$air_temp, c(1/4, 3/4), na.rm=TRUE)
q_dew <-quantile(okecie_47s$dew_point, c(1/4, 3/4), na.rm=TRUE)
q_pres <-quantile(okecie_47s$atmos_pres, c(1/4, 3/4), na.rm=TRUE)
q_rh <-quantile(okecie_47s$RH, c(1/4, 3/4), na.rm=TRUE)

# Dodanie informacji o kwantylach do jednego obiektu
desc_okecie_q <- cbind(q_wd, q_ws, q_hgt, q_vis, q_temp, q_dew, q_pres, q_rh)

# Nadanie tych samych nazw kolumn statystykom opisowym i kwantylom
colnames(desc_okecie) <- c("Kierunek wiatru","Prędkość wiatru","Podstawa chmur", "Widoczność",
                           "Temperatura", "Punkt rosy", "Ciśnienie atm.", 
                           "Wilgotność względna")

colnames(desc_okecie_q) <- c("Kierunek wiatru","Prędkość wiatru","Podstawa chmur", "Widoczność",
                           "Temperatura", "Punkt rosy", "Ciśnienie atm.", 
                           "Wilgotność względna")

# Połączenie statystyk opisowych wraz z kwantylami
desc_okecie <- rbind(desc_okecie, desc_okecie_q)

# Nadanie nazw statystykom opisowym
rownames(desc_okecie) <- c("Liczba", "Średnia", "Odchylenie st.", "Minimum", "Maksimum", 
                           "Amplituda", "Kwantyl 1/4", "Kwantyl 3/4")

# Transpozycja ramki danych
desc_okecie <- as.data.frame(t(desc_okecie))

Na podstawie samej tabeli ze statystykami opisowymi ciężko wnioskować o generalnych trendach poszczególnych parametrów. Niemniej jednak w przypadku temperatury czy prędkości wiatru te informacje są przydatne. Czasem natomiast też znajdziemy tutaj wartości “spodziewane” – np. że minimalny kierunek wiatru wynosi 0, a jego amplituda 360 stopni kątowych. Wartościowe są też informacje o liczebności danych. Możemy zauważyć, że pomiarów dotyczących ciśnienia atmosferycznego jest zdecydowanie mniej – to wynika z tego, że na początku obserwacje te były raportowane co kilka godzin (to się zmieniło dopiero w 1999 roku). Podobnie było z wysokością podstawy chmur.

Niekiedy warto przytaczać osobno wartości statystyk opisowych dla poszczególnych lat pomiarowych lub pór roku, aby wydobyć więcej informacji z zestawu danych i sformułować dodatkowe wnioski.

Przydatne są też metody graficzne, za pomocą których można zwizualizować rozkłady zmiennych lub nawet częściowo zależności między nimi i ich zmienność czasową (np. za pomocą wykresów pudełkowych). Przygotujmy kilka wykresów na przykładzie temperatury powietrza.

Histogram

Histogram pokazuje jak rozkładają się wartości danej zmiennej. Bardzo ciekawie to będzie wyglądać w przypadku danych meteorologicznych z podziałem np. na pory roku. Od tego momentu będziemy wspomagać się pakietem openair, który umożliwia, między innymi, szybkie dodanie informacji o porach roku do naszych danych, a także dodatkową wizualizację danych.

Dodajmy do naszego zbioru danych informację o porze roku za pomocą funkcji cutData, a potem skorzystajmy z funkcji hist z pakietu graphics do wyrysowania rozkładu temperatur powietrza w zimie oraz w lecie.

# Wczytanie biblioteki i dodanie informacji o porach roku do ramki danych
library(openair)
okecie_47s <- cutData(okecie_47s, type="season")

# Wykreślenie kolejno dwóch histogramów nakładających się na siebie dla lata i zimy
hist(okecie_47s$air_temp[okecie_47s$season == "summer (JJA)"], breaks=10, xlim=c(-50,50), 
     col=rgb(1,0,0,0.5), xlab="Temperatura powietrza [°C]", ylab="Częstość występowania", 
     main="Rozkład temperatury powietrza w latach 1973-2019 na stacji Okęcie\nz podziałem na lato i zimę")
hist(okecie_47s$air_temp[okecie_47s$season == "winter (DJF)"], breaks=10, xlim=c(-50,50), 
     col=rgb(0,0,1,0.5), xlab="Temperatura powietrza [°C]", ylab="Częstość występowania", 
     main="Rozkład temperatury powietrza w latach 1973-2019 na stacji Okęcie\nz podziałem na lato i zimę", add = T)
legend("topright", legend=c("Lato", "Zima"), col=c(rgb(1,0,0,0.5), 
                                                      rgb(0,0,1,0.5)), pt.cex=2, pch=15)

Otrzymane rozkłady są bardzo czytelne, widać wyraźnie różnicę między temperaturami na stacji Okęcie rejestrowanymi w zimie oraz w lecie. Możemy zauważyć, w lecie dominowały temperatury poniżej 25°C, a w zimie – do 5°C powyżej zera.

Wykres pudełkowy

Dzięki wykresom pudełkowym możemy zwizualizować wartości najważniejszych statystyk i łatwo porównać różne zmienne bądź jedną zmienną podzieloną według zadanego warunku. Inna ich nazwa to wykresy “ramka-wąsy”, bo dokładnie tak wyglądają.

Sporządźmy wykres pudełkowy dla pomiarów temperatury z podziałem na różne pory roku w analizowanym okresie.

ggplot(okecie_47s, aes(x=season, y=air_temp, fill=season)) + 
  geom_boxplot(alpha=0.5) +
  stat_boxplot(geom ='errorbar',width = 0.2) + 
  stat_summary(fun.y=mean, geom="point", shape=18, size=4, color="black", fill="red") +
  theme_minimal() +
  theme(legend.position="none", 
        plot.title = element_text(size=14, face="bold", hjust = 0.5), 
        axis.title=element_text(size=12),
        axis.text=element_text(size=12, color="black")) +
  labs(x = "Pora roku",
       y = "Temperatura powietrza [°C]",
       title = "Wykresy pudełkowe dla temperatury powietrza w latach 1973-2019 na stacji Okęcie\n
       z podziałem na pory roku")

Pozioma kreska w “pudełku” to mediana, a czarny romb symbolizuje średnią wartość danych – dodano go ręcznie, aby pokazać, że te dwie wartości nie są tym samym. Końce pudełka są wyznaczone przez tzw. rozstęp kwartylowy (IQR), a więc różnicę między kwantylem rzędu 3/4 (Q3) a 1/4 (Q1). W pudełku znajduje się zatem połowa wszystkich obserwacji, a poza nią pozostałe wartości. Czarnymi kropkami oznaczono wartości odstające, które są:

mniejsze niż Q1 - 1.5 * IQR,
większe niż Q3 + 1.5 * IQR.

Z powyższego wykresu możemy wywnioskować, że, rzecz jasna, najwyższe temperatury obserwowane są w lecie a najmniejsze w zimie. Co ciekawe, temperatury na jesień są nieco wyższe niż te na wiosnę, chociaż intuicyjnie wydaje się być na odwrót. Warto tutaj podkreślić, że podział na pory roku jest umowny i wyznaczany na podstawie pełnych miesięcy, np. wiosna to wartości z marca, kwietnia oraz maja. Duża ilość obserwacji odstających wynika z charakteru danych meteorologicznych i z faktu, że są to dane godzinowe o wysokich wahaniach. Prawdopodobnie średnie dobowe skupiałyby się bardziej wokół mediany.

Wykres z linią trendu

Trendy w serii czasowej danych (time series) można przedstawić na kilka sposobów. Dzisiaj nauczymy się sporządzać wykresy przedstawiające tzw. średnią kroczącą. Jest to średnia wartość obserwacji dla danej liczby okresów (dla nas będą to kolejne lata). Jest ona przypisywana do ostatniego roku branego pod uwagę w obliczeniach. Pozwala ona również na wygładzanie zmienności danych o tak dużych wahaniach jak temperatura powietrza. Możemy ją obliczyć wykorzystując funkcję ma z pakietu forecast.

Z powodu dużej zmienności temperatury (w ciągu dnia i miesiący) na samym początku zagregujemy godzinowe wartości temperatury powietrza do średnich miesięcznych, a następnie zamienimy ramkę danych na serię czasową za pomocą funkcji ts z pakietu stats.

Teraz sporządzimy trzy wykresy obrazujące średnie miesięczne temperatury powietrza w badanym okresie oraz dodamy kolorową linię ze średnią kroczącą – 1-roczną, 5-letnią oraz 10-letnią.

# Wczytanie biblioteki "forecast" oraz dodanie kolumn z informacją o roku i miesiącu do ramki danych.
library(forecast)
okecie_47s <- cutData(okecie_47s, type="year")
okecie_47s <- cutData(okecie_47s, type="month")

# Obliczenie średniej miesięcznej temperatury dla każdego miesiąca w roku i zapisanie jej do nowego obiektu.
okecie_47s %>% select(air_temp, year, month) %>%
  group_by(year, month) %>%
  summarise(mean_temp = mean(air_temp, na.rm=TRUE)) -> t_okecie_mean

# Zamiana średnich temperatur na serię czasową
t_okecie_ts <- ts(t_okecie_mean[3], frequency = 12, start = c(1973,1), end = c(2019,12))

# Ustawienie wielu wykresów na jednym rysunku
par(mfrow = c(3,1))

# Wykreślenie średnich miesięcznych temperatury wraz ze średnią kroczącą dla 1 roku, 5 lat oraz 10 lat
plot.ts(t_okecie_ts, col="darkgray", main = "1-roczna średnia krocząca temperatury powietrza", 
        xlab="Rok", ylab="Temperatura powietrza [°C]")
lines(ma(t_okecie_ts, order = 12), col = "red", lwd=3)
grid(col = "lightgray", lty = "dotted",
     lwd = par("lwd"), equilogs = TRUE)

plot.ts(t_okecie_ts, col="darkgray", main = "5-letnia średnia krocząca temperatury powietrza",
        xlab="Rok", ylab="Temperatura powietrza [°C]")
lines(ma(t_okecie_ts, order = 60), col = "orange", lwd=3)
grid(col = "lightgray", lty = "dotted",
     lwd = par("lwd"), equilogs = TRUE)

plot.ts(t_okecie_ts, col="darkgray", main = "10-letnia średnia krocząca temperatury powietrza",
        xlab="Rok", ylab="Temperatura powietrza [°C]")
lines(ma(t_okecie_ts, order = 120), col = "darkgreen", lwd=3)
grid(col = "lightgray", lty = "dotted",
     lwd = par("lwd"), equilogs = TRUE)

Powyższe wykresy ujawniają, że dla stacji Okęcie w ciągu ostatnich 24 lat istnieje niewielki, jednak zauważalny dodatni trend średniej rocznej temperatury powietrza (im lepsze wygładzenie średniej kroczącej tym łatwiej to zauważyć). Patrząc na średnie miesięczne można również stwierdzić, że w miesiącach letnich średnie temperatury powietrza coraz częściej przekraczają 20°C w ostatnich latach.

Róża wiatrów

Róża wiatrów to powszechnie stosowany wykres obrazujący w sposób przejrzysty częstość wystąpień (w procentach) wiatru o danej prędkości z danego sektora kierunku wiatru. Może on przybierać nieco zróżnicowane formy, ale zawsze jest łatwy do interpretacji. Żeby przyspieszyć tworzenie tego wykresu skorzystamy z pakietu openair i funkcji windRose. Pakiet openair jest bardzo obszerny i oprócz dokumentacji dostępnej na CRAN’ie, autor opracował również wyczerpujący manual zawierający wiele przykładów zastosowań poszczególnych funkcji.

Dosyć częstą wadą wykresów, które są dostępne w niektórych bibliotekach jest pewna ograniczona możliwość ingerencji w jego ostateczny wygląd. Na szczęście wykresy z openaira są dosyć schludne i można je wykorzystać bez konieczności wprowadzania dużych modyfikacji.

Wróćmy do wykresu róży wiatrów. Sporządzimy go dzieląc dane na cztery pory roku w celu potencjalnego ujawnienia różnic w sezonach.

windRose(okecie_47s, ws="ws", wd="wd", type="season", paddle="false", layout=c(4,1),
         par.settings=list(fontsize=list(text=18)), 
         key = list(labels = c("0-2", "2-4", "4-6", "6-24")),
         key.footer = "(m/s)",
         statistic =  list("fun"=length,"unit" = "%","scale" = "all", 
                           "lab" = "Częstość wystąpień wiatru z danego kierunku [%]", 
                           "fun2" = function(x) signif(mean(x, na.rm = TRUE), 3), 
                           "lab2" = "mean","labcalm" = function(x) round(x, 1)))

Na pierwszy rzut oka można zauważyć, że dla stacji Okęcie przez cały rok utrzymuje się dominujący kierunek wiatru z zachodu (typowe dla Polski), częściowo też z południowego wschodu (oprócz lata). Widoczne są też różnice w średnich prędkościach wiatru (mean pod wykresem) oraz częstotliwością wystąpień cisz wiatrowych (calm). W lecie średnia prędkość wiatru jest najniższa i wynosi 3.35 m/s (przy 7.4% ciszy wiatrowej), natomiast w zimie jest najwyższa – 4.29 m/s (3.1% ciszy wiatrowej).

Zależność liniowa

Przejdźmy do pobieżnego sprawdzenia, czy pomiędzy parametrami występują zależności. Posłużymy się najbardziej popularnym sposobem – obliczenia ile wynosi współczynnik korelacji liniowej Pearsona.

Macierz korelacji

Macierz korelacji liniowej Pearsona pozwala wstępnie sprawdzić występowanie zależności między analizowanymi zmiennymi. Mówiąc prosto, współczynnik ten może przyjmować wartości od -1 do 1. Im jest mniejszy, tym silniejsza zależność ujemna między parametrami, a im jest bliższy 1, tym silniejsza zależność dodatnia. Dodatkowo możemy sprawdzić, czy otrzymane wartości są istotne statystycznie (przyjmiemy poziom istotności równy 0.05).

Korzystamy z funkcji corrplot z pakietu corrplot do wyrysowania macierzy korelacji. Będziemy też potrzebować pakietu Hmisc do policzenia współczynników.

# Wczytanie odpowiednich bibliotek
library(corrplot)
library(Hmisc)

# Wybranie interesujących nas kolumn z zestawu danych
okecie_47s_d <- okecie_47s[2:9]

# Nadanie kolumnom odpowiednich nazw, które będą wyświetlane na wykresie
colnames(okecie_47s_d) <- c("Kierunek wiatru","Prędkość wiatru","Podstawa chmur", "Widoczność",
                           "Temperatura", "Punkt rosy", "Ciśnienie atm.", 
                           "Wilgotność względna")
# Definiowanie indywidualnej palety kolorów w zależności od współczynnika korelacji (kodowanie hex)
col <- colorRampPalette(c("#BB4444", "#EE9988", "#FFFFFF", "#77AADD", "#4477AA"))

# Policzenie odpowiednich współczynników za pomocą funkcji 'corr'
corr <- rcorr(as.matrix(okecie_47s_d))

# Zdefiniowanie wektora ze współczynnikami korelacji
corr_r <- corr$r

# zdefiniowanie wektora z wartościami p-value
corr_p <- corr$P

# Wykreślenie macierzy korelacji
corrplot(corr_r, method = "color", col = col(200),  
         type = "upper", order = "hclust", addCoef.col = "black", diag = FALSE,
         tl.col = "black", tl.srt = 45, p.mat = corr_p, sig.level = 0.05)

To co należy na wstępie zauważyć to fakt, że wszystkie wartości współczynnika korelacji są statystycznie istotne na przyjętym poziomie istotności α = 0.05) – gdyby nie były, to byłyby przekreślone na wykresie. Druga najważniejsza obserwacja dotyczy bardzo silnej dodatniej korelacji między temperaturą powietrza a temperaturą punktu rosy. Jest to akurat rzecz naturalna, wynikająca z definicji i zależności fizycznych tych dwóch parametrów - teoretycznie temperatura punktu rosy nie może wynosić więcej niż temperatura powietrza (ten stan przekraczałby wartość 100% wilgotności względnej). Kolejne ciekawe wnioski dotyczą widoczności - jest ona większa przy wyższych temperaturach, wyżej położonej podstawie chmur i mniejszej wilgotności powietrza (a więc w warunkach niesprzyjających tworzeniu się mgieł).

Natomiast sceptycznie i z ograniczonym zaufaniem należy podchodzić do wyników dotyczących wilgotności względnej – przypominam, że jest to wartość obliczona na podstawie temperatury powietrza i temperatury punktu rosy, jasne jest zatem, że będzie się charakteryzować zależnością od tych parametrów.

Wykres rozrzutu

Wykres rozrzutu ukazuje zależności między dwoma parametrami, umieszczając je kolejno na osi X oraz Y. Z ułożenia punktów możemy wstępnie wnioskować o występowaniu zależności a także o jej postaci. Pokażmy dwa wykresu rozrzutu ilustrujące zależności, które zidentyfikowaliśmy w macierzy korelacji wcześniej – między temperaturą powietrza a temperaturą punktu rosy. Przedstawimy go za pomocą metody hexbin, która grupuje punkty i liczy ich ilość. Wykonamy też drugi wykres rozrzutu między trzema zmiennymi, gdzie trzecia zmienna – widoczność – będzie przedstawiona na wykresie za pomocą skali kolorystycznej.

# Wybranie odpowiednich kolumn
okecie_47s_d <- okecie_47s[2:9]

# wykreślenie wykresu rozrzutu typu 'hexbin'
scatterPlot(okecie_47s, x = "air_temp", y = "dew_point", method = "hexbin", col= "Blues", 
            xlab = "Temperatura powietrza [°C]", main="Wykres rozrzutu z liczbą wystąpień", 
            key.footer = "Liczba",
            ylab = "Temperatura punktu rosy [°C]")

Powyższy wykres obrazuje bardzo wyraźnie to co wynika z zależności fizycznych temperatury powietrza i temperatury punktu rosy. Można zauważyć, że w istocie temperatura punktu rosy nie jest większa niż temperatura powietrza, ale może być niższa. Z teorii: im większa różnica między tymi parametrami, tym wilgotność względna powietrza jest niższa. Wyraźnie widoczny jest stożkowy kształt punktów na wykresie, sugerujący, że przy wyższych temperaturach zachodzi większa rozbieżność między analizowanymi parametrami. Podpierając się teorią można powiedzieć, że przy wyższych temperaturach powietrza częściej zdarza się, że wilgotność względna powietrza jest niższa. I ten fakt jest również potwierdzony przez wartość współczynnika korelacji między temperaturą powietrza a wilgotnością względną (-0.51).

scatterPlot(okecie_47s, x = "air_temp", y = "RH", z="visibility", col= "Purples", 
            xlab = "Temperatura powietrza [°C]", main = "Wykres rozrzutu z trzecią zmienną", 
            key.footer="Widoczność\n[m]",
            ylab = "Wilgotność względna [%]")

Wykresy rozrzutu można wzbogacić o trzecią zmienną. Oprócz zależności między temperaturą powietrza a wilgotnością względną przedstawiono również widoczność za pomocą koloru fioletowego. Wprost z wykresu wynika, że z reguły większa widoczność związana jest z wyższymi temperaturami i niższą wilgotnością względną powietrza – dosyć oczywisty wniosek.

5. Projekt do wykonania

Wybierz stację pomiarową (niekoniecznie z Polski), zlokalizuj ją na mapie i opisz jej położenie wspomagając się metadanymi z bazy ISD NOAA.
Dla wybranej stacji sprawdź kompletność danych w bazie ISD NOAA w kolejnych latach (1901-2019).
Wybierz dziesięciolecie charakteryzujące się największą kompletnością danych, pobierz je, a następnie na ich podstawie sporządź:

Podsumowanie ze statystykami najważniejszych parametrów w kolejnych latach (średnia, minimalna i maksymalna wartość, wartości kwartyli, odchylenie standardowe i amplituda),
Sezonowe róże wiatrów,
Wykres zmienności temperatury (na podstawie średnich rocznych) z linią trendu.

Sprawdź czy istnieją istotne zalezności liniowe między parametrami (za pomocą współczynnika korelacji Pearsona) i sporządź wykres rozrzutu dla najsilniej skorelowanych parametrów.
Skomentuj uzyskane wyniki i sporządź odpowiednie wnioski.

Gotowy projekt w formacie *.pdf wraz ze skryptem *.R należy przesłać w ciągu dwóch tygodni od przeprowadzonych zajęć na adres mailowy prowadzącego.

S.A. del Greco, J.N. Lott, K. Hawkins, et al.: Surface data integration at NOAA’s National Climatic Data Center: Data format, processing QC, and product generation↩
J.N. Lott: The Quality Control of the Integrated Surface Hourly database↩
D.C. Carslaw, K. Ropkins (2012): openair — An R package for air quality data analysis, Environmental Modelling & Software, 27–28(0), 52–61. DOI: 10.1016/j.envsoft.2011.09.008 ↩
A. Szulecka, R. Oleniacz, M. Rzeszutek (2017): Functionality of openair package in air pollution assessment and modeling – a case study of Krakow, Environmental Protection and Natural Resources, 28(2), 22-27. DOI: 10.1515/OSZN-2017-0009 ↩
D.B. Stephenson (2005): Data analysis methods in weather and climate research, on-line course.↩