A lato było cieplejsze...
Czy w Boże Narodzenie zawsze brakowało śniegu? Czy klimat się ociepla? Jak ciepło przepływa przez Polskę? Dzisiaj zajmiemy się pogodą. Sporo ciekawych zagadnień przed nami:
- Skąd wziąć dane o pogodzie w Polsce?
- Pory roku na przestrzeni lat - czy wiosna przychodzi szybciej? A może lat trwa krócej?
- zobaczymy procentowy podział dni w roku na poszczególne pory roku (ile było dni wiosennych, a ile zimowych)
- czy wiosna przychodzi coraz szybciej - jak kształtowała się temperatura w marcu na przestrzeni lat?
- Czy klimat się ociepla?
- jak wyglądały różnice w temperaturze na przestrzeni lat
- czy poszczególne miesiące są coraz cieplejsze lub zimniejsze?
- podobnie dla kolejnych dni - jak wyglądała liczba dni cieplejszych i zimniejszych niż średnia?
- Czy jest różnica miedzy temperaturami w kraju?
- przygotujemy animowane mapki temperatury dla 2017 roku
- A kiedyś to było…
- Boże Narodzenie ze śniegiem
- w Lany Poniedziałek było cieplej
Pobranie danych
Zanim przejdziemy do analizy potrzebujemy danych o pogodzie. Zgodnie z ustawą o dostępie do danych publicznych IMGW udostępnia dane historycznie na stronie https://dane.imgw.pl/. Dane aktualne są dostępne poprzez API, ale dane historyczne spakowane już niestety nie. Dostępne są w postaci spakowanych plików CSV, podzielone na kilkanaście katalogów. Do wszystkiego da się jednak dobrać. Aby skorzystać z tych zasobów należy każdy plik (archiwum) pobrać, rozpakować i wyciągnąć z niego dane. Aby nie operować na dużej liczbie plików przepiszemy je do bazy danych (w tym przypadku będzie to plikowa baza SQLite, chociaż może być dowolna - mySQL, PostgreSQL). Skrypt przepisujący dane znajdziecie na moim githubie (jak i pozostałe poniżej zaprezentowane kody).
Przegląd danych
Zaczniemy od wczytania potrzebnych pakietów i podłączeniu się do bazy danych ze zgromadzonymi informacjami:
library(tidyverse)
library(lubridate)
library(ggridges)
library(ggrepel)
library(forcats)
library(DBI)
library(jsonlite)
library(kknn)
dbase <- dbConnect(RSQLite::SQLite(), "imgw.sqlite")
Teraz możemy zacząć pracę.
średnia temperatura dzienna ze wszystkich stacji, dzień po dniu
Zaczniemy od pobrania temperatury (uśrednionej po wszystkich stacjach z danego dnia).
mTemps <- dbGetQuery(dbase,
"SELECT Rok, Miesiac, Dzien, AVG(MeanTemp) AS mTemp
FROM imgw
WHERE Rok <> 2018
GROUP BY Rok, Miesiac, Dzien;") %>%
# robimy sobie pole z datą
mutate(data = make_date(Rok, Miesiac, Dzien)) %>%
# średnia temperatura dzienna ze wszystkich lat
group_by(Miesiac, Dzien) %>%
mutate(mTemp_all = mean(mTemp)) %>%
ungroup() %>%
# różnica pomiędzy temperaturą danego dnia a średnią ze wszystkich lat tego dnia
# dodatnie = cieplej niż średnia
mutate(dTemp = mTemp - mTemp_all) %>%
# w zależności od średniej dziennej temperatury przypisujemy porę roku
mutate(pora_roku = case_when(
.$mTemp >= 15 ~ "lato",
.$mTemp >= 5 ~ "wiosna/jesień",
.$mTemp >= 0 ~ "przedwiośnie/przedzimie",
.$mTemp < 0 ~ "zima",
TRUE ~ "coś innego" # to nie ma prawa się trafić :)
))
Jak widzicie z bazy wyciągamy dane odpowiednio skonstruowanym zapytaniem SQL. To bardziej efektywne niż wciągnięcie całych danych do pamięci w R i ich obróbka funkcjami z pakietu dplyr. Pakiet jest świetny, ale serwer bazodanowy jest bardziej wydajny. Przy danych rzędu kilkuset tysięcy wierszy nie ma to większego znaczenia, jednak kiedy mamy wiersze idące w miliony i przy okazji kilkadziesiąt kolumn to robi się problem (w R). Po pobraniu danych do pamięci w powyższym kodzie budujemy sobie datę (przyda się gdzieś później), liczymy średnią temperaturę dla każdego dnia po wszystkich latach (czyli średnią dla na przykład 20 marca bez względu na rok). Dodatkowo przypisujemy klimatyczne pory roku według średniej dziennej temperatury:
- temperatura poniżej zera stopni Celsiusza = zima
- temperatura 0-5°C - przedwiośnie lub przedzimie
- 5-15°C = wiosna lub jesień
- temperatura powyżej 15°C to lato
Liczymy też odchylenie temperatury danego dnia od temperatury uśrednionej po wszystkich latach. Tym razem przygotowanie danych poszło w miarę szybko, zobaczmy jak wygląda
pora roku na przestrzeni lat
# kolorki
pory_roku_palette <- c("lato" = "#d7191c", "wiosna/jesień" = "#fdae61", "przedwiośnie/przedzimie" = "#abd9e9", "zima" = "#2c7bb6")
mTemps %>%
ggplot() +
geom_tile(aes(make_date(2000, month(data), day(data)), Rok, fill = pora_roku)) +
scale_fill_manual(values = pory_roku_palette) +
scale_y_reverse() +
scale_x_date(date_labels = "%d/%m", date_breaks = "1 month") +
theme(legend.position = "bottom") +
labs(title = "Pora roku określona na podstawie temperatury dobowej",
x = "Dzień w roku", y = "Rok",
fill = "Pora roku")
Widzimy, że lato (czerwone punkty; właściwie letnie dni) z roku na rok zaczyna się delikatnie coraz szybciej - różnica jest na poziomie pojedynczych dni. Bywały lata kiedy ponad 15°C było już na początku maja (jakoś około 1993 roku).
Ciekawostka - biała linia przy 1 marca to oczywiście 29 lutego.
W każdym razie widać tutaj spory szum - uśrednijmy te dane na poziomie tygodni:
mTemps %>%
mutate(week = week(data)) %>%
group_by(week, Rok) %>%
summarise(mTemp = mean(mTemp)) %>%
ungroup() %>%
mutate(pora_roku = case_when(
.$mTemp >= 15 ~ "lato",
.$mTemp >= 5 ~ "wiosna/jesień",
.$mTemp >= 0 ~ "przedwiośnie/przedzimie",
.$mTemp < 0 ~ "zima",
TRUE ~ "cos innego"
)) %>%
ggplot() +
geom_tile(aes(week, Rok, fill = pora_roku)) +
scale_fill_manual(values = pory_roku_palette) +
scale_y_reverse() +
scale_x_continuous(breaks = seq(1, 55, 5)) +
theme(legend.position = "bottom") +
labs(title = "Pora roku określona na podstawie temperatury dobowej",
x = "Tydzień w roku", y = "Rok",
fill = "Pora roku")
Około 1990 roku mieliśmy długą wiosnę (od jakiegoś 11 do 23 tygodnia). Czy jest coraz cieplej? W odpowiedzi na to pytanie może pomóc procentowy podział dni w roku na poszczególne pory roku:
mTemps %>%
group_by(Rok, pora_roku) %>%
summarise(n = n()) %>%
ungroup() %>%
group_by(Rok) %>%
mutate(p = 100*n/sum(n)) %>%
ungroup() %>%
filter(Rok != 2018) %>%
ggplot() +
geom_smooth(aes(Rok, p, color = pora_roku), se = FALSE) +
theme(legend.position = "bottom") +
labs(title = "Udział procentowy dni z określoną (na podstawie średniej\ndobowej temperatury) porą roku na przestrzeni lat 1951-2017",
y = "% dni w danej porze roku", x = "", color = "Pora roku")
Widzimy, że linia fioletowa (odpowiedzialna za udział dni ze średnią temperaturą dobową poniżej 0°C) spada na przestrzeni lat z około 20% (1/5 roku) do 15% w 2017 roku (1/6 roku). Tym samym rośnie udział dni letnich (linia czerwona - ostatnio jest ich nawet 1/4 w roku). Można więc powiedzieć, że mamy coraz więcej cieplejszych dni w roku, a tym samym, że klimat się ociepla. 60 lat to nie jest co prawda długi okres jeśli chodzi o badanie zmian klimatycznych, ale lepsze to niż 10 lat (w ciągu ostatnich 10 lat zmiany nie są tak bardzo widoczne).
Temperatura w marcu
Prześledźmy te zmiany na przykładzie samego marca:
mTemps %>%
filter(Miesiac == 3) %>%
# średnia temperatura dzienna ze wszystkich lat
ggplot() +
geom_point(aes(make_date(2000, month(data), day(data)), mTemp, color = pora_roku)) +
geom_line(aes(make_date(2000, month(data), day(data)), mTemp_all), alpha = 0.8) +
facet_wrap(~Rok, ncol = 10) +
scale_x_date(date_labels = "%d", date_breaks = "7 days") +
theme(legend.position = "bottom") +
labs(title = "Średnie dobowe temperatury w marcu w porównaniu\ndo średniej z lat 1951-2017",
y = "Średnia dobowa temperatura °C", x = "", color = "Pora roku")
Kolor punktu to pora roku przypisana na podstawie temperatury, czarna cienka linia to średnia ze wszystkich lat. Początkowe lata to temperatura poniżej średniej, końcowe - powyżej. W dużym uproszczeniu. W 2013 wiosna przyszłą późno, początek lat 1990 to sytuacja odwrotna, co widzieliśmy już wyżej.
Czy klimat się ociepla?
Z reguły sprawdzenie czy klimat się ociepla liczone jest jako odchyłki temperatury od średniej na przestrzeni lat. Zbudujmy taki histogram dla każdego roku - ile dni było cieplejszych niż średnia wieloletnia dla danego dnia?
mTemps %>%
ggplot() +
geom_histogram(aes(dTemp, fill = as.factor(Rok)), binwidth = 0.5, show.legend = FALSE) +
geom_vline(xintercept = 0, color = "red") +
facet_wrap(~Rok, ncol = 10) +
labs(title = "Rozkład różnicy temperatur dziennych w stosunku\ndo średniej dla danego dnia z lat 1951-2017",
x = "Różnica temperatury °C", y = "Liczba dni")
Im więcej mamy w słupkach po prawej stronie czerwonych linii tym więcej dni cieplejszych. Wyszło nam, że mamy coraz więcej dni letnich a coraz mniej zimowych, ale czy powyższy wykres to potwierdza? Bez policzenia dokładnych statystyk opisujących poszczególne rozkłady trudno jednoznacznie powiedzieć (na oko) czy jest cieplej czy nie. Być może zajmiemy się tym w innym wpisie
Te same dane możemy narysować w inny sposób:
mTemps %>%
# co 5 lat + ostatni rok extra
filter(Rok %% 5 == 0 | Rok == 2017) %>%
mutate(Rok = fct_rev(as.factor(Rok))) %>%
ggplot() +
geom_density_ridges(aes(dTemp, Rok, fill = as.factor(Rok)), show.legend = FALSE) +
geom_vline(xintercept = 0, color = "red") +
scale_x_continuous(limits = c(-10, 7.5)) +
labs(title = "Rozkład różnicy temperatur dziennych w stosunku\ndo średniej dla danego dnia z lat 1951-2017",
x = "Różnica temperatury °C", y = "Rok")
I znowu: czy jednoznacznie można powiedzieć, że z prawej strony w poszczególnych latach mamy większą część góry?
Zagregujmy dane na poziomie miesięcy (tu to już w ogóle machinacje ze statystyką - średnia z odchyłek od średniej…) i zobaczmy, czy były miesiące zdecydowanie cieplejsze (lub zimniejsze):
mTemps %>%
group_by(Rok, Miesiac) %>%
summarise(m_dTemp = mean(dTemp)) %>%
ungroup() %>%
mutate(Miesiac = factor(Miesiac, levels = 12:1)) %>%
ggplot() +
geom_tile(aes(Rok, Miesiac, fill = m_dTemp), color = "gray10") +
scale_fill_gradient2(low = "blue", mid = "#ffffaa", high = "red", midpoint = 0) +
theme(legend.position = "bottom") +
labs(title = "Miesiące cieplejsze i zimniejsze w stosunku\ndo średniej z lat 1951-2017",
x = "Rok", y = "Miesiąc", fill = "Średnia miesięczna różnica temperatury w °C ")
Wygląda na to, że styczeń i luty w latach 1955, 1963, jakieś 1985 to były zimy stulecia (średnio chłodniej o ponad 5°C niż średnia 60-letnia). Wyraźnie cieplejszych miesięcy letnich nie widać.
Wróćmy na poziom dni i podobnie jak poprzednio dla liczby dni danej pory roku policzmy liczbę dni cieplejszych i zimniejszych niż wieloletnia średnia:
mTemps %>%
mutate(cieplej = dTemp > 0) %>%
group_by(Rok, cieplej) %>%
summarise(n = n()) %>%
ungroup() %>%
group_by(Rok) %>%
mutate(cieplej_p = 100*n/sum(n)) %>%
ungroup() %>%
ggplot() +
geom_col(aes(Rok, cieplej_p, fill = cieplej)) +
geom_hline(yintercept = 50) +
scale_fill_manual(values = c("blue", "red"), labels = c("Zimniej", "Cieplej")) +
theme(legend.position = "bottom") +
labs(title = "Podział dni cieplejszych i zimniejszych\nniż średnia z lat 1951-2017",
x = "Rok", y = "Udział %", fill = "Cieplej?")
Podział biega wokoło 50%, im bliżej współczesności tym chyba jednak więcej czerwonego (dni cieplejszych). Zobaczmy to w inny sposób:
mTemps %>%
mutate(cieplej = dTemp > 0) %>%
group_by(Rok, cieplej) %>%
summarise(n = n()) %>%
ungroup() %>%
group_by(Rok) %>%
mutate(cieplej_p = 100*n/sum(n)) %>%
ungroup() %>%
ggplot() +
geom_smooth(aes(Rok, cieplej_p, color = cieplej), method = "loess") +
scale_y_continuous(limits = c(0, 100)) +
scale_color_manual(values = c("blue", "red"), labels = c("Zimniej", "Cieplej")) +
theme(legend.position = "bottom") +
labs(title = "Podział dni cieplejszych i zimniejszych\nniż średnia z lat 1951-2017",
x = "Rok", y = "Udział %", color = "Cieplej?")
Potwierdzałoby się po raz kolejny - po 1990 roku w Polsce mamy cieplej. Czyżby piosenka Cztery pokoje (Polacy są tak agresywni, a to dlatego, że nie ma słońca nieomal przez siedem miesięcy w roku, a lato nie jest gorące, tylko zimno i pada, zimno i pada na to miejsce w środku Europy) nie była prawdziwa?
Czy jest różnica między temperaturami w kraju?
Oczywiście, że jest ale sprawdźmy to. Poszukajmy najpierw stacji, które mają najwięcej pomiarów:
# z jakich stacji mamy najwięcej pomiarów?
dbGetQuery(dbase,
"SELECT Nazwa_stacji, Kod_stacji, COUNT(*) AS n
FROM imgw
GROUP BY Kod_stacji ORDER BY n DESC LIMIT 10;") %>%
mutate(Nazwa_stacji = iconv(Nazwa_stacji, from = "latin2"))
| Nazwa stacji | Kod stacji | Liczba pomiarów |
|---|---|---|
| KRAKÓW-OBSERWATORIUM | 250190390 | 24503 |
| KÓRNIK | 252170210 | 24503 |
| WARSZAWA-BIELANY | 252200150 | 24503 |
| PUŁTUSK | 252210050 | 24503 |
| LIDZBARK WARMIŃSKI | 254200080 | 24503 |
| PUCZNIEW | 251190050 | 24473 |
| TOMASZÓW LUBELSKI | 250230070 | 24472 |
| SKIERNIEWICE | 251200030 | 24472 |
| PUŁAWY | 251210120 | 24472 |
| WIELICHOWO | 252160230 | 24472 |
24503 codziennych pomiarów to jakieś 67 lat, tak więc pierwsze pięć stacji z powyższej listy nadaje się idealnie.
stacje <- tribble(
~Nazwa_stacji, ~Kod_stacji,
"Kraków", 250190390,
"Kórnik", 252170210,
"Warszawa", 252200150,
"Pułtusk", 252210050,
"Lidzbark Warmiński", 254200080)
tab <- dbGetQuery(dbase,
"SELECT Rok, Miesiac, Dzien, Kod_stacji, MeanTemp
FROM imgw
WHERE Kod_stacji IN (250190390, 252170210, 252200150, 252210050, 254200080);") %>%
mutate(data = make_date(Rok, Miesiac, Dzien)) %>%
left_join(stacje, by = "Kod_stacji") %>%
# kolejność stacji od północy (mniej więcej)
mutate(Nazwa_stacji = fct_relevel(Nazwa_stacji, "Lidzbark Warmiński", "Pułtusk", "Warszawa", "Kórnik", "Kraków"))
Przypisaliśmy kodom stacji ręcznie nazwy (trochę je upraszczając) i pobraliśmy dane o średniej temperaturze zanotowanej w tych stacjach. Złączyliśmy dane z nazwami stacji. Teraz policzmy na którym miejscu (od 1 do 5) pod względem temperatury była każda z wybranych stacji każdego dnia:
tab %>%
group_by(data) %>%
# od najniższej temperatury
arrange(MeanTemp) %>%
mutate(poz = 1:n()) %>%
ungroup() %>%
# ile razy dana stacja była najcieplejsza itp
count(Nazwa_stacji, poz) %>%
group_by(Nazwa_stacji) %>%
mutate(p = 100*n/sum(n)) %>%
ungroup() %>%
ggplot() +
geom_tile(aes(Nazwa_stacji, poz, fill = p), show.legend = FALSE, color = "gray50") +
geom_text(aes(Nazwa_stacji, poz, label = sprintf("%.1f%%", p))) +
scale_y_reverse() +
scale_fill_distiller(palette = "RdYlBu") +
labs(title = "Bieguny ciepła i zimna według stacji meteo",
y = "Pozycja w kolejności od najzimniejszych (1)\ndo najcieplejszych (5)", x = "")
W prawie połowie dni Lidzbark Warmiński był stacją z najniższą temperaturą. Na drugiej stronie mamy Kraków, który jest najbardziej na południu (spośród wybranych stacji). Zwróćcie uwagę, że nie ma tutaj stacji w górach (Śnieżka, Kasprowy Wierch). Wyszło to trochę przypadkowo (nie ma pełnej historii pomiarów dla 67 lat), ale dobrze że tak się stało - dzięki temu porównujemy stacje położone mniej więcej na podobnym terenie.
Mapki 2017
To będzie fajna zabawa - zrobimy animowaną mapę zmiany temperatury na przestrzeni jednego roku, dzień po dniu. Wybierzemy sobie rok 2017. Pomysł powstał na bazie tekstu o regionalizmach w Niemczech, skąd przede wszystkim zaczerpnąłem informacje o bibliotece kknn. Tam mechanizm jest bardziej złożony. Ale po kolei. Potrzebujemy danych ze wszystkich stacji, ale czy są stacje które mają komplet danych? Poszukajmy tych stacji, które dokonały pomiarów w 2017 roku:
stacje2017 <- dbGetQuery(dbase,
"SELECT Kod_stacji, COUNT(*) AS n
FROM imgw
WHERE Rok = 2017
GROUP BY Kod_stacji;") %>%
# tylko te, z pełnym rokiem pomiarów
filter(n == 365) %>%
# na potrzebny JOINa:
mutate(Kod_stacji = as.character(Kod_stacji))
Mamy listę stacji, ale do mapy potrzebujemy ich położenia. To nie było łatwe (dane na stronie IMGW ich nie posiadają, a szkoda), ale udało się znaleźć stronę z mapą stacji (i aktualnymi warunkami pogodowymi). A skoro na tej stronie są zaznaczone na mapie punkty, to muszą być gdzieś ich współrzędne. Po chwili poszukiwań w kodzie znalazło się miejsce z danymi w formacie JSON. Zatem bierzemy:
# położenie stacji
stacje_lokalizacje <- fromJSON("http://monitor.pogodynka.pl/api/map/?category=meteo") %>%
select(Kod_stacji = i, long = lo, lat = la, Nazwa_stacji = n) %>%
mutate(Kod_stacji = as.character(Kod_stacji))
Teraz pozostaje wyszukanie tych stacji, które mają dane za pełny 2017 rok oraz dla których mamy położenie:
stacje2017 <- left_join(stacje2017, stacje_lokalizacje, by = "Kod_stacji") %>%
filter(!is.na(long))
# kontury Polski
mapa <- map_data("world") %>% filter(region == "Poland")
Zobaczmy w których miejscach leżą stacje?
ggplot() +
geom_polygon(data = mapa, aes(long, lat), color = "black", fill = "white") +
geom_point(data = stacje2017, aes(long, lat)) +
coord_map() +
theme_void()
Mamy sporo punktów blisko Tart, w miarę równomiernie rozłożone punkty po prawej stronie Wisły. Szkoda, że nie ma danych z okolic Szczecina - to zaburzy obliczenia.
Najpierw pobieramy z naszych danych temperaturę w 2017 roku w wybranych stacjach pomiarowych:
temperatura2017 <- dbGetQuery(dbase,
paste0("SELECT Kod_stacji, Miesiac, Dzien, MeanTemp
FROM imgw
WHERE Rok = 2017 AND Kod_stacji IN (", paste(stacje2017$Kod_stacji, collapse = ", "), ");")) %>%
mutate(Kod_stacji = as.character(Kod_stacji))%>%
left_join(stacje_lokalizacje, by = "Kod_stacji") %>%
mutate(data = make_date(2017, Miesiac, Dzien))
Teraz czas na przygotowanie danych dla kknn. Idea polega na tym, żeby dla każdego punktu na jakiejś siatce poszukać punktu najbliższego dla któego mamy dane i taką samą wartość przypisać badanemu punktowi siatki. Trochę taka ekstrapolacja. Dla danych kategorycznych (jak w przykładzie o regionalizmach) przypisalibyśmy prawdopodobieństwo przynależności do danej grupy i na jego podstawie odpowiednio nasycili kolor punktu siatki. Ale dla zmiennej ciągłej po prostu bierzemy wartość (też można by ją policzyć na podstawie na przykład kilku najbliższych punktów z jakąś średnią ważoną). Budujemy więc siatkę punktów równomiernie rozsianych po Polsce:
# siatka punktów w Polsce
poland_grid <- expand.grid(long = seq(min(mapa$long), max(mapa$long), 0.05),
lat = seq(min(mapa$lat), max(mapa$lat), 0.05))
W tym miejscu dobrze byłoby sprawdzić, które z punktów siatki leżą w Polsce i odciąć te, które są poza granicami. Do tego celu przydatna będzie funkcja gWithin() z pakietu rgeos - więcej szczegółów w którymś poprzednim wpisie. Nie będziemy liczyć temperatury dla całej siatki od razu, a dzień po dniu. A żeby przygotować animację zapiszemy serię obrazków, które później (już poza R) połączymy w jedno. Dlatego też poniższa funkcja od razu (dla wskazanego dnia) interpoluje temperaturę i zapisuje gotowa mapę.
save_temp_map <- function(f_data) {
# dane ze wszystkich stacji tego dnia
temp_dzien <- temperatura2017 %>%
filter(data == f_data)
# dane potrzebne do obliczenia wartości pozostałych punktów w siatce
temp_dzien <- temp_dzien %>%
select(long, lat, MeanTemp)
# obliczamy wartości temperatury dla wszystkich punktów siatki
temp_dzien_kknn <- kknn(MeanTemp~., temp_dzien, poland_grid, distance = 1, kernel = "gaussian")
# dodajemy obliczone temperatury do siatki punktów
temp_dzien_wynik <- poland_grid %>%
mutate(MeanTemp = fitted(temp_dzien_kknn))
# rysujemy całą mapę
plot <- ggplot(temp_dzien_wynik) +
# warstwa z temperaturami
geom_point(aes(long, lat, color = MeanTemp)) +
# warstwa z konturami Polski
geom_polygon(data = mapa, aes(long, lat, group = group), color = "black", fill = NA) +
# stala skala kolorow dla wszystkich dni
scale_color_gradient2(low = "blue", mid = "#ffffaa", high = "red", midpoint = 0,
limits = c(min(temperatura2017$MeanTemp), max(temperatura2017$MeanTemp))) +
coord_map() +
theme_void() +
theme(legend.position = "bottom") +
labs(title = f_data,
subtitle = sprintf("Średnia temperatura w kraju:%*.1f°C", 5, mean(temp_dzien$MeanTemp)),
color = "Średnia temperatura dobowa [°C]: ")
# zapisanie wykresu na dysk
ggsave(sprintf("mapki/%03d.png", yday(f_data)), plot,
# rozmiar obrazka - 1920x1080
width = 19.20, height = 10.80, dpi = 100)
}
Dla wszystkich kolejnych dni z 2017 roku generujemy mapę i zapisujemy jako plik PNG:
seq(as_date("2017-01-01"), as_date("2017-12-31"), by = "day") %>% lapply(save_temp_map)
Teraz trzeba złączyć te 365 plików w animację. Potrzebny będzie na przykład ImageMagick (opis jak zainstalować go na Ubuntu znajdziesz tutaj) i wywołane z shella:
convert -delay 50 -loop 0 *.png animation.gif
lub ffmpeg użyte:
ffmpeg -framerate 5 -i %03d.png -c:v libx264 -r 30 -pix_fmt yuv420p out.mp4
Wygenerowałem to (trochę trwa), zobaczyłem i… Pogoda zmienia się zbyt szybko, nie widać ładnych przepływów ciepłego albo zimnego powietrza przez kraj. Lepsze byłyby dane godzinowe. Nie mamy aż tak dokładnych, ale można skorzystać z danych z pomiarem trzy razy na dobę (o godzinach 6, 12 i 18). W takim przypadku odpowiednio należy dostosować loadery, pobrać dane (jest ich odpowiednio więcej) i przetworzyć całość jeszcze raz. Odpowiedni skrypt pobierający dane znajdziecie na GitHubie, skrypt robiący z nich mapę też tam jest. Wynik zaś wygląda następująco (polecam powiększyć): Niestety tylko trzy pomiary na dobę i mała liczba stacji na terenie całego kraju (braku na północno-zachodniej stronie) psuje efekt. W każdym razie idea słuszna :)
Śnieg w Wigilie
Często przy okazji spotkań przy bożonarodzeniowym stole usłyszeć można a kiedyś to dopiero były zimy, nie to co teraz!. Sprawdźmy. Zobaczmy jak wyglądała średnia temperatura w okresie Bożego Narodzenia i jaka była średnia grubość pokrywy śnieżnej.
# średnia temperatura dzienna ze wszystkich stacji, dzień po dniu
pogoda_bn <- dbGetQuery(dbase,
"SELECT Rok, Miesiac, Dzien, AVG(WysSniegu) AS WysSniegu, AVG(MeanTemp) AS MeanTemp
FROM imgw
WHERE Rok <> 2018 AND Miesiac = 12 AND Dzien IN (22, 23, 24, 25, 26, 27, 28)
GROUP BY Rok, Miesiac, Dzien;") %>%
# robimy sobie pole z datą
mutate(data = make_date(Rok, Miesiac, Dzien))
pogoda_bn %>%
group_by(Rok) %>%
summarise(snieg = mean(WysSniegu),
temp = mean(MeanTemp)) %>%
ungroup() %>%
ggplot() +
geom_col(aes(Rok, snieg), fill = "lightblue") +
geom_smooth(aes(Rok, snieg), color = "blue", se = FALSE) +
geom_point(aes(Rok, temp), color = "red", alpha = 0.7) +
geom_smooth(aes(Rok, temp), color = "red", se = FALSE) +
labs(title = "Śnieg i temperatura w okresie Bożego Narodzenia (22-28 grudnia)",
x = "Rok", y = "Grubość pokrywy śnieżnej w cm (słupki, niebieska linia)\nŚrednia temperatura dobowa w °C (punkty, czerwona linia)")
Rzeczywiście wygląda na to, że śniegu coraz mniej (ale jednak średnio w całym kraju - pamiętajcie o górach! - jakiś jest), temperatura coraz większa (od 2010 roku to wręcz powyżej zera).
No to wiemy, że kiedyś to były zimy, a co z Wielkanocą?
Czy Lany Poniedziałek zawsze był zimny?
Kiedyś to było gorąco, latało się w krótkich gaciach, z wiadrami pełnymi wody, teraz tylko pada i nikt nie ma ochoty na polewanie się wodą - tak mówią starzy górale (podobno). Zadanie podobne do poprzedniego, ale przecież Wielkanoc jest świętem ruchomym! Jest jednak sposób na wyznaczenie daty Wielkanocy. W ślad z stroną skorzystamy z metody Meeusa/Jonesa/Butchera przedstawiona przez Jeana Meeusa w jego książce Astronomical Algorithms w 1991 roku. Metoda polega na serii dość prostych obliczeń:
data_Wielkanocy <- function(rok) {
a <- rok %% 19
b <- floor(rok/100)
c <- rok %% 100
d <- floor(b/4)
e <- b %% 4
f <- floor((b+8)/25)
g <- floor((b-f+1)/3)
h <- (19*a + b - d - g +15) %% 30
i <- floor(c/4)
k <- c %% 4
l <- (32+2*e+2*i-h-k) %% 7
m <- floor((a + 11*h + 22*l)/451)
p <- (h + l - 7*m + 114) %% 31
dzien <- p + 1
miesiac <- floor((h + l - 7*m + 114)/31)
return(make_date(rok, miesiac, dzien))
}
# wyznaczamy daty Lanych Poniedziałków
lane_poniedzialki <- tibble(data = data_Wielkanocy(1951:2017) + days(1))
Mając daty Lanych Poniedziałków pobieramy dla nich pogodę - temperaturę i sumę opadów:
# pobieramy temperaturę i sumę opadów dla wszystkich dni
pogoda_wielkanoc <- dbGetQuery(dbase,
"SELECT Rok, Miesiac, Dzien, AVG(SumaOpadow) AS SumaOpadow, AVG(MeanTemp) AS MeanTemp
FROM imgw
WHERE Rok <> 2018
GROUP BY Rok, Miesiac, Dzien;") %>%
# robimy sobie pole z datą
mutate(data = make_date(Rok, Miesiac, Dzien))
# JOINem dat Lanych Poniedziałków z pogodą zostawiamy tylko pogodę Lanych Poniedziałków
pogoda_wielkanoc <- left_join(lane_poniedzialki, pogoda_wielkanoc, by = "data")
Jak to wyglądało?
pogoda_wielkanoc %>%
ggplot() +
geom_col(aes(Rok, SumaOpadow ), fill = "lightblue") +
geom_smooth(aes(Rok, SumaOpadow ), color = "blue", se = FALSE) +
geom_point(aes(Rok, MeanTemp), color = "red", alpha = 0.7) +
geom_smooth(aes(Rok, MeanTemp), color = "red", se = FALSE) +
labs(title = "Opady i temperatura w Lany Poniedziałek",
x = "Rok", y = "Suma opadów w mm (słupki, niebieska linia)\nŚrednia temperatura dobowa w °C (punkty, czerwona linia)")
Deszcz mniej więcej na tym samym poziomie, ale temperatura w ostatnich latach rzeczywiście niższa. Starzy górale mają zatem dobrą pamięć ;)
Dane pogodowe (zebrane dzień po dniu) to ciekawy przykład danych, które można potraktować jako szeregi czasowe. I zapewne pod tym kątem przyjrzymy się im w kolejnym wpisie.