Введение

Недавно открыл для себя платформу Streamlit и был впечатлен простотой интеграции в питоновский проект. По детски, очень радовался тому что контроллеры на дашборде напрямую меняют питоновские переменные. И вот для тестирования решил поиграть с одной из тем которая мне очень интересна – климат. Начал с самого простого параметра который можно проанализировать – температуру воздуха с метеостанции города в котором я живу - Алматы (Казахстан). Интересно было посмотреть эффект глобального изменения климата на отдельно взятый город.

Для анализа скачал данные за 100 лет (1922-2022). Для начала использовал только месяцы с ноября по апрель включительно, во-первых, потому что в несвободное от работы время я работаю с морским льдом и зимний климат мне понятней, а во-вторых потому что в данных была пропущена неделя наблюдений в августе, в одном из годов.

Один из плюсов Streamlit, то что дашборды можно бесплатно публиковать в Streamlit Cloud, что я и сделал:

https://share.streamlit.io/yevkad/alaclim/main/app.py

Пройдя по этой ссылке, можно поиграть с интерактивными графиками. А в этом посте я поподробней покажу что и как я делал (с кодом) и опишу некоторые свои наблюдения.

Инструменты

Все было сделано на Python и следующими библиотеками:

Streamlit - платформа для быстрого создания DS/ML веб-приложений
Pandas - манипуляции с таблицами
Numpy - использовал некоторые вспомогательные функции
Plotly - Все интерактивные графики были построены с помощью этой библиотеки
SciPy - для подбора распределения вероятностей
pyMannKendall - для Теста Манна-Кендалла, который оценивает статистическую значимость тренда

import streamlit as st
import pandas as pd
import numpy as np
import plotly.express as px
import plotly.figure_factory as ff
import plotly.graph_objects as go
import scipy
import pymannkendall as mk

from datetime import datetime

Данные

Основная часть данных была взята с NOAA Climate Data Center, к слову - один из любимых источников открытых данных, наравне с европейским Copernicus.

Почему-то было много пропусков после 2006го года, заполнил их данными rp5 - там трехчасовые данные с метеостанций доступны для большинства локаций с 2005го года. На rp5 данные в местном времени, предварительно перевел время в UTC и усреднил подневно.

Данные подготавливаются функцией представленной ниже. Декоратор st.cahce для того чтобы данные не обрабатывались каждый раз при обновлении контроллеров.

@st.cache
def get_winter_data(fl):
		#initial columns: date, at - daily mean air temperature
    df_winter=pd.read_csv(fl,parse_dates=['date'])
    df_winter['fd']=np.abs(np.minimum(0,df_winter['at'])) #freezing days
    df_winter['month']=df_winter['date'].dt.strftime('%b')
    
    #relative date for timeseries.
    df_winter['date_rel']=pd.to_datetime(df_winter['date'].dt.strftime('2016-%m-%d'),
                                         format='%Y-%m-%d')
    df_winter.loc[df_winter['mon']>8,
                  'date_rel']=df_winter['date_rel'] - pd.DateOffset(years=1)

    df_winter.loc[df_winter['date']<datetime(1972,6,1), 'period']='1922-1972'
    df_winter.loc[df_winter['date']>datetime(1972,6,1), 'period']='1972-2022'
    df_winter['dt']=df_winter['date'].dt.strftime('%Y%m').astype(int)
    
    df_fdd=df_winter.groupby(['seas'])['fd'].sum().reset_index(name='fdd')
    
    return df_winter,df_fdd

В результате получаем два датафрейма df_winter и df_fdd. У df_winter следующие колонки:

