Мат. моделирование / ДЗЗ
Практика 11

Прогнозирование развития территории. Плагин MOLUSCE

О чём эта тема

Как по двум разновременным картам покрытия земель обучить модель и построить прогноз развития территории — от подготовки данных MODIS Land Cover до обучения нейросети и проверки прогноза в плагине MOLUSCE.

Аннотация

Тема разбирает полный цикл прогнозирования изменений землепользования средствами плагина MOLUSCE. Сначала готовятся входные данные: скачивание и обрезка карт MODIS Land Cover (MCD12Q1) за несколько лет, работа с форматом HDF, приведение цифровой модели рельефа (ЦМР) и рассчитанной по ней крутизны склона к тому же разрешению и охвату, что и карты состояний. Далее в самом MOLUSCE загружаются растры начального и конечного состояний и дополнительные пространственные переменные, строится статистика и матрица переходов между классами, обучается модель (чаще всего — нейронная сеть, многослойный персептрон) с диагностикой по кривым обучения, запускается симуляция будущего состояния территории и, наконец, проводится валидация прогноза по трём вариантам коэффициента каппа. После изучения темы студент сможет подготовить согласованные пространственные данные для MOLUSCE, обучить модель прогноза и отличить успешное обучение от переобучения по виду кривых обучения.

Пререквизиты
  • Практика 6 — плагин SCP, работа с растровыми классификациями
  • Лекция 5 — оценка точности классификации, коэффициент каппа
Мотивация

Прогноз развития территории — это не гадание, а количественная модель, обученная на прошлой динамике: там, где лес уже сокращался в пользу застройки, модель предскажет продолжение этой тенденции. Коэффициент каппа при валидации показывает, насколько этому прогнозу можно доверять, — а не просто «выглядит правдоподобно».

Практика 11. Прогнозирование развития территории при помощи нейросетей. Плагин MOLUSCE

Презентация Данные к занятию

Знакомство с плагином MOLUSCE

MOLUSCE (Modules for Land Use Change Simulations) — это расширение для QGIS, предназначенное для анализа изменений землепользования на основе данных за два временных среза и построения прогноза будущих изменений.

Плагин позволяет использовать различные алгоритмы моделирования, включая методы на основе нейросетей.

Для установки плагина выполните следующие шаги:

  1. Перейдите в меню:

Модули → Управление модулями

  1. В строке поиска введите:

MOLUSCE

  1. Нажмите Установить.

После установки плагин будет доступен в меню Модули.

пиктограмма установленного плагина MOLUSCE на панели инструментов
Рис. П11.1 — пиктограмма установленного плагина MOLUSCE на панели инструментов

Подготовка исходных данных

В качестве входных данных расширению MOLUSCE требуются:

  • Карта исходного состояния (например, Начальная_дата)
  • Карта конечного состояния (например, Начальная_дата + N)
  • Карта для валидации (Начальная_дата + 2N) — используется для оценки качества прогноза

Каждая карта — это растровое изображение, в котором каждый пиксель содержит код состояния территории. Примеры:

  • 1 = лес
  • 2 = поле
  • 3 = городская застройка
  • и т.д.
пример растровой карты состояния территории — каждому пикселю присвоен код класса покрытия земель
Рис. П11.2 — пример растровой карты состояния территории — каждому пикселю присвоен код класса покрытия земель

N — глубина прогноза (измеряется в днях, неделях, месяцах, годах в зависимости от задачи).

Требования к входным данным

Все входные растр-карты должны иметь:

  • одинаковое пространственное разрешение
  • одинаковый географический охват
  • одинаковый размер в пикселях
требования MOLUSCE к пространственным характеристикам входных растров
Рис. П11.3 — требования MOLUSCE к пространственным характеристикам входных растров

Скачивание данных MODIS Land Cover (Collection 6)

