Предсказание горимости территорий методами машинного обучения на основе географических данных08.07.2022 17:32

Автор статьи: Ярослав Найчук, выпускник OTUS.

Статья написанна Ярославом на основе выпускного проекта, в котором показано, как методы ML можно применять в самых разных областях.

Добрый день! Сегодня я хотел бы рассказать о том, как с помощью средств машинного обучения можно предсказать горимость регионов. Что же такое горимость? В данной работе она рассматривается как отношение площади гарей за год к площади исследуемого региона.

Целью работы не ставилась идентификация самих выгоревших площадей, напротив, это попытка составить оптимальную модель предсказания, используя простые географические данные региона.

В работе рассмотрены Атырауская, Западно-Казахстанская области Республики Казахстан, а также Астраханская, Саратовская, Волгоградская области и Республика Калмыкия России.

Данные по горимости и некоторые другие географические данные любезно предоставлены авторами данной работы: Шинкаренко С.С., Дорошенко В.В., Берденгалиева А.Н. Динамика площади гарей в зональных ландшафтах юго-востока европейской части России // Известия Российской академии наук. Серия географическая. 2022. Т. 86. № 1. С. 122–133. DOI 10.31857/S2587556622010113.

Какие данные использовались ещё?

• Климатические данные. Взяты с сайта https://crudata.uea.ac.uk/cru/data/hrg/

• Данные по поголовью скота — https://stat.gov.kz/region/list, https://www.gks.ru/dbscripts/munst/munst.htm

• NDVI — Вега-science http://sci-vega.ru/
  Loupian E., Burtsev M., Proshin A., Kashnitskii A., Balashov I., Bartalev S., Konstantinova A., Kobets D., Radchenko M., Tolpin V., Uvarov I. 2022. Usage Experience and Capabilities of the VEGA-Science System // Remote Sensing. Vol. 14. No. 1. P. 77. DOI: 10.3390/rs14010077

• Границы муниципальных образований.

Рассмотрим признаки более детально:

• flammability — динамика годовой горимости по регионам. Наша целевая переменная.

• NDVI — динамика нормализованного относительного индекса растительности по регионам. Данный признак рассчитан как сумма пикселей со значениями  NDVI в пределах региона.

• precipitation — динамика кол-ва осадков.

• t_max — среднегодовая динамика максимальных суточных температур.

• livestock — динамика поголовья скота, представленная в условных единицах.

• region — исследуемый район

Что же такое NDVI? NDVI (Normalized Difference Vegetation Index) — нормализованный относительный индекс растительности. Он отражает плотность растительности в определенной точке изображения, которая равна разнице интенсивностей отраженного света в красном и инфракрасном диапазоне, деленной на сумму их интенсивностей. Формула NDVI выглядит так:

Исходя из этого индекса можно сделать предположение, что чем меньше значение NDVI, тем биомасса более склонна к горению. Как видно из рисунка выше, сочные части растения имеют высокий индекс, они более влажные; отмершие или отсохшие части имеют более низкий индекс и более склонны к горению.

Изначально некоторые данные были представлены не в форме табличных данных, а в форме растровых изображений с географической привязкой (GeoTIFF). Для представления данных в нужной форме была использована программа QGIS для работы с геоданными. Расчеты по растрам проводились с помощью инструмента «Калькулятор растров». Далее для расчета статистики по регионам, был использован инструмент «Зональная статистика». На получившийся растр как бы накладывался слой с полигонами муниципальных районов, и для каждого района мы получили годовую статистику.  (см. рис)

В итоге проведенных расчетов получается векторный слой со следующей структурой  таблицы атрибутов (см рис.).

Каждая строка здесь — это полигон района, а поля »2001»,»2002» и др. — некоторые средние значения на конец года. Эта таблица не совсем подходит для анализа данных, т.к. следует представить каждый год и каждый регион отдельной строкой. Была написана следующая функция на Python, преобразующая данные в нужный вид