date - дата
at - среднедневная температура воздуха
seas - сезон, строка в виде начального и конечного года (1922-1923)
month - месяц в формате строчном формате (Nov, Dec, etc). От нее можно было бы избавиться и оптимизировать обработку, но руки пока не дошли
mon - месяц в числовом формате (11, 12 и тд.)
period - фиксированный период по 50 лет с 1922 по 1972 и 1972 по 2022
fd - freezing days, градусы ниже нуля в абсолюте
dt - полукостыль, год-месяц в виде целого числа в формате YYYYMM, для группирования нескольких сезонов в одну группу функцией np.digitize (пример будет в следующих секциях)
rel_date - дата с заменой года для всех дней на 2015 до января и 2016 с января, для того что бы можно было агрегировать 100 лет за каждый день.

У df_fdd всего 100 рядов, по ряду за каждый сезон и следующие колонки:

seas - сезон
fdd - Freezing Degree Days - аккумулированные градусы ниже нуля в абсолюте. Более подробное описание будет ниже.

Временной ряд

Временной ряд был составлен выведением на график агрегированных за каждый день с ноября по апрель столетних температур (максимум, минимум, среднее) и возможностью сравнивать эти параметры с выбранными периодами.

Если визуально сравнить первую половину периода со второй, то видно, что в основном средние температуры с 1922 по 1972 ниже столетнего среднего значения, а с 1972 по 2022 выше.

Сравнение средних температур за первые и последние 50 лет со столетними средними.

Распределения

Далее были построены гистограммы с подбором распределения вероятностей для двух периодов по 50 лет: 1922-1972 и 1972-2022. В падающем меню можно выбрать тип распределения - Нормальное или Гамма. В данном случае визуально понятно, что Нормальное больше подходит, но в будущем добавлю метрику "goodness of fit".

Левая часть распределения 1972-2022 "сужается" в право, по сравнению с периодом 1922-1972. Это показывает, что за последние 50 лет количество наблюдений экстремально холодных температур уменьшилось. Это не противоречит выведенными параметрами распределения - среднее (Mean) увеличилось на 0.6 градуса, а среднеквадратическое отклонение (SD) уменьшилось по сравнению с начальным периодом.

Распределение температур за два периода: 1922-1972 и 1972-2022

Усатые Ящики

Помимо сравнения двух фиксированных периодов в виде гистограмм и распределений описанных в предыдущей секции, также была добавлена возможность сравнения распределения выбранных групп периодов (1 год, 5 лет, 10 лет и тд). Наиболее подходящей визуализацией для этого, по моему мнению, являются графики, которые я называл всегда боксплотами, и только сегодня выяснил что русское называние для них - "ящик с усами".

Для тех кто не сталкивался с такими графиками, подробное описание "ящика с усами" можно почитать, например, здесь, но грубо и приближенно, вот как можно его интерпретировать:

"Ящик" содержит 50% наблюдений
Линия внутри ящика - медиана, или 2-й квартиль - выше и ниже по 50% наблюдений
Нижний и верхний "ус" - 25% нижних и верхних наблюдений соответственно
Точки - выбросы (наблюдения вне 2.7 среднеквадратических отклонений)
Нижняя граница "ящика" - 1-й квартиль - ниже него 25% наблюдений, выше 75%
Верхняя граница "ящика" - 3-й квартиль - выше него 25% наблюдений, ниже 75%

Сгруппировав по длительным периодам, например 20 лет, как на примере ниже, видно что, медиана и 1-й квартиль стабильно возрастает с периода который включает в себя 1970е, а 3-й квартиль в последнем периоде. Также как на распределениях описанных выше, изменение зим по большей части выражается в уменьшении "холодов", а в более поздние периоды еще и в увеличении "тепла".

Бокс плот температур воздуха, сгрупированные по 20 лет.

Вот часть кода, где группируются температуры по годам и выводятся в график:

# Years grouping:
nmax=df_fdd.shape[0]
yr_start=df_winter['year'].min()

# list of only divisible groups
group_lst=[i for i in range(1,nmax+1,1) if nmax % i==0]