MODIS Land Cover — это глобальный набор данных о земном покрове, предоставляемый продуктами:

  • MCD12Q1 — ежегодные растровые карты покрытия земель с пространственным разрешением 500 метров
  • MCD12C1 — агрегированные карты на крупной сетке

Они относятся к Collection 6 и охватывают весь земной шар, предоставляя информацию о типах землепользования по различным международным классификациям (например, IGBP, UMD, LAI и др.).

анимация: карты покрытия земель MODIS Land Cover за разные годы
Рис. П11.4 — анимация: карты покрытия земель MODIS Land Cover за разные годы

Данные MCD12Q1 могут быть использованы в качестве:

  • карты исходного состояния территории (например, за 2022 год)
  • карты конечного состояния территории (например, за 2023 год)

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

Таким образом, эти данные идеально подходят для анализа изменений покрытия территории во времени и построения прогноза в плагине MOLUSCE.

Дополнительно можно загрузить данные за третий временной период (год = начальная дата + 2N) для валидации модели, полученной с помощью MOLUSCE.

Где искать и как скачивать

Для поиска и загрузки данных используйте платформу Earthdata Search — открытый онлайн-инструмент для поиска и скачивания данных дистанционного зондирования. Он предоставляет доступ к более чем 400 миллионам записей из разных спутниковых миссий.

интерфейс Earthdata Search — поиск данных ДЗЗ по инструментам, датам и территории
Рис. П11.5 — интерфейс Earthdata Search — поиск данных ДЗЗ по инструментам, датам и территории
  1. Зарегистрируйтесь на Earthdata Search
  2. В вкладке Instruments Найдите MODIS
  3. Выберите интересующие годы, регион и формат данных
  4. Скачайте растр из коллекции MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V061 при необходимости преобразуйте его в нужную проекцию или обрежьте по нужной области

Проблемы открытия hdf файлов

Продукт MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V061 распространяется в формате HDF.

HDF (Hierarchical Data Format) — это универсальный формат для хранения и организации больших массивов научных данных. Он поддерживает сложную структуру данных, включая множество внутренних слоёв (subdatasets), что делает его мощным, но не всегда удобным в повседневной работе.

  • Формат HDF не является основным для большинства пользователей ГИС.
  • QGIS поддерживает работу с HDF, но может потребоваться:
  • обновление библиотек (например, GDAL),
  • ручной выбор нужного слоя при открытии файла,
  • или предварительная конвертация в более привычный формат (GeoTIFF и др.).

Для просмотра внутренних слоёв HDF-файла можно использовать инструмент командной строки OSGeoShell. Введите команду:

gdalinfo MCD12Q1.A2015001.h20v03.061.2022166031305.hdf
вывод команды gdalinfo: список вложенных слоёв (subdatasets) внутри HDF-файла
Рис. П11.6 — вывод команды gdalinfo: список вложенных слоёв (subdatasets) внутри HDF-файла

В результате будет выведен список всех вложенных слоёв (subdatasets), содержащихся в файле. Обычно нас интересует слой, соответствующий одной из классификаций землепользования, например:

SUBDATASET_1_NAME=HDF4_EOS:EOS_GRID:"MCD12Q1...":MOD12Q1:LC_Type1

Рекомендуется извлекать именно этот слой (например, LC_Type1), так как он соответствует классификации IGBP.

После определения нужного слоя его можно импортировать в QGIS как обычный растр или использовать gdal_translate для преобразования в GeoTIFF:

gdal_translate "HDF4_EOS:EOS_GRID:MCD12Q1...:MOD12Q1:LC_Type1" output.tif

Такой файл уже будет легко обрабатываться в QGIS и совместим с плагином MOLUSCE.

нужный слой добавлен в QGIS как обычный растровый слой
Рис. П11.7 — нужный слой добавлен в QGIS как обычный растровый слой

Настройка отображения

Откройте MCD12_User_Guide_V6.pdf или перейдите по ссылке

Для работы с данными Collection 6 MODIS Land Cover (MCD12Q1 и MCD12C1) используйте следующие классы, соответствующие классификации FAO-Land Cover Classification System (LCCS2):

