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:
|
1 2 3 4 5 6 7 8 9 10 |
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).
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
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
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
# 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:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
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:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
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:
|
1 2 3 4 5 6 7 8 9 10 11 |
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?
|
1 2 3 4 5 6 7 |
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:
|
1 2 3 4 5 6 7 8 9 10 |
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):
|
1 2 3 4 5 6 7 8 9 10 11 |
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:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
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:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
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:
|
1 2 3 4 5 6 |
# 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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
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:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
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:
|
1 2 3 4 5 6 7 8 9 |
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:
|
1 2 3 4 |
# 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:
|
1 2 3 4 5 |
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?
|
1 2 3 4 5 |
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:
|
1 2 3 4 5 6 7 |
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:
|
1 2 3 |
# 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ę.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 |
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:
|
1 |
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:
|
1 |
convert -delay 50 -loop 0 *.png animation.gif |
lub ffmpeg użyte:
|
1 |
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ć):
Ś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.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
# ś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ń:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 |
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:
|
1 2 3 4 5 6 7 8 9 10 11 |
# 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?
|
1 2 3 4 5 6 7 8 |
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.
Świetna robota, świetna.
Okazuje się, że dane potwierdzają powszechne odczucia, jednak wygląda również na to, że występują jakieś ~30-letnie cykle, około roku 1985 aura była z grubsza podobna jak dziś a amplituda zmian się powtarza.
Fajny wpis Łukasz :)
Niestety ale potwierdziłeś dwa mity które mnie interesowały: zimny lany poniedziałek i ciepłe Boże Narodzenie.
Pamiętam z dzieciństwa jak biegałem po ulicy z woreczkami napełnionymi wodą. W ostatni lany poniedziałek założyłem zimowy płaszcz :(
Błąd w poleceniu 'rsqlite_send_query(conn@ptr, statement)’:no such table: imgw
A jakie tabele masz w bazie?
dzięki Łukaszu. już opanowane.
Dzięki za artykuł. Dowiedziałem się dzięki Tobie, że historyczne dane meteorologiczne dostępne są jako informacja publiczna, bez konieczności web scrapingu.
Mała wskazówka: zamiast zapytania SQL można napisać:
dbase %>%
tbl(„imgw”) %>%
filter(Rok != 2018) %>%
group_by(Rok, Miesiac, Dzien) %>%
summarise(mTemp = mean(MeanTemp))
Jak pisze Hadley, „surprisingly, this sequence of operations never touches the database” (https://db.rstudio.com/dplyr/). Zamiast tego, dbplyr generuje zapytanie SQL (można je podejrzeć przez show_query()), które jest wysyłane do bazy dopiero gdy użyjesz collect(). Dlatego efekt jest dokładnie taki sam: obróbka danych następuje w bazie, nie zaś w R. Dostajesz jednocześnie wydajność bazy i wygodę tidyverse. Taki R-owy ORM :)
Wszystko się zgadza i właściwie zawsze piszę w ten sposób. Tutaj chodziło mi o pokazanie, że można inaczej – brawo za wyłapanie!
Ale nie zawsze to działa, zależy od bazy danych. Grupowanie i wybieranie kolumn działa zazwyczaj bez problemu, ale już filtrowanie niekoniecznie – wszystko zależy od funkcji jakich używamy i czy da się je przełożyć na SQLa oraz czy baza te funkcje obsługuje. To potrafi być upierdliwe.
Szczerze mówiąc, ja też ostatnio coraz częściej kończę na gołym SQL. Przypomina mi się dyskusja, która kilka lat temu się przetoczyła przez sieć – czy ORM ma rację bytu (głośne swego czasu artykuły: https://maetl.net/talks/rise-and-fall-of-orm, http://solnic.eu/2015/09/18/ditch-your-orm.html). Ciekawe, czy taka dyskusja się rozwinie w tidyverse.