Rozmowy o pogodzie. Taki small talk.
Z potwierdzeniem liczbowym faktów znanych z życia codziennego. Oraz – kiedy na urlop?
Dane pobieramy ze strony wunderground.com, gdzie każdy dzień ma swój URL (przykład), dzięki czemu możemy dość łatwo przejść dzień po dniu w pętli. Tym razem będzie bez gotowego kodu, a jedynie wskazówka co jest potrzebne:
- biblioteka rvest
- po wczytaniu strony dla konkretnego dnia trzeba znaleźć tabelkę (inspektor HTMLa w Chrome bardzo przydatny do tego celu)
- z pomocą html_table() pobieramy dane z HTMLa do tabeli
- dla każdego kolejnego dnia robimy to samo, łącząc w jednej tabeli wszystkie dane
Tak przygotowane dane warto zapisać lokalnie i później korzystać już z lokalnego pliku.
Pobrałem dane dla całych lat 2000 – 2016, dla Warszawy. Po uporządkowaniu typów danych i wyczyszczeniu ich z braków lub błędów możemy zobaczyć jak wyglądają. Na przykład temperatura:
|
1 2 3 |
dane %>% ggplot() + geom_point(aes(Time, Temp), color="gray") + geom_smooth(aes(Time, Temp), size=2, color="red") |
Widać wyraźną sezonowość. Co ciekawe – średnia temperatura w Warszawie przez te 17 lat to 9.2 stopnia Celsiusza. Szczerze mówiąc sądziłem, że jest to 2-3 stopnie więcej.
To samo możemy zrobić dla innych parametrów – WindSpeed, Pressure, Humidity, Conditions. Oszczędzę wykresów, odpowiedni kawałek kodu:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 |
dane %>% ggplot() + geom_point(aes(Time, WindSpeed)) + geom_smooth(aes(Time, WindSpeed)) dane %>% ggplot() + geom_point(aes(Time, Pressure)) + geom_smooth(aes(Time, Pressure)) dane %>% ggplot() + geom_point(aes(Time, Humidity)) + geom_smooth(aes(Time, Humidity)) dane %>% ggplot() + geom_point(aes(Time, Conditions)) |
Aby powiedzieć czy pogoda się zmienia musimy mieć punkt odniesienia. Można porównywać dane okres do okresu – pogoda jest zdecydowanie sezonowa, co widać na wykresie z temperaturą, okres zmienności to rok. Można więc porównać dane rok do roku, ale co nam z tego przyjdzie, skoro warunki w różnych latach mogą się zmienić? Lepiej porównać dane do uśrednionych wartości po całym okresie. Tak też sprawdza się czy mamy do czynienia z globalnym ociepleniem – należy wziąć średnią długoterminową (na przykład z całego wieku), a potem porównać odstępstwa od tej średniej w poszczególnych latach.
Nie mamy tutaj danych z całego wieku, ale średnia z prawie 20 lat powinna coś pokazać w perspektywie zmian w ciągu dnia.
Policzmy średnie i odstępstwa od nich dla poszczególnych godzin:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
godziny <- dane %>% # rozdzielamy datę na rok - miesiąc - dzień - godzinę mutate(hour=hour(Time), year=year(Time), month=month(Time), day=day(Time)) %>% # dla każdego dnia liczymy średnią dzienną parametrów group_by(month, day) %>% mutate(mTemp = mean(Temp, na.rm = TRUE), mWindSpeed = mean(WindSpeed, na.rm = TRUE), mPressure = mean(Pressure, na.rm = TRUE), mHumidity = mean(Humidity, na.rm = TRUE)) %>% ungroup() %>% # dla każej godziny liczymy odchylenie od średniej dziennej group_by(hour) %>% mutate(dTemp = Temp - mTemp, dWindSpeed = WindSpeed - mWindSpeed, dPressure = Pressure - mPressure, dHumidity = Humidity - mHumidity) %>% ungroup() |
Dodatkowo oznaczymy pory roku i faktoryzyjemy miesiące – pierwsze dla agregacji, drugie dla ładniejszych wykresów.
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 |
godziny$season <- ifelse(godziny$month >= 3, "Spring", "Winter") godziny$season <- ifelse(godziny$month >= 6, "Summer", godziny$season) godziny$season <- ifelse(godziny$month >= 9, "Autumn", godziny$season) godziny$season <- ifelse(godziny$month == 12, "Winter", godziny$season) # nazwy pór roku jako faktory godziny$season <- factor(godziny$season, levels = c("Spring", "Summer", "Autumn", "Winter")) # nazwy miesięcy jako faktory godziny$month <- factor(godziny$month, levels = 1:12, labels = month.name) |
Zobaczmy jak wyszedł podział miesięcy na pory roku? Czy nie pomyliliśmy się gdzieś w warunkach ifelse()? I ile próbek mamy w każdym z miesięcy:
|
1 |
table(godziny$month, godziny$season) |
| Spring | Summer | Autumn | Winter | |
|---|---|---|---|---|
| January | 0 | 0 | 0 | 27492 |
| February | 0 | 0 | 0 | 26392 |
| March | 27943 | 0 | 0 | 0 |
| April | 27029 | 0 | 0 | 0 |
| May | 27808 | 0 | 0 | 0 |
| June | 0 | 27835 | 0 | 0 |
| July | 0 | 29301 | 0 | 0 |
| August | 0 | 28561 | 0 | 0 |
| September | 0 | 0 | 28283 | 0 |
| October | 0 | 0 | 29537 | 0 |
| November | 0 | 0 | 28498 | 0 |
| December | 0 | 0 | 0 | 29396 |
Dobieramy paletę – polecam serwis ColorBrewer 2.0 – miesiące z danej pory roku w jednym odcieniu.
|
1 2 3 4 |
pal_month <- c("#3182bd", "#08519c", "#74c476", "#31a354", "#006d2c", "#fd8d3c", "#e6550d", "#a63603", "#9e9ac8", "#756bb1", "#54278f", "#6baed6") pal_season <- c("#08519c", "#006d2c", "#a63603", "#54278f") |
Teraz już możemy sprawdzić jak wygląda zmiana temperatury w ciągu dnia. Najpierw globalnie:
|
1 2 3 4 |
godziny %>% ggplot() + geom_boxplot(aes(hour, dTemp, group=hour)) + geom_hline(yintercept = 0, color="red") |
Widać na pierwszy rzut oka, że w dzień jest cieplej niż w nocy, najzimniej nad ranem (5-6). Ale czy miesiąc ma jakieś znaczenie?
|
1 2 3 4 5 |
godziny %>% ggplot() + geom_smooth(aes(hour, dTemp, group=month, color=month), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_month) |
Widać że najcieplejszy moment przesuwa się – zimą (niebieskie linie) najcieplej jest około południa, latem – już po południu (około 15-16). Gdyby dodać do tego wykresu informację o wysokości Słońca nad horyzontem (albo czas od wschodu Słońca w danym dniu) to przesunięcie stałoby się widoczne od razu. Po raz kolejny mamy potwierdzenie w danych świata jaki znamy.
Widać też większą rozpiętość temperatury – w lecie różnice pomiędzy najcieplejszą a najzimniejszą godziną mogą sięgać nawet 10 stopni. Lepiej widać to przy agregacji do pór roku, zamiast miesięcy:
|
1 2 3 4 5 |
godziny %>% ggplot() + geom_smooth(aes(hour, dTemp, group=season, color=season), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_season) |
Podobne wykresy można narysować dla innych parametrów:
Wiatr:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
godziny %>% ggplot() + geom_boxplot(aes(hour, dWindSpeed, group=hour)) + geom_hline(yintercept = 0, color="red") godziny %>% ggplot() + geom_smooth(aes(hour, dWindSpeed, group=month, color=month), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_month) godziny %>% ggplot() + geom_smooth(aes(hour, dWindSpeed, group=season, color=season), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_season) |
Ciśnienie:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 |
godziny %>% ggplot() + geom_boxplot(aes(hour, dPressure, group=hour)) + geom_hline(yintercept = 0, color="red") godziny %>% ggplot() + geom_smooth(aes(hour, dPressure, group=month, color=month), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_month) godziny %>% ggplot() + geom_smooth(aes(hour, dPressure, group=season, color=season), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_season) |
Wilgotność powietrza:
|
1 2 3 4 5 6 7 8 |
godziny %>% ggplot() + geom_boxplot(aes(hour, dHumidity, group=hour)) + geom_hline(yintercept = 0, color="red") godziny %>% ggplot() + geom_smooth(aes(hour, dHumidity, group=month, color=month), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_month) |
Spójrzmy na wykres danych o wilgotności zagregowany do pór roku:
|
1 2 3 4 5 |
godziny %>% ggplot() + geom_smooth(aes(hour, dHumidity, group=season, color=season), se=F) + geom_hline(yintercept = 0, color="red") + scale_color_manual(values = pal_season) |
Znowu widać większą rozpiętość w lecie, najmniejszą w zimie. Widać też, że latem wilgotność po południu jest o około 25 punktów procentowych mniejsza niż nad ranem. Słoneczko przyświeciło i wysuszyło poranną rosę. Potwierdzenie faktów.
Teraz zobaczmy na dane z pola Conditions. Na poniższym wykresie mamy obraz danych zagregowanych (z 43 do 9 cech, przy okazji agregacji przetłumaczyłem opisy na polski). Wartości na wykresie to udział procentowy danych warunków w dniu. Można więc nazwać to prawdopodobieństwem wystąpienia określonego warunku w ciągu dnia.
Co widzimy?
- burze występują głównie latem, deszcz mniej więcej po równo (czasem pada, czasem nie pada)
- zachmurzenie w Warszawie jest stosunkowo jednakowe, najbardziej słonecznym miesiącami są maj i sierpień
- z tych dwóch punktów wynika, że urlopy planowane na sierpień dają największe prawdopodobieństwo “dobrej pogody” (przejrzyste niebo w 55%, chmury w 34%, niecałe 6% na deszcz, nieco ponad 1% na burzę, 3% mgły ale raczej w nocy i nad ranem)
- na grad możemy trafić na wiosnę i jesienią – prawdopodobieństwo jednak jest niewielkie (0.03%)
- w ciepłych miesiącach szansę na mgłę są trzykrotnie mniejsze niż np. w listopadzie
- śnieg – co dość oczywiste – pada w zimie, najwięcej w styczniu
Widać to wszystko również na wykresie z warunkami zagregowanymi do miesięcy:
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
dane %>% mutate(month=month(Time)) %>% count(month, Conditions) %>% ungroup() %>% group_by(month) %>% mutate(p=100*n/sum(n)) %>% ungroup() %>% left_join(CondDict_df, by=c("Conditions"="eng")) %>% filter(pol != c("b/d"), !is.na(pol)) %>% ggplot() + geom_bar(aes(month, p, fill=pol), stat="identity", show.legend = F) + facet_grid(pol~., scales = "free_y") + scale_x_continuous(breaks = 1:12, labels = month.name) + labs(x="Miesiąc", y="%") + theme(strip.text.y = element_text(angle=0), axis.text.x = element_text(angle=30, hjust=1)) |
Samodzielnie proszę zrobić to samo w podziale na poszczególne godziny. Pada bardziej po południu czy nad ranem? Podpowiedź – agregujemy po godzinach.
Wyszło mi, że burze są popołudniami, grad pada wciągu dnia, a mgły występują raczej nocą.
Sprawdziliśmy jak wyglądają odchylenia od średnich dziennych poszczególnych czynników w ciągu dnia, ale jak wyglądają wartości uśrednione po wszystkich latach?
Najpierw policzmy średnią dzienną (ze wszystkich godzin w ciągu dnia) dla poszczególnych czynników:
|
1 2 3 4 5 6 7 8 |
dane_day_means <- dane %>% mutate(year=year(Time), month=month(Time), day=day(Time)) %>% group_by(year, month, day) %>% summarise(mTemp = mean(Temp, na.rm = TRUE), mWindSpeed = mean(WindSpeed, na.rm = TRUE), mPressure = mean(Pressure, na.rm = TRUE), mHumidity = mean(Humidity, na.rm = TRUE)) %>% ungroup() |
Teraz średnią z całego miesiąca:
|
1 2 3 4 5 6 7 8 |
dane_month_means <- dane %>% mutate(year=year(Time), month=month(Time)) %>% group_by(year, month) %>% summarise(mTemp_m = mean(Temp, na.rm = TRUE), mWindSpeed_m = mean(WindSpeed, na.rm = TRUE), mPressure_m = mean(Pressure, na.rm = TRUE), mHumidity_m = mean(Humidity, na.rm = TRUE)) %>% ungroup() |
A teraz zobaczmy te dane – temperatura:
|
1 2 3 4 5 6 7 8 |
ggplot() + geom_point(data=dane_day_means, aes(make_date(year, month, day), mTemp), color="gray") + geom_line(data=dane_month_means, aes(make_date(year, month, 15), mTemp_m), color="red", size=1) + labs(x="Data") |
Widać po raz kolejny sezonowość. Widać, że zima 2006 roku była najzimniejsza, a lato 2007 najcieplejsze jeśli chodzi o średnią dzienną, chociaż rekordy ciepła padały w 2015 roku. Może ktoś pamięta czy tak było? Ja szczerze mówiąc nie pamiętam i nigdy do tego nie przywiązywałem wagi – nie pamiętam czy w marcu zeszłego roku padał śnieg (w razie potrzeby można sięgnąć do danych).
Podobnie możemy zrobić dla innych czynników, chociaż poza wilgotnością dane już nie są tak wyraźnie sezonowe.
|
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 |
# predkosc wiatru ggplot() + geom_point(data=dane_day_means, aes(make_date(year, month, day), mWindSpeed), color="gray") + geom_line(data=dane_month_means, aes(make_date(year, month, 15), mWindSpeed_m), color="red", size=1) + labs(x="Data") # cisnienie ggplot() + geom_point(data=dane_day_means, aes(make_date(year, month, day), mPressure), color="gray") + geom_line(data=dane_month_means, aes(make_date(year, month, 15), mPressure_m), color="red", size=1) + labs(x="Data") # wilgotnosc ggplot() + geom_point(data=dane_day_means, aes(make_date(year, month, day), mHumidity), color="gray") + geom_line(data=dane_month_means, aes(make_date(year, month, 15), mHumidity_m), color="red", size=1) + labs(x="Data") |
A jak wygląda temperatura (i inne czynniki) dla danego miesiąca, po uśrednieniu wartości ze wszystkich lat?
|
1 2 3 4 5 6 7 8 |
dane_year_means <- dane %>% mutate(month=month(Time)) %>% group_by(month) %>% summarise(mTemp_y = mean(Temp, na.rm = TRUE), mWindSpeed_y = mean(WindSpeed, na.rm = TRUE), mPressure_y = mean(Pressure, na.rm = TRUE), mHumidity_y = mean(Humidity, na.rm = TRUE)) %>% ungroup() |
Zobaczmy:
|
1 2 3 4 |
ggplot(dane_year_means) + geom_line(aes(month, mTemp_y), color="gray") + geom_point(aes(month, mTemp_y), color="black", size=2) + scale_x_continuous(breaks = 1:12, labels = month.name) |
Średnio rzecz biorąc – lato niezbyt gorące jak to twierdzi Staszewski w “4 pokojach”. Faktycznie – średnio 20 stopni Celsiusza w lipcu to nie jest szał. Zim nie mamy też super mroźnych (jakieś -3 w styczniu).
Ciekawe są inne parametry – prędkość wiatru – najsłabiej wieje w lecie:
|
1 2 3 4 |
ggplot(dane_year_means) + geom_line(aes(month, mWindSpeed_y), color="gray") + geom_point(aes(month, mWindSpeed_y), color="black", size=2) + scale_x_continuous(breaks = 1:12, labels = month.name) |
Ciśnienie:
|
1 2 3 4 |
ggplot(dane_year_means) + geom_line(aes(month, mPressure_y), color="gray") + geom_point(aes(month, mPressure_y), color="black", size=2) + scale_x_continuous(breaks = 1:12, labels = month.name) |
Ten listopad wygląda podejrzanie – należałoby mu się przyjrzeć, czy nie ma tam jakichś błędnych danych. Wykres pokazuje jakieś wielkie rozpiętości, ale proszę zwróćcie uwagę na skalę osi Y – różnice są na poziomie pojedynczych hektopaskali. Żadne (0.1%) właściwie, tak samo jak dla wyżej pokazanego wiatru – sądzę, że nikt nie czuje różnicy pomiędzy wiatrem wiejącym z prędkością 12 a 13 km/h.
Przy wilgotności powietrza jest już nieco inaczej – rozpiętość jest na poziomie nawet 20 punktów procentowych – w lecie (Słoneczko) suszy bardziej.
|
1 2 3 4 |
ggplot(dane_year_means) + geom_line(aes(month, mHumidity_y), color="gray") + geom_point(aes(month, mHumidity_y), color="black", size=2) + scale_x_continuous(breaks = 1:12, labels = month.name) |
Na koniec zobaczmy czy widać efekt cieplarniany? Porównajmy średnie miesięczne temperatury (czyli: średnia z dni danego miesiąca, z uwzględnieniem roku) do średniej temperatury w miesiącu z 17 lat (czyli na przykład średnia dla wszystkich dni kwietnia, ze wszystkich kwietni).
|
1 2 3 4 5 6 7 8 9 10 |
dane_deltas <- left_join(dane_month_means, dane_year_means, by="month") %>% mutate(dTemp = mTemp_m - mTemp_y, dWindSpeed = mWindSpeed_m - mWindSpeed_y, dPressure = mPressure_m - mPressure_y, dHumidity = mHumidity_m - mHumidity_y) ggplot(dane_deltas) + geom_hline(yintercept = 0, color="red") + geom_point(aes(year, dTemp)) + geom_smooth(aes(year, dTemp)) |
Na upartego można powiedzieć, że w ostatnich trzech latach widać jakiś wzrost. Średnią długookresową powinniśmy policzyć dla kilkudziesięciu lat (np. wszystkie kwietnie XX wieku) i porównać z średnią miesięczną z danego roku (średnia z kwietnia 2016). Jeśli chcesz sprawdzić czy efekt globalnego ocieplenia jest widoczny w danych – poszukaj danych i pobaw się nimi trochę. Instrukcję właśnie masz za sobą.
nie skumałem jednego. rozumiem rozpiętość temperatury latem jako największą. ale dlaczego na wykresach spada do do poniżej zero?
Na którym wykresie? Może po prostu były jakieś nocne przymrozki?