Код Класс Описание
1 Пустынные земли ≥60% площади не покрыто растительностью (песок, скалы, почва) или вечными снегами/льдом (<10% растительности).
2 Постоянные снега и льды ≥60% площади покрыто снегом и льдом ≥10 месяцев в году.
3 Водные объекты ≥60% площади занято постоянными водными объектами.
9 Городские и застроенные земли ≥30% площади — искусственные поверхности (здания, асфальт, транспорт).
10 Густые леса Покрытие деревьями >60% (высота крон >2м).
20 Редкие леса Покрытие деревьями 10-60% (высота крон >2м).
25 Мозаика лесов и пахотных земель Мозаика мелкомасштабных посевов (40-60%) с естественным древесным покровом (>10%).
30 Естественные травянистые земли Преобладают травянистые растения (<2м). Покрытие ≥10%.
35 Мозаика травянистых и пахотных земель Мозаика посевов (40-60%) с естественной кустарниковой/травянистой растительностью.
36 Травянистые пахотные земли Преобладают культивируемые травянистые растения (<2м). Покрытие ≥60%.
40 Кустарниковые земли Покрытие кустарниками >60% (высота 1-2м).
255 Неклассифицировано Не классифицировано (отсутствие входных данных).

Для настройки легенды Слои -> ПКМ -> Свойства -> Стиль -> Одноканальное псевдоцвектное - > Подписи и значения из таблицы pdf

карта покрытия земель с настроенной легендой по кодам классов
Рис. П11.8 — карта покрытия земель с настроенной легендой по кодам классов

Обрезание растров

Скачанные карты покрытия земель и открытый векторный слой границы территории (например, Рязанская область.gpkg) используются, чтобы вырезать из глобального растра только нужный регион: Растр → Извлечение → Обрезать растр по маске. Повторите обрезку отдельно для каждого временного среза (например, 2022, 2023 и 2024 года).

путь к инструменту: Растр → Извлечение → Обрезать растр по маске
Рис. П11.9 — путь к инструменту: Растр → Извлечение → Обрезать растр по маске

В диалоге инструмента укажите исходный слой (например, растр 2023 года) и слой маски — векторную границу интересующей территории:

диалог «Обрезать растр по маске»: исходный слой — растр 2023 года, слой маски — граница Рязанской области
Рис. П11.10 — диалог «Обрезать растр по маске»: исходный слой — растр 2023 года, слой маски — граница Рязанской области

Типичная ошибка

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

Работа с ЦМР

Помимо самих карт состояний, для модели полезны дополнительные пространственные переменные — например, рельеф местности. Откройте цифровую модель рельефа SRTM.tif.

SRTM (Shuttle Radar Topography Mission) — радиолокационная топографическая миссия шаттла, международный проект по созданию цифровой модели высот Земли методом радарной топографической съёмки поверхности.

Обрежьте ЦМР той же маской и тем же инструментом Извлечение → Обрезать растр по маске, что и карты состояний.

Изменение размера растра

Обрезанная ЦМР и обрезанные карты состояний, как правило, имеют разное пространственное разрешение и охват (SRTM — 30 м, MODIS Land Cover — 500 м), поэтому их нужно привести к одной сетке пикселей. Откройте Растр → Калькулятор растров:

  1. Дважды щёлкните по обрезанной ЦМР в списке каналов, чтобы вставить её в выражение.
  2. Одним щелчком выберите (не вставляя в выражение) любой из обрезанных классифицированных растров состояний.
  3. Нажмите «Использовать охват выделенного слоя» — калькулятор возьмёт пространственный охват и размер ячейки от классифицированного растра.
  4. Укажите, куда сохранить новый (пересчитанный) растр ЦМР.
калькулятор растров: выражение содержит обрезанную ЦМР, охват и размер выходного растра взяты у классифицированного растра состояния
Рис. П11.11 — калькулятор растров: выражение содержит обрезанную ЦМР, охват и размер выходного растра взяты у классифицированного растра состояния