def parse_csv(path, col_name, start_index=0):
    df = pd.read_csv(path, decimal=',')
    list_of_df = []
    for i in range(len(df)):
        new_df = df.iloc[i,:].to_frame(name=col_name).reset_index()
        if col_name == 'flammability':
            new_df = df.iloc[i,:].to_frame(name="flammability").reset_index()
            new_df['region'] = new_df.iloc[1,1]
            new_df['OSM_ID'] = new_df.iloc[0,1]
            new_df = new_df.iloc[3:,:]
            list_of_df.append(new_df)
        else:
            new_df['OSM_ID'] = new_df.iloc[0,1]
            new_df = new_df.iloc[1:,:]
            list_of_df.append(new_df)
    
    df = pd.concat(list_of_df).reset_index()
    df = df.rename(columns={"index": "year"})
    df = df.drop(['level_0'], axis=1)
    return df

В итоге получаем таблицу следующего вида:  

Получившиеся датасеты сцеплялись по полю OSM_ID функцией pd.merge с использованием itertools.

import functools as ft
df = ft.reduce(lambda left, right: pd.merge(left, right, on=['OSM_ID','year']), list_of_data)

Также были проведены EDA и предобработка данных, добавились новые признаки.

• Значения NDVI и осадков были смещены на год назад. Это сделано потому, что исходные данные отражают информацию на конец года. Пропущенные значения, которые образовывались из-за смещения, были заполнены средним по региону.

• Был добавлен новый признак, который отражал динамику количества осадков на начало года, с января по апрель (precip_beggining_y)

• Муниципальные районы (категориальная переменная) были переделаны в новые признаки: широта и долгота их географического центра (N — координаты широты, E — координаты долготы)

• Целевая переменная была прологарифмирована .

• Nan«ы заполним значением »-9999», пропущенные данные присутствуют только в данных по поголовью.

• Поля «year» и поля с ненужной информацией были удалены.

Теперь посмотрим на матрицу корреляций признаков

Целевая переменная имеет слабую корреляцию с другими признаками. Также можно заметить, что координаты широты обратно коррелируют с динамикой максимальных суточных температур, что встраивается в картину закона о «Широтной зональности» (когда наблюдается закономерное изменение физико-географических процессов, компонентов и комплексов геосистем от экватора к полюсам).

Посмотрим также распределение целевой переменной.

Здесь более темные цвета отражают большую горимость

Попробуем протестировать классические модели машинного обучения. Для начала оценим линейные модели с регуляризацией и с применением различных scaler«ов. Качество можно увидеть в таблице ниже.

SVM со следующей сеткой показал следующие результаты:

param_grid = {
    'C': [1, 5, 100, 200, 1000],
    'epsilon': [0.01, 0.1, 0.05, 0.0003, 1, 0.2, 5, 10],
    'gamma': [0.001, 0.1, 1, 5, 10, 100]
}

R2: 0.465

Best params:   {'C': 100, 'epsilon': 0.05, 'gamma': 0.1}

CatBoost «из коробки» с разными scaler«ами показал:

LightGBM со значениями по умолчанию RMSE был равен 0.644, а R2: 0.563

После использования библиотек градиентного бустинга, всё же хочется посмотреть на «feature_importance». Здесь мы имеем следующие значения:  

1. previous_NDVI: 24.23

2. previous_precipitation: 20.7

3. precip_beggining_y: 15.91

4. N: 10.7

5. E: 9.82

6. t_max: 9.3

7. livestock: 9.18

Первое место занимает показатель NDVI, смещенный на год назад. Это можно объяснить тем, что на горимость влияет то, насколько много отмершей биомассы накопилось за прошлый год; также немаловажны осадки за прошлый год и осадки на начало года.

Любопытное наблюдение: для модели KNeighborsRegressor для k=7, модель показала результат R-квадрат равный 0.519. Как такая простая модель могла почти догнать градиентный бустинг? Всё объясняется тем, что метод k-ближайших соседей хорошо подходит для решения географических задач, ведь с большой вероятностью можно сказать, что природные условия некоего географического объекта сходны с его соседями.

Заключение: наилучший результат показал CatBoost со значениями по умолчанию и использованием MinMaxScaler«а, с R2 равным 0.601. Отдельного внимания заслуживает метод k-ближайших соседей, как простая и менее ресурсозатратная альтернатива градиентному бустингу (в рамках данной задачи).

В конце статьи хочу порекомендовать всем заинтересованным двухдневный бесплатный интенсив от OTUS, на котором вас ждет глубокое погружение в область рекомендательных систем. Уже в первый день вы изучите наиболее популярные подходы, которые существуют для формирования рекомендаций и реализуете один из них своими руками. Зарегистрироваться на интенсив можно по ссылке ниже.

© Habrahabr.ru