step=st.selectbox('Group years ', group_lst, index=group_lst.index(10))
date_range=np.array([datetime(yr_start+i,6,1).strftime('%Y%m')
for i in range(step,nmax+1,step)
]).astype(int)
date_bin_n=np.array([f'{yr_start+i}-{yr_start+step+i}' for i in range(0,nmax,step)])
df_winter_g=df_winter.copy()
#assigning bin name to date range
df_winter_g['date_bin']=date_bin_n[np.digitize(df_winter_g['dt'], date_range)]
#month filtering
mon_lst=df_winter['month'].unique()
st_mon_lst = st.multiselect("Months Used", mon_lst, default=mon_lst)
df_mon=df_winter_g[df_winter_g['month'].isin(st_mon_lst)]

fig_box = px.box(df_mon, x="date_bin", y="at",
labels={'seas':'Winter Season',
'at':'Air Temperature (\u00B0C)',
'date_bin':'Date Range'},)
st.write(fig_box)

Суровость зим по FDD

Сумма градусодней мороза (Freezing Degree Days - FDD) рассчитывается как абсолютная сумма всех отрицательных температур воздуха за период – обычно за зимний сезон. На всякий случай формула FDD:

$FDD=\sum_i^N|min(0,at_i)|$

Где at - среднедневная температура воздуха, N - число дней в сезоне.

FDD часто используется для примерного определения прироста толщины льда, уровня промерзания почвы, и может служить индексом суровости зим: чем меньше FDD тем мягче зима.

Отобразив на графике FDD за каждый сезон и положив линию тренда, то можно увидеть тенденцию снижения суровости зим, особенно с 1970х годов.

Тренд был посчитан в NumPy как полиномиальная кривая (с возможностью управлять порядком).

poly_deg=st.slider('Degree of a Polynomial Trend',
                    min_value=1, max_value=20, value=1, step=1)

trend_x = np.arange(df_fdd.shape[0])
trend_fit = np.polyfit(trend_x, df_fdd['fdd'], poly_deg)
trend_fit_fn = np.poly1d(trend_fit)

Помимо визуальной оценки, тренд и его статистическая значимость была оценена тестом Манна-Кенадалла при помощи библиотеки pyMannKendall вот так:

mk_res=mk.original_test(df_fdd['fdd']) # results of Mann-Kendall Test

И вывести результат..

trend_mk=mk_res.trend # Trend: Increasing, Decreasing or No Trend
pval_mk=mk_res.p      # P-Value of MK Test

if pval_mk<0.05: # alpha=0.05
    p_str='**statistically significant**'
    alpha_compr='smaller'
else:
    p_str='**not statistically significant**'
    alpha_compr='greater'

#Markdown string generated to present MK Test results of FDD Trend
mk_str=f'''Mann-Kendall Test shows that FDD trend is **{trend_mk}**.
\n\n**P-value** is **{pval_mk:.3e}**, which is **{alpha_compr}**
than **0.05**, meaning that the trend in the data is {p_str}'''

И в таком виде результаты теста выводятся на страницу:

Ну и напоследок можно посмотреть указанное число наиболее мягких и наиболее суровых зим ранжированных по FDD. В пятерку наиболее мягких попали сезоны последних 30 лет. В пятерку наиболее суровых попали зимы не позже 1970х.

Заключение

Анализ столетних наблюдений температуры воздуха, показал что зимы в Алматы становятся мягче.

Весь дашборд уложился в менее чем 300 строк питоновского файла, значительную часть из которых занимает пояснительный текст. На составление основной части веб-приложения ушло несколько часов, ну и потом периодически что-то добавлял. Для сравнения, до знакомства со Streamlit я делал визуализацию с Bokeh+Flask, и на составление одного графика, когда я только начал изучать Bokeh, ушло полдня.

В будущем буду расширять дашборд новыми элементами и локациями. Добавлю сравнение средних максимумов и минимумов за период. Еще планирую добавить другие параметры, такие как ветер и осадки. Ну и если руки дойдут, может и моделирование какое-нибудь.

Визуализация и анализ зимних температур Алматы за последние сто лет на Streamlit