Efekt cieplarniany – czy widać go w danych?

Efekt cieplarniany to całkiem proste zjawisko, polegające na systematycznym ogrzewaniu objętości. Czyli, że jest jakaś objętość — np. planeta Ziemia, a my systematycznie rozpalamy na niej coraz więcej „ognisk” które ją ogrzewają. Nie byłoby to nic szczególnego, gdybyśmy mieli pod ręką drugą chłodną planetę, na którą możemy oddać nadmiar naszego ciepła. Tak jak np. pojemnik z lodem mógłby przyjąć ciepło z pomieszczenia, ogrzewając się. Kosmiczna próżnia jest jednak całkiem niezłym izolatorem termicznym. Skoro nic w niej nie ma, to nie jest to ani chłodne, ani gorące, czyli nie może zbytnio przekazywać ciepła.

Mamy jeszcze opcję wypromieniowania energii — rozgrzany do czerwoności pręt powoduje, że „czujemy” ciepło, mimo iż nie zalewa nas gorące powietrze. Wyemitowanie więc znacznych ilości takiego promieniowania w kosmos mogłoby nas trochę ochłodzić. Z tego, co kojarzę jednak, mamy takie sprytne coś jak górne warstwy atmosfery, które mają nas właśnie chronić przed m.in. takim promieniowaniem. Zdaje się więc, że całkiem sporo tego promieniowania jest odbijane z powrotem. Istnieje hipoteza (globalne ocieplenie), która mówi, że poprzez coraz większą aktywności przemysłową doprowadzamy do sytuacji, że wytwarzamy więcej ciepła, niż możemy zutylizować.

Efekt cieplarniany — dane

Ja zupełnie się nie znam na tym, jak działa efekt cieplarniany i zdecydowanie polecam zgłębienie tematu na własną rękę. Wiem jednak, że morał z tej sytuacji jest prosty — jeśli efektu nie ma, to mamy mniej więcej taką samą temperaturę. Jeśli efekt istnieje, to jest coraz cieplej. Czy możemy spróbować samemu ocenić który scenariusz jest bardziej prawdopodobny? Myślę, że możemy, szczególnie że mamy do dyspozycji dane meteorologiczne zebrane i opracowane przez Instytut Meteorologii i Gospodarki Wodnej.

Danych tych po raz pierwszy miałem użyć jako danych niezależnych przy budowie jakości modelu powietrza (aktualny najnowszy artykuł znajdziesz tutaj), jednakże temat mi się nieco przeciągnął. Zamiast więc czekać na św. Nigdy, zabierzmy się za nie tu i teraz!

Pomysł

Dane z IMGW są dostępne przez API, jak i w paczkach zip. Standardem jest, że przez API wystawiane są dane aktualne, które są wykorzystywane na bieżąco, a w paczkach zip mamy dane historyczne. Nie inaczej jest też tutaj.

Mamy więc dane meteorologiczne, hydrologiczne i aktynometryczne. Skorzystamy z tych pierwszych, bo one nam bezpośrednio pokażą, jak jest gorąco. W tym katalogu mamy podkatalogi z różnymi poziomami agregacji — miesięcznym, dobowym i terminowym. Z tego, co zrozumiałem, to terminowa agregacja ma największą rozdzielczość. Przyda się więc komuś, kto chce się mocniej wgryźć w temat. Ja natomiast zadowolę się miesięczną agregacją — każda stacja pomiarowa będzie mieć uśrednione wartości dla każdego miesiąca.

W ten sposób będziemy mogli wyliczyć sobie komentarz na: „Panie, kiedyś to w marcu było śniegu po kolana, a w lipcu to samochody przez miasto nie jeździły, bo asfalt się topił!”. Kawałek kodu Pythona i już będziemy wiedzieć, o co chodzi. A przy okazji może uda nam się zweryfikować hipotezę, że właśnie stosujemy w praktyce efekt cieplarniany.

Pobieranie danych

Zaczynamy od pobrania danych. Na szczęście IMGW udostępniło dane w postaci bardzo ładnie skatalogowanych plików. Możemy je pobrać np. za pomocą programu wget z poziomu basha:

