План исследования
Более детальное описание тут.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
from scipy import stats
import operator
%matplotlib inline
# import warnings
# warnings.filterwarnings("ignore")
ORIGIANL_DATA_FILE = './../data/project_brazil/sudeste.csv'
PROCESSED_DATA_FILE = './../data/project_brazil/processed.csv'
RANDOM_STATE = 14
Имеется набор данных почасовых наблюдений погодных условий от 122 погодных станций Юго-Восточного региона Бразилии. В регион входят штаты: Рио-де-Жанейро, Сан-Паулу, Минас-Жерайс и Эспириту-Санту. Данные собираются с 2000 года (однако не все станции начали работы по сбору именно с этого года). Источником данных является INMET (Национальный Метеорологический Институт - Бразилия). Описание и ссылку на скачивание данных можно найти на странице Kaggle Datasets.
Описания признаков в датасете
wsid (Numeric) - Уникальный идентификатор погодной станции
wsnm (String) - Наименование погодной станции
elvt (Numeric) - Высота расположения погодной станции
lat (Numeric) - Широта месторасположения погодной станции (град.)
lon (Numericw) - Долгота месторасположения погодной станции (град.)
inme (String) - INMET-код станции для области
city (String) - Город
prov (String) - Штат (провинция)
mdct (DateTime) - Дата и время наблюдения
date (DateTime) - Дата наблюдения
yr (Numeric) - Год наблюдения (2000-2016)
mo (Numeric) - Месяц наблюдения (1-12)
da (Numeric) - День наблюдения (1-31)
hr (Numeric) - Час наблюдения (0-23)
prcp (Numeric) - Количество осадков за последний час (мм)
stp (Numeric) - Давление воздуха (мгновенное) (гПа)
smax (Numeric) - Максимальное давление воздуха за последний час (гПа)
smin (Numeric) - Минимальное давление воздуха за последний час (гПа)
gbrd (String) - Солнечное излучение (кДж/м2)
temp (Numeric) - Температура воздуха (мгновенная) (град. Цельсия)
dewp (Numeric) - Температура точки росы (мгновенная) (град. Цельсия)
tmax (Numeric) - Максимальная температура за последний час (град. Цельсия)
dmax (Numericw) - Максимальная температура точки росы за последний час (град. Цельсия)
tmin (Numeric) - Минимальная температура за последний час (град. Цельсия)
dmin (Numeric) - Минимальная температура точки росыза последний час (град. Цельсия)
hmdy (Numeric) - Относительная влажность (мгновенная) (%)
hmax (Numeric) - Максимальная относительная влажность за последний час (%)
hmin (Numeric) - Минимальная относительная влажность за последний час (%)
wdsp (String) - Скорость ветра (м/с)
wdct (Numeric) - Направление ветра (град) (0-360)
gust (String) - Порывы ветра (м/с)
Данный датасет может быть использован для опробования моделей машинного обучения в области прогнозирования погоды. Явно выделенной целевой переменной нет, а на Kaggle предлагается прогнозировать количество осадков или температуру. Мы попробуем спрогнозировать осадки, т.е. целевой меткой будет являться признак prcp. Подробнее последует далее.
Необходимо учесть тот факт, что данные по осадкам хранятся для того часа, для которого актуальны и прочие наблюдения параметров. Для прогнозирования осадков на час, следующий за наблюдаемыми данными, в тренировочном датасете надо выполнить сдвиг значений целевой переменной на один шаг "вперед".
Учитывая, что размер датасета значительный (1.72 Гб), сделаем предобработку файла данных, чтобы не загружать весь датасет в память. В простом предварительном взгляде на данные было выявлено, что примерно в 86% исходных данных отсутствует значение признака **prcp**. Выполним фильтрацию таких данных, оставив для дальнейшего анализа и обработки только корректные.
Также стоит отметить, что, несмотря на присутствие в данных признаков с датой и временем, датасет не стоит расматривать как упорядоченный (какого-то тренда в погоде с течением времени не может быть в такой короткий временной отрезок, как 10-20 лет). Соответственно, эти признаки несут категориальный характер.
%%time
with open(ORIGIANL_DATA_FILE, 'r', encoding='utf-8') as in_file, open(PROCESSED_DATA_FILE, 'w', encoding='utf-8') as out_file:
total_lines = 0
processed_lines = 0
print('Обработано строк:', end=' ', flush=True)
# read full set
for line in in_file:
total_lines += 1
# write dataset header line
if total_lines == 1:
out_file.write(line)
continue
# progress indication
if total_lines % 1000000 == 0:
print('{}M'.format(total_lines // 1000000), end=' ', flush=True)
# serch 'prcp' field
prcp_field = line.split(',')[14]
if not prcp_field:
continue
# write line
processed_lines += 1
out_file.write(line)
print()
print('Всего обработано строк (с заголовком): {}'.format(total_lines))
print('Отфильтрованные строки: {}'.format(processed_lines))
Загрузим предварительно обработанные данные.
df = pd.read_csv(PROCESSED_DATA_FILE)
print('Размер предварительно обработанного датасета: {}'.format(df.shape))
df.head()
Посмотрим на основную информацию о датасете.
df.info()
Посмотрим на характеристики распределения категориальных признаков. Сразу приведем целочисленные признаки wsid, yr, mo, da, hr (которые являются категориальными) к строковому типу для корректной обработки методом describe().
int_features = ['wsid', 'yr', 'mo', 'da', 'hr']
for feat in int_features:
df[feat] = df[feat].astype('str')
df.describe(exclude=['float64']).T
В выборке представлены 4 провинции и 177 городов. Аномальных значений в дате и времени нет. По описанию и практически совпадающим числам уникальных значений признаков **wsid**, **wsnm** и **inme** можно предположить, что данные признаки дублируют друг друга. Проверим это предположение.
Идея состоит в следующем. Используя метод *groupby()*, можно сгруппировать значения одного признака по уникальным значениям другого. Затем с помощью метода *describe()* можно посмотреть на число уникальных значений второго признака в группах. Если на всех позициях стоят единицы, то можно утверждать о полном соответствии в значениях двух признаков.
wsid_groupedby_wsnm_descr = df.groupby('wsnm')['wsid'].describe()
wsid_wsnm_non_unique = wsid_groupedby_wsnm_descr[wsid_groupedby_wsnm_descr['unique'] != 1].shape[0]
print('Число несоответствий между признаками wsid и wsnm: {}'.format(wsid_wsnm_non_unique))
wsid_groupedby_inme_descr = df.groupby('inme')['wsid'].describe()
wsid_inme_non_unique = wsid_groupedby_inme_descr[wsid_groupedby_inme_descr['unique'] != 1].shape[0]
print('Число несоответствий между признаками wsid и inme: {}'.format(wsid_inme_non_unique))
Рассмотрим несоответствие детальней. Найдем значение inme той станции, по которой происходят расхождения, и посмотрим на уникальные значения координат, соответствующие ей.
non_unique_inme = wsid_groupedby_inme_descr[wsid_groupedby_inme_descr['unique'] != 1].index[0]
print(df[df['inme'] == non_unique_inme]['lat'].value_counts())
print()
print(df[df['inme'] == non_unique_inme]['lon'].value_counts())
Видно, что координаты этой проблемной станции принимают во всем датасете лишь одно значение. Значит, предположение о дублировании друг друга признаками wsid, wsnm, inme оказывается верным.
Выделим вещественные признаки и посмотрим на характеристики их распределения.
df.describe(include=['float64']).T
Можно сделать предварительные выводы по признакам:
Посомтрим на количество пропусков данных для признаков и отобразим их в порядке убывания.
null_counts = df.isnull().sum()
null_counts[null_counts != 0].sort_values(ascending=False)
Из описания датасета можно увидеть, что признаки mdct (дата и время наблюдения) и date (дата наблюдения) избыточны, их полностью заменяют следующие признаки: date, yr, mo, da, hr (дата, год, месяц, день, час соответственно). Сконвертируем тип признаков mdct и date во внутренний тип pandas.datetime и проверим, что значения этих признаков соответствуют друг другу.
df['mdct'] = df['mdct'].apply(pd.to_datetime)
df['date'] = df['date'].apply(pd.to_datetime)
def datetime_accordance_check(row):
return ((row['yr'] == row['mdct'].year) & (row['yr'] == row['date'].year) &
(row['mo'] == row['mdct'].month) & (row['mo'] == row['date'].month) &
(row['da'] == row['mdct'].day) & (row['da'] == row['date'].day) &
(row['hr'] == row['mdct'].hour))
df.apply(datetime_accordance_check, axis=1).value_counts()
Видим, что рассматриваемые признаки полностью соответствуют друг другу. Избыточные признаки можно будет удалить из датасета.
Рассмотрим корреляцию вещественных признаков, исключая целевой. Отобразим распределение и топ-20 коррелирующих пар признаков по двум методам: Пирсона и Спирмана (оценивающие силу линейной и монотонной взаимосвязей соответственно).
numeric_columns = [col for col in df if (df[col].dtype == 'float64')]
pearson_corr = df[numeric_columns].corr(method='pearson')
spearman_corr = df[numeric_columns].corr(method='spearman')
no_target_unique_pearson_corr = dict()
no_target_unique_spearman_corr = dict()
for (row, r_name) in enumerate(pearson_corr.index):
for (col, c_name) in enumerate(pearson_corr.columns[row+1:]):
if('prcp' in [r_name, c_name]):
continue
no_target_unique_pearson_corr['{}-{}'.format(r_name, c_name)] = abs(pearson_corr.iloc[row, col])
for (row, r_name) in enumerate(spearman_corr.index):
for (col, c_name) in enumerate(spearman_corr.columns[row+1:]):
if('prcp' in [r_name, c_name]):
continue
no_target_unique_spearman_corr['{}-{}'.format(r_name, c_name)] = abs(spearman_corr.iloc[row, col])
pearson_corr_df = pd.DataFrame(sorted(no_target_unique_pearson_corr.items(), key=operator.itemgetter(1), reverse=True))
pearson_corr_df.columns = ['', 'pearson corr']
spearman_corr_df = pd.DataFrame(sorted(no_target_unique_spearman_corr.items(), key=operator.itemgetter(1), reverse=True))
spearman_corr_df.columns = ['', 'spearman corr']
corr_df = pd.concat([pearson_corr_df, spearman_corr_df], axis=1)
corr_df.hist()
plt.show()
corr_df = pd.concat([pearson_corr_df.head(20), spearman_corr_df.head(20)], axis=1)
corr_df
Видим, что оба метода на топ-результатах дают очень схожие показатели. Можно предположить, что корреляция между этими признаками есть, и взаимосвязь линейна. Но зависимость не объясняется физической природой признаков.
Рассмотри корреляцию с целевым признаком.
target_pearson_corr = dict(pearson_corr.apply(abs).loc['prcp', :])
target_pearson_corr.pop('prcp')
target_spearman_corr = dict(spearman_corr.apply(abs).loc['prcp', :])
target_spearman_corr.pop('prcp')
target_pearson_corr_df = pd.DataFrame(sorted(target_pearson_corr.items(), key=operator.itemgetter(1), reverse=True))
target_pearson_corr_df.columns = ['', 'pearson corr']
target_spearman_corr_df = pd.DataFrame(sorted(target_spearman_corr.items(), key=operator.itemgetter(1), reverse=True))
target_spearman_corr_df.columns = ['', 'spearman corr']
target_corr_df = pd.concat([target_pearson_corr_df, target_spearman_corr_df], axis=1)
target_corr_df
Тут мы уже не наблюдаем схожесть в значениях коэффициентов корреляции. На оценку по методу Пирсона скорее всего повлияла ненормальность распределения, а по методу Спирмана вероятно, что оценки завышены. В любом случае, можно говорить о некоторой монотонной зависимости целевого признака от ряда других, некоторые объясняются физической природой признаков (например, осадки от влажности и точки росы).
Наконец, посмотрим на распределение целевого признака. Построим график в разных масштабах оси ординат.
_, axes = plt.subplots(nrows=1, ncols=3, figsize=(15, 5))
sns.distplot(df['prcp'], ax=axes[0])
sns.distplot(df['prcp'], ax=axes[1])
axes[1].set_ylim(0, 0.05)
sns.distplot(df['prcp'], ax=axes[2])
axes[2].set_ylim(0, 0.01)
plt.show()
Очевидно, что распределение целевого признака далеко от номрального. Тем не менее, проведем статистические тесты на нормальность и скошенность распределения. Для оценки критерия Шапиро-Уилка в модуле scipy.stats существует метод shapiro(). У него есть ограничение на размер выборки - 5000 (в документации указано, что при выборках большего размера точность p-value не гарантируется). Проведем 10 случайных взятий подвыборки заданного размера и найдем среднее значение p-value на этих разбиениях.
shapiro_p_values = []
for split_num in range(10):
np.random.seed(RANDOM_STATE + split_num)
rand_split = np.random.randint(0, df.shape[0], size=5000)
shapiro_p_values.append(stats.shapiro(df.loc[rand_split, 'prcp'])[1])
print('Среднее значение p-value критерия Шапиро-Уилка на случайных разбиениях: {}'.format(np.mean(shapiro_p_values)))
print('Значение p-value теста на скошенность распределения: {}'.format(stats.skewtest(df['prcp'])[1]))
Полученное среднее p-value критерия Шапиро-Уилка меньше уровня значимости в 0.05, значит распределение не нормальное. Тест на скошенность распределения также не прошел.
Посмотрим внимательней на правый "хвост" распределения.
df[df['prcp'] > 75]['prcp'].value_counts().iloc[-25:]
Осадки интенсивностью в районе 100 мм/час - это очень сильные ливни. И хоть такие наблюдения значимы, но их слишком мало в общей массе для какого-то вклада в модель. Выполним фильтрацию нашей целевой переменной по 99% квантилю, чтобы избавиться от таких экстремальных единичных значений.
prcp_quantile = df['prcp'].quantile(0.99)
prcp_filtered = df[df['prcp'] < prcp_quantile]['prcp']
plt.figure(figsize=(15, 5))
sns.distplot(prcp_filtered)
После фильтрации по квантилю обнаружилась следующая особенность целевого признака - он принимает дискретные значения (видимо, связанные с ограничениями на измерения в погодных станциях).
Учитывая все эти факты, принимаем решение перейти от казалось бы предполагаемой задачи регрессии к задаче многоклассовой классификации. Подчерпнув информацию о принятой классификации интенсивности дождей (точнее, об отсутствии единой таковой), выделим "на глаз" следующие классы целевой переменной. - Отсутствие осадков: $prcp = 0.0$ - Небольшой дождь: $prcp\in(0.0; 1.2]$ - Умеренный дождь: $prcp\in(1.2; 2.5]$ - Сильный дождь: $prcp > 2.5$
Начнем визуальный анализ с подгруппы категориальных признаков, отражающих временные характеристики: yr, mo, da, hr. Построим гистограммы распределений.
# переведем тип обратно из str
int_features = ['wsid', 'yr', 'mo', 'da', 'hr']
for feat in int_features:
df[feat] = df[feat].astype('int64')
time_features = ['yr', 'mo', 'da', 'hr']
_, axes = plt.subplots(nrows=2, ncols=2, figsize=(14, 8))
for feat_idx, feat in enumerate(time_features):
ax = axes[feat_idx // 2, feat_idx % 2]
df[feat].hist(ax=ax)
ax.set_title(feat)
Видны следующие факты:
Построим графики зависимости среднего числа осадков от года и месяца
_, axes = plt.subplots(nrows=1, ncols=2, figsize=(14, 4))
df[['prcp', 'yr']].groupby('yr').mean().plot(ax=axes[0])
df[['prcp', 'mo']].groupby('mo').mean().plot(ax=axes[1])
Видны следующие факты:
Посмотрим на распределение данных по погодным станциям. В первичном анализе данных мы выяснили, что три варианта идентификаторов станций (wsid, wsnm, inme) дублируют друг друга. Будем испольовать признак wsid, как целочисленный.
wsid_int = df['wsid'].astype('int64')
print('Уникальные занчения признака wsid: {}'.format(np.unique(wsid_int)))
wsid_int.hist(bins=len(np.unique(wsid_int)))
Количество данных от разных станций в датасете сильно разнится.
Рассмотрим теперь вещественные признаки. Из первичного анализа данных мы знаем, что в ряде этих признаков присутствуют пропуски. Заполним их соответствующими средними значениями.
real_features = [feat for feat in df if (df[feat].dtype == 'float64')]
real_features = list(set(real_features) - {'lat', 'lon', 'prcp'})
df_real_features = df[real_features].copy()
df_real_features.fillna(df_real_features.mean(), inplace=True)
Построим гистограммы распределений и "ящики с усами".
_, axes = plt.subplots(nrows=6, ncols=3, figsize=(16, 20))
for feat_idx, feat in enumerate(df_real_features):
ax = axes[feat_idx // 3, feat_idx % 3]
#df_real_features[feat].hist(ax=ax)
sns.distplot(df_real_features[feat], ax=ax)
#ax.set_title(feat)
_, axes = plt.subplots(nrows=6, ncols=3, figsize=(16, 20))
for feat_idx, feat in enumerate(df_real_features):
ax = axes[feat_idx // 3, feat_idx % 3]
sns.boxplot(df_real_features[feat], ax=ax)
Видны следующие факты:
Рассмотрим теперь графически корреляцию вещественных признаков между собой и с целевым признаком по методу Спирмана.
df_real_features_with_target = pd.concat([df_real_features, df['prcp']], axis=1)
corr = df_real_features_with_target.corr(method='spearman')
# маска для отсечения верхнего треугольника корреляционной матрицы
corr_mask = np.zeros_like(corr, dtype=np.bool)
corr_mask[np.triu_indices_from(corr_mask)] = True
plt.figure(figsize=(10, 10))
sns.heatmap(corr, mask=corr_mask, annot=True, fmt='0.1f', linewidths=0.5, cmap='YlGnBu')
plt.show()
Как было отмечено ранее, видна достаточно высокая монотонная корреляция между всеми признаками. Это не очень хорошо, но возможно это связано с большим числом нулевых значений во всех признаках. Также есть признаки, коррелирующие по своей природе. Например, мгновенные, максимальные и минимальные значения какого-то показателя явно будут сильно коррелирующими.
Построим графики попарных зависимостей коррелирующих по своей природе признаков. Это мгновенные, минимальные и максимальные значения для: давления, температуры, температуры точки росы, влажности.
pairplot_features = ['stp', 'smax', 'smin']
sns.pairplot(df_real_features_with_target[pairplot_features])
plt.show()
pairplot_features = ['temp', 'tmax', 'tmin']
sns.pairplot(df_real_features_with_target[pairplot_features])
plt.show()
pairplot_features = ['dewp', 'dmax', 'dmin']
sns.pairplot(df_real_features_with_target[pairplot_features])
plt.show()
pairplot_features = ['hmdy', 'hmax', 'hmin']
sns.pairplot(df_real_features_with_target[pairplot_features])
plt.show()
Исключая из графиков нулевые выбросы, явно просматривается линейный тренд, для всех групп, кроме значений влажности, как это и обнаружилось в первичном анализе данных.
Построим теперь графики попарных зависимостей между целевым признаком и следующими (тут мы исключаем направление ветра и минимальные/максимальные значения температуры, влажности, давления, температуры точки росы): elvt, stp, gbrd, temp, dewp, hmdy, wdsp, gust. Для читаемости разделим построение графиков на 2 шага.
pairplot_features = ['prcp', 'elvt', 'stp', 'gbrd', 'temp']
sns.pairplot(df_real_features_with_target[pairplot_features])
plt.show()
pairplot_features = ['prcp', 'dewp', 'hmdy', 'wdsp', 'gust']
sns.pairplot(df_real_features_with_target[pairplot_features])
plt.show()
Явного линейного тренда не видно ни на одной из попарных зависимостей
Подводя итог первичному анализу данных, сведем воедино все найденные факты и предположения.
Для выбора метрики сначала посмотрим на баланс классов целевой переменной. Выведем долю каждого класса.
prcp = df['prcp']
def prcp_transform(val):
if val == 0.0:
return 0
elif (val > 0.0) & (val <= 1.2):
return 1
elif (val > 1.2) & (val <=2.5):
return 2
elif (val > 2.5):
return 3
prcp = prcp.apply(prcp_transform)
prcp.value_counts(normalize=True)
Можно видеть, что классы сильно несбалансированы. Обычной выбором метрики для задачи классификации с сильным дисбалансом классов является метрика **ROC-AUC**. В отличие от более простых метрик (например, доли верных ответов), ROC-AUC учитывает как TPR (True Positive Rate), так и FPR (False Positive Rate). Значит, она не чувствительна к дисбалансу классов. Также эта метрика позволяет дать оценку качества классификации, основываясь на вероятностых предположениях принадлежности к классу, не привязываяь к какому-то конкретному порогу классификации.
Т.к. классификация в нашей задаче многоклассовая, то необходимо выбрать способ усреднения метрики. После изучения возможных способов в документации sklearn, выберем метод **'macro'** (итоговая оценка метрики усредняется без учета веса класса, что дает более пессимистичную оценку, чем методы 'micro' и 'weighted').
Для нашей задачи будем рассматривать следующие модели:
Подготовим датасет для дальнейшей работы моделей. Начнем с целевой переменной. Выполним фильтрацию всего датасета по 99% квантилю целевой перемнной. Выполним сдвиг ее значений на одну позицию вперед, чтобы прогнозировать осадки на час, следующий за наблюдаемым, и выкинем из датасета первую строчку, у которой значение признака prcp станет пустым.
prcp_quantile = df['prcp'].quantile(0.99)
df = df[df['prcp'] < prcp_quantile]
df['prcp'] = df['prcp'].shift(1)
df.drop(index=0, inplace=True)
Выполним кодирование значений целевой переменной и выделим ее в отдельный объект.
y = df['prcp'].apply(prcp_transform).astype('int')
Создадим словарь обработанных признаков, для дальнейшего их объединения в итоговый датасет.
processed_features = dict()
Т.к. одной из моделей у нас является линейная, то категориальные признаки нам необходимо закодировать, а численные - отмасштабировать.
Обозначим схему преобразования данных. В нашей задаче у нас отсутствует тестовая выборка, валидацию моделей будем проводить на отложенной части выборки. Но преобразовывать признаки мы будем на всей имеющейся выборке сразу. В общем случае это неправильно, т.к. из валидационной части может "просочиться" информация в обучающую часть (например, учтется масштаб всех значений признаков, будут закодированы все значения категориальных признаков, и т.д.). Мы же закроем глаза на это допущение.
Признак wsid. Является категориальным, выполним One-Hot-Encoding с выбрасыванием первой колонки (для того, чтобы не возникало зависимости в новых категориальных признаках). Признаки wsnm, inme не используем как дублирующие.
wsid_ohe = pd.get_dummies(df['wsid'], prefix='wsid_ohe', drop_first=True)
processed_features['wsid_ohe'] = wsid_ohe
Признаки elvt, lat, lon являются численными, выполним масштабирование.
stnd_scaler = StandardScaler()
elvt_scaled = stnd_scaler.fit_transform(df[['elvt']])
processed_features['elvt_scaled'] = pd.DataFrame(elvt_scaled, columns=['elvt_scaled'])
lat_scaled = stnd_scaler.fit_transform(df[['lat']])
processed_features['lat_scaled'] = pd.DataFrame(lat_scaled, columns=['lat_scaled'])
lon_scaled = stnd_scaler.fit_transform(df[['lon']])
processed_features['lon_scaled'] = pd.DataFrame(lon_scaled, columns=['lon_scaled'])
Признаки city, prov не будем использовать в моделях, как малоинформативные (географическую информацию мы уже добавили в виде широты и долготы), но потенциально "раздувающие" датасет (они категориальны и придется кодировать).
Как было обозначено ранее, из признаков даты и времени мы будем использовать только mo (месяц). Пояснения см. в п. 4. Несмотря на целочисленность признака, он является конечно же категориальным, так что выполним One-Hot-Encoding.
mo_ohe = pd.get_dummies(df['mo'], prefix='mo_ohe', drop_first=True)
processed_features['mo_ohe'] = mo_ohe
Для stp, temp, dewp (давления, температуры воздуха, температуры точки росы) возьмем только мгновенные их значения (максимальные и минимальные значения линейно зависимы, см. п. 4). Воспользуемся масштабированием. Пустые значения заменим средними.
stp_scaled = stnd_scaler.fit_transform(df[['stp']])
processed_features['stp_scaled'] = pd.DataFrame(stp_scaled, columns=['stp_scaled'])
temp_scaled = stnd_scaler.fit_transform(df[['temp']].fillna(df[['temp']].mean()))
processed_features['temp_scaled'] = pd.DataFrame(temp_scaled, columns=['temp_scaled'])
dewp_scaled = stnd_scaler.fit_transform(df[['dewp']].fillna(df[['dewp']].mean()))
processed_features['dewp_scaled'] = pd.DataFrame(dewp_scaled, columns=['dewp_scaled'])
Для gbrd (солнечное излучение) заполним пропуски и выполним масштабирование.
gbrd_scaled = stnd_scaler.fit_transform(df[['gbrd']].fillna(df[['gbrd']].mean()))
processed_features['gbrd_scaled'] = pd.DataFrame(gbrd_scaled, columns=['gbrd_scaled'])
Показатели относительной влажности (мгновенная, максимальная и минимальная) используем все, т.к. линейной зависимости в них не наблюдалось, и данные могут быть полезными для прогнозирования осадков. Воспользуемся масштабированием. Пустые значения заменим средними.
hmdy_scaled = stnd_scaler.fit_transform(df[['hmdy']])
processed_features['hmdy_scaled'] = pd.DataFrame(hmdy_scaled, columns=['hmdy_scaled'])
hmin_scaled = stnd_scaler.fit_transform(df[['hmin']].fillna(df[['hmin']].mean()))
processed_features['hmin_scaled'] = pd.DataFrame(hmin_scaled, columns=['hmin_scaled'])
hmax_scaled = stnd_scaler.fit_transform(df[['hmax']].fillna(df[['hmax']].mean()))
processed_features['hmax_scaled'] = pd.DataFrame(hmax_scaled, columns=['hmax_scaled'])
Показатели ветра wdsp, gust (скорость, порывы) отмасштабируем и заполним пропуски. Направление ветра wdct разобьем на 10-градусные сектора, 36-ой сектор заменим нулевым, и выполним One-Hot-Encoding.
wdsp_scaled = stnd_scaler.fit_transform(df[['wdsp']].fillna(df[['wdsp']].mean()))
processed_features['wdsp_scaled'] = pd.DataFrame(wdsp_scaled, columns=['wdsp_scaled'])
gust_scaled = stnd_scaler.fit_transform(df[['gust']].fillna(df[['gust']].mean()))
processed_features['gust_scaled'] = pd.DataFrame(gust_scaled, columns=['gust_scaled'])
wdct_processed = df['wdct'].fillna(df[['wdct']].mean())
wdct_processed = wdct_processed.apply(lambda x: x // 10).apply(lambda x: x if x != 36.0 else 0.0)
wdct_ohe = pd.get_dummies(wdct_processed, prefix='wdct_ohe', drop_first=True)
processed_features['wdct_ohe'] = wdct_ohe
Наконец, соберем все признаки воедино.
for name_df, part_df in processed_features.items():
processed_features[name_df] = part_df.reset_index(drop=True)
X = pd.concat(processed_features.values(), axis=1)
print('Размер обработанной выборки: {}'.format(X.shape))
Выделим обучающую и валидационную части выборки. Ранее уже упоминалось, что упорядочивание выборки по времени не имеет смысла, так что воспользуемся случайным разбиением и выделим для валидации 30% всей выборки. Укажем параметр stratify для сохранения баланса классов в новых выборках.
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, shuffle=True, stratify=y, random_state=RANDOM_STATE)
print('Размер обучающей выборки:', X_train.shape)
print('Размер обучающего вектора целевой переменной:', y_train.shape)
print('Размер валидационной выборки:', X_test.shape)
print('Размер валидационного вектора целевой переменной:', y_test.shape)
Классификация осадков по интенсивности была подобрана "на глаз", основываясь на классификации из статьи.
2. http://www.zodchii.ws/books/info-756.html
"Количество солнечной энергии, поступающей за год на 1 м2 поверхности Земли, изменяется приблизительно от 3000 МДж/м2 на севере до 8000 MДж/м2 в наиболее жарких пустынных местах"
Что эквивалентно почасовым значениям: 340 кДж/м2 для северных регионов, 910 кДж/м2 для жарких пустынных регионов