A lato było cieplejsze…

Czytanie zajmie Ci około 13 minut

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:

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).

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

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:

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:

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:

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?

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:

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):

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:

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:

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:

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.

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:

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:

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:

Teraz pozostaje wyszukanie tych stacji, które mają dane za pełny 2017 rok oraz dla których mamy położenie:

Zobaczmy w których miejscach leżą stacje?

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:

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:

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ę.

Dla wszystkich kolejnych dni z 2017 roku generujemy mapę i zapisujemy jako plik PNG:

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:

lub ffmpeg użyte:

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.

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ń:

Mając daty Lanych Poniedziałków pobieramy dla nich pogodę – temperaturę i sumę opadów:

Jak to wyglądało?

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.

8 komentarzy do wpisu „A lato było cieplejsze…”

  1. Ś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.

  2. 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 :(

  3. 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.

Dodaj komentarz