%%bash
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1951_1955/1951_1955_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1956_1960/1956_1960_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1961_1965/1961_1965_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1966_1970/1966_1970_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1971_1975/1971_1975_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1976_1980/1976_1980_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1981_1985/1981_1985_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1986_1990/1986_1990_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1991_1995/1991_1995_m_k.zip
wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/1996_2000/1996_2000_m_k.zip

for i in {2001..2019}
do
   wget -c https://dane.imgw.pl/data/dane_pomiarowo_obserwacyjne/dane_meteorologiczne/miesieczne/klimat/$i/${i}_m_k.zip
done

Gdy już mamy te dane w plikach zip, to najsensowniej byłoby je rozpakować do jednego katalogu. Tutaj znowu użyjemy basha:

!unzip -o '*.zip' -d dane/ 

W tym momencie masz już na dysku odpowiednie pliki miesięczne z temperaturą.

Wczytanie i wybór stacji pomiarowej

Jak to zwykle z danymi bywa, można je przetwarzać na nieskończenie wiele sposobów. Ja postanowiłem przejść po wszystkich plikach i z racji tego, że mają taką samą strukturę, skleić je w jedną ramkę danych:

names = ["kod_stacji", "nazwa_stacji", "rok", "miesiac", "absolutna_temperatura_maksymalna_C", "status_pomiaru_TMAX", 
         "srednia_temperatura_maksymalna_C", "status_pomiaru_TMXS", "absolutna_temperatura_minimalna_C", 
         "status_pomiaru_TMIN", "srednia_temperatura_minimalna_C", "status_pomiaru_TMNS", 
         "srednia_temperatura_miesięczna_C", "status_pomiaru_STM", "minimalna_temperatura_przy_gruncie_C", 
         "status_pomiaru_TMNG", "miesieczna_suma_opadow_mm", "status_pomiaru_SUMM", "maksymalna_dobowa_suma_opadow_mm", 
         "status_pomiaru_OPMX", "pierwszy_dzien_wystspienia_opadu_maksymalnego", 
         "ostatni_dzien_wystapienia_opadu_maksymalnego", "maksymalna_wysokosc_pokrywy_snieznej_cm", 
         "status_pomiaru_PKSN", "liczba_dni_z_pokrywa_sniezna", "liczba_dni_z_opadem_deszczu", 
         "liczba_dni_z_opadem_sniegu"]

full_df = pd.DataFrame()

years_list = ["1951_1955", "1956_1960", "1961_1965", "1966_1970", 
              "1971_1975", "1976_1980", "1981_1985", "1986_1990", 
              "1991_1995", "1996_2000"] + list(range(2001, 2020))

for year in years_list:
    tmp_df = pd.read_csv("dane/k_m_d_{}.csv".format(year), encoding='cp1250', names = names)
    full_df = full_df.append(tmp_df, sort=False)

Okazuje się, że nie wszystkie stacje pomiarowe mają tyle samo pomiarów. Niektóre pewnie zaczęły działać wcześniej, niektóre już pewnie są rozebrane, a niektóre pewnie miały dłuższe awarie. Postanowiłem więc „na ślepo” wybrać którąś z puli tych, co mają najwięcej pomiarów:

most_observations = full_df["kod_stacji"].value_counts().head(n=1).index[0]
target = full_df[full_df["kod_stacji"] == most_observations]

Aktualnie jest to stacja WARSZAWA-BIELANY.

Wizualizacje

Jedną z ciekawszych wizualizacji jest wizualizacja, która polega na narysowaniu linii obrazującej średnią temperaturę miesięczną w ciągu roku. Jeśli narysujemy taką linię dla każdego roku, to otrzymamy łuk. A jeśli dodatkowo dla każdego roku dobierzemy inny kolor, to od razu będziemy mogli określić czy coś się zmienia. Jeśli łuk będzie „gruby”, to zmiany są duże, a jeśli kolory będą „uporządkowane” to zmiana ta jest monotoniczna, czyli po prostu rosnąca lub malejąca. Zobaczmy, jak jest u nas:

ax = sns.lineplot(x="miesiac", y="srednia_temperatura_miesięczna_C", hue="rok", data=target)

 

