Pogoda w Warszawie
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:
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:
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:
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.
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:
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.
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:
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?
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:
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:
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:
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:
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:
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:
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:
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:
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:
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.
# 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?
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:
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:
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:
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.
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).
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ą.