Расчёт крутизны по ЦМР

Крутизна склона — ещё одна пространственная переменная, которая помогает модели учитывать рельеф при прогнозировании изменений. Найдите инструмент «Крутизна» на панели Инструменты анализа (Вид → Панели → Инструменты анализа), в группе Растр → Анализ рельефа:

поиск инструмента «Крутизна» на панели «Инструменты анализа»: Растр → Анализ рельефа → Крутизна
Рис. П11.12 — поиск инструмента «Крутизна» на панели «Инструменты анализа»: Растр → Анализ рельефа → Крутизна

Передайте пересчитанную (приведённую к нужному разрешению) ЦМР во вкладку «Слой рельефа» и запустите расчёт. Результат — растр крутизны склона на всю территорию:

результат расчёта крутизны склона по ЦМР для территории Рязанской области
Рис. П11.13 — результат расчёта крутизны склона по ЦМР для территории Рязанской области

При необходимости повторите обрезку и приведение к общему разрешению для результата — все пространственные переменные, которые пойдут в MOLUSCE, должны совпадать по охвату, разрешению и размеру в пикселях.

Загрузка данных и обучение модели

Откройте плагин через пиктограмму на панели инструментов или через меню Растр → Molusce → Molusce. Интерфейс плагина состоит из нескольких вкладок: Входные данные, Корреляции, Изменение территорий, Моделирование, Симуляция, Валидация, Сообщения.

интерфейс плагина MOLUSCE: вкладка «Входные данные» со списком доступных растровых слоёв и полями для растра начальных/конечных состояний
Рис. П11.14 — интерфейс плагина MOLUSCE: вкладка «Входные данные» со списком доступных растровых слоёв и полями для растра начальных/конечных состояний

На вкладке «Входные данные»:

  1. Выберите растр более раннего года и нажмите «Растр начальных состояний».
  2. Выберите растр более позднего года и нажмите «Растр конечных состояний».
  3. Добавьте рассчитанную крутизну (и другие пространственные переменные) кнопкой «Добавить» в список «Пространственные переменные».
  4. Нажмите «Проверить совпадение пространственных характеристик» — только после успешной проверки («Пространственные характеристики совпадают») становятся доступны остальные вкладки.

Типичная ошибка

Если проверка совпадения пространственных характеристик не проходит, дело почти всегда в том, что один из слоёв (чаще всего — дополнительная переменная вроде крутизны) не был приведён к тому же охвату и разрешению, что растры состояний, — см. раздел «Изменение размера растра» выше.

Изменение территорий

На вкладке «Изменение территорий» нажмите «Обновить таблицы» — будут построены «Статистика по классам» (площадь и доля каждого класса на начальную и конечную дату) и «Матрица переходов» (сколько площади перешло из одного класса в другой):

вкладка «Изменение территорий»: таблица статистики по классам и матрица переходов между классами
Рис. П11.15 — вкладка «Изменение территорий»: таблица статистики по классам и матрица переходов между классами

Эта информация уже сама по себе полезна для анализа динамики территории — до всякого прогнозирования. Полученные таблицы можно скопировать из контекстного меню ячейки (вместе с заголовками строк и столбцов) и вставить в отчёт. Кнопка «Создать карту изменений» сохраняет растр, в котором каждый тип перехода между классами закодирован отдельным цветом.

Моделирование: обучение модели для предсказания будущих изменений

На вкладке «Моделирование» доступны четыре типа модели: нейронные сети (многослойный персептрон — используется чаще всего), значимость признака, многокритериальная оценка и логистическая регрессия.

вкладка «Моделирование»: метод «Нейронные сети (многослойный перцептрон)» и параметры обучения
Рис. П11.16 — вкладка «Моделирование»: метод «Нейронные сети (многослойный перцептрон)» и параметры обучения