Jak się okazuje, miesiące zimowe są nieźle wymieszane. Coś tam się dzieje niedobrego, jednakże trudno jednoznacznie powiedzieć, o co chodzi. Tzn. raz był wtedy pewnie metry śniegu — średnia -11 C, a raz ludzie chodzili bez kurtek — średnia +5 C.

Miesiące letnie z kolei mają zdecydowanie mniejszy rozrzut, bo mieszczą się od +15 C do +24 C. Tutaj jednak kolory zdają się być uporządkowane. Co zdaje się pokazywać, że z roku na rok jest coraz cieplej.

Zobaczmy dokładniej, w jaki sposób zmieniał się wybrany miesiąc w tej stacji pomiarowej na przestrzeni lat:

target_month = full_df[(full_df["kod_stacji"] == most_observations) & (full_df["miesiac"] == 6)]
ax = sns.lineplot(x="rok", y="srednia_temperatura_miesięczna_C", data=target_month)

Widzimy, że do lat 90 średnia temperatura miesiąca się wahała. A później systematycznie zaczęła rosnąć. Być może wynikało to ze zmiany ustrojowej i innego podejścia do badań naukowych? Może zmienił się sprzęt pomiarowy na trwalszy i dokładny? A może po prostu tak wygląda to zjawisko?

Model linowy

Na powyższym wykresie mamy jedną zmienną niezależną (rok) i jedną zmienną zależną (średnia temperatura w miesiącu). Nie da się więc tutaj doszukać jakiejś kosmicznej tajemnicy. Jeśli jednak mamy przesłanki, że nasze zjawisko jest liniowe, to możemy pokusić się o wyznaczenie parametru tego zjawiska metodą najmniejszych kwadratów:

X = target_month["rok"]
y = target_month["srednia_temperatura_miesięczna_C"]

model = sm.OLS(y, X).fit()
predictions = model.predict(X) 

model.summary()

Dostajemy tutaj parametr dla zmiennej rok: 0,0089. Oznacza to, że według tego modelu, bierzemy sobie rok i mnożymy przez tę wartość. Dla roku 2019 otrzymamy: 2019 * 0,0089 = 17,9691. A dla roku 1959 byłoby to: 1959 * 0,0089 = 17,4351. Czyli podbicie modelowej średniej temperatury o 0,534 C w ciągu 60 lat. Nie wygląda to aż tak źle na wykresie:

Model liniowy ma jednak spore ograniczenia. Przede wszystkim jest liniowy. Czyli jest linią prostą. A mało co w naturze jest linią prostą. Oznaczałby on też, że w czasach narodzenia Jezusa, na terenie Warszawy średnia temperatura wynosiła 0 C. Byłby to więc chyba koniec epoki lodowcowej. A zdaje mi się, że epoka ta skończyła się o wiele wcześniej w tym miejscu. Tak więc należy na niego spoglądać z dużym przymrużeniem oka.

Konkluzja

Wiesz już jak zabrać się za dane z IMGW. I jak samemu spróbować zbadać efekt cieplarniany w Polsce. Oczywiście, zaproponowany przeze mnie przykład pozostawia wiele do życzenia:

• Dane temperaturowe mierzone są w miastach. A miasta są zawsze nieco cieplejsze niż najbliższa ich okolica. Mamy więc tutaj efekt miejskiej wyspy ciepła. Warto byłoby się rozejrzeć za bardziej „obiektywnymi” stacjami pomiarowymi.
• Mamy mało danych pomiarowych. W IMGW mamy je od 1951. Zdecydowanie byliśmy już wtedy w erze przemysłowej. Fajnie byłoby mieć takie dane z ostatnich czterech tysięcy lat. Jest to niestety niemożliwe ;-(.
• Model, który zbudowałem, dotyczył jednej stacji i jednego miesiąca. Trzeba byłoby coś takiego zrobić dla wszystkich stacji i miesięcy albo jakoś sprytnie to zrobić całościowo.
• Model liniowy jest naiwny. Warto byłoby znaleźć coś sensowniejszego.

Tak czy owak, temat jest jak najbardziej rozwojowy i do obserwowania.

Jeśli interesuje Cię jakiś temat – nie musi być związany z tym artykułem – to zostaw mi sygnał tutaj. Dzięki!

Pełny kod przykładu użytego w artykule znajduje się tutaj.