Перед обучением задаются:

  • тип выборки — полная (все данные, медленно), случайная (быстрее, но может пропустить редкие переходы) или стратифицированная (стремится равномерно охватить все сценарии перехода), и количество примеров в выборке;
  • количество соседних пикселей, учитываемых при анализе (обычно 0 или 1);
  • коэффициент скорости обучения — чем меньше значение, тем осторожнее (и медленнее) обучается модель;
  • максимальное количество итераций — слишком высокое значение может привести к переобучению;
  • количество скрытых слоёв — уровень сложности модели.

Нажмите «Обучить нейронную сеть». На графике строятся кривая обучения (Train) и кривая ошибки на контрольной выборке (Validation) — по их поведению можно судить, получилась ли модель:

три сценария кривых обучения: плавное совместное снижение (модель хороша), хаотичные пики без улучшения (гипотеза/данные не подходят), расхождение кривых — переобучение
Рис. П11.17 — три сценария кривых обучения: плавное совместное снижение (модель хороша), хаотичные пики без улучшения (гипотеза/данные не подходят), расхождение кривых — переобучение

Типичная ошибка

Если кривая ошибки на контрольной выборке перестаёт снижаться или начинает расти, пока кривая обучения продолжает падать, — это переобучение: модель запомнила обучающую выборку, а не научилась общей закономерности. Решение — уменьшить количество скрытых слоёв или увеличить количество примеров в выборке, а не продолжать обучение дольше.

При успешном обучении обе кривые плавно снижаются, а показатель «Текущая каппа-статистика» составляет около 0,8 и выше.

Симуляция состояния территории на заданное время вперёд

После обучения модель можно применить для построения прогноза. На вкладке «Симуляция» задайте число итераций симуляции (количество периодов, на которые строится прогноз; по умолчанию 1) и путь для сохраняемых файлов, затем нажмите «Пуск»:

вкладка «Симуляция»: пути сохранения карты потенциальных переходов, карты уверенности в прогнозе и результата симуляции, число итераций
Рис. П11.18 — вкладка «Симуляция»: пути сохранения карты потенциальных переходов, карты уверенности в прогнозе и результата симуляции, число итераций

Помимо самого результата симуляции можно дополнительно сгенерировать:

  • карту потенциальных переходов — для каждого возможного перехода отдельный растр со значениями 0–100%, показывающий вероятность этого перехода;
  • карту уверенности в прогнозе — значения 0–100%, где низкие значения означают, что переход, скорее всего, не попал в обучающую выборку.

Валидация результатов

Валидацию можно провести, если есть карта истинных значений для того периода, на который был построен прогноз (например, если прогноз строился на 2024 год по данным 2022–2023, а фактическая карта 2024 года тоже доступна). На вкладке «Валидация» укажите карту истинных значений и карту-результат симуляции, затем нажмите «Рассчитать каппа-статистики»:

вкладка «Валидация»: доля правильных ответов и три вида каппа-статистики (общая, гистограммы, положения), график ошибок при разных пространственных масштабах
Рис. П11.19 — вкладка «Валидация»: доля правильных ответов и три вида каппа-статистики (общая, гистограммы, положения), график ошибок при разных пространственных масштабах

Помимо общего процента правильных ответов, MOLUSCE считает три варианта коэффициента каппа: общую (согласие модели с эталоном в целом), каппа гистограммы (совпадение только по долям классов, без учёта их взаимного положения на карте) и каппа положения (насколько точно совпадает именно расположение классов). При включённой опции «Карта ошибок» результат также содержит растр с тремя типами пикселей:

  • постоянный — класс пиксела не менялся, и прогноз это верно предсказал;
  • пустой — модель верно предсказала изменение;
  • ошибка — прогноз не совпал с истинным значением.

Список использованных источников

  1. Официальная документация плагина MOLUSCE в репозитории плагинов QGIS.
  2. MODIS Land Cover Type (MCD12Q1) — User Guide.
Контрольные вопросы