Статьи

Интерполяция в QGIS, часть 2

ГИС Туториалы
Время прочтения: 6 минут

Введение

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

О данных

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

Подготовка

Кригинг не входит в основной набор функций QGIS, поэтому нам нужно подготовиться, чтобы произвести вычисления.
Во-первых, для кригинга потребуется плагин Smart-map, который вы можете установить так же, как и большинство модулей для QGIS: заходим во вкладку Модули, затем Управление модулями, ищем плагин по названию и устанавливаем.
Во-вторых, может так выйти, что недостаточно просто установить плагин, чтобы все работало. Все плагины QGIS используют какой-то код, написанный на каком-либо языке программирования. Smart-Map использует для кригинга библиотеку Scikit-learn, которую она устанавливает для работы плагина. Но бывает так, что по какой-то причине у нее не получается это сделать — это может быть либо из-за версии QGIS, либо из-за других установленных плагинов, либо из-за того, что вы заклеили пластырем вашу веб-камеру (или у вас нет веб-камеры), и при установке плагина ему не хватило силы вашей улыбки (это, между прочим, главная составляющая профилактики любых багов и глюков open-source программного обеспечения).
В общем, лучше всего превентивно провести мероприятие, которое позволит плагину корректно работать. Для этого вручную установим библиотеку Scikit-learn для нашего QGIS.
Чтобы сделать это, нам нужно открыть приложение OSGeo4W Shell — его можно найти по поиску на компьютере либо в папке с кучей ярлыков на рабочем столе, которая появляется при установке QGIS.
Вот так она выглядит
В открывшейся панели пишем вот такой код и нажимаем enter:
python -m pip install scikit-learn
Если все произошло успешно, то вам напишут текст, примерно как у меня.

А вот и кригинг!

Демонстрировать кригинг будем на том же слое из 301 точки, что и использовали в предыдущем уроке. И теперь важная часть о кригинге в QGIS: использовать большие массивы данных в нем не получится =(
Если ваш объем данных слишком большой, то вы получите либо ошибку на одном из этапов, которая приведет к зависанию программы, либо будете бесконечно смотреть на надпись типа “Progress 0%”. Я связываю это с тем, что для многих действий QGIS не может использовать мощности более одного ядра процессора вашего компьютера (одна из причин, почему некоторые операции с помощью SQL или Python могут стать привлекательнее). У нас, кстати, есть курсы и по первому и по второму. Есть даже вариант для тех, кто уже что-то знает про питон, как я! Например, то что у них красивая чешуя, и они ползают в тропиках, и… наверное… простите, вернемся к уроку.
Исходный датасет точек с высотами
Теперь приступаем к кригингу: на панели вашего QGIS ищем иконку плагина Smart Map и смело жмем на нее.
В правом верхнем углу можно поменять язык на бразильский (португальский). В поле Input Layer выбираем слой с точками и не забываем, чтобы он был в метрической системе координат для корректности вычислений. В поле Z выбираем поле с нужной нам переменной, после чего смело нажимаем Import.
После импорта точек они появятся в таблице ниже. Красным выделены точки, которые по какой-то причине не смогли импортироваться, но потеря нескольких из них не критична, так что можно идти дальше.
После того, как точки импортированы переходим на вторую слева вкладку “Grid”.
Нам выдается простенькая картинка наших точек со значением. Тут самое важное — выставить размер пикселя по X и по Y, сделать это можно в полях слева сверху. Справа программа вам автоматически посчитает количество колонок и линий пикселей. Я выставил точность 100 метров. Будьте внимательны, по умолчанию может стоять очень маленькая точность, из-за которой программа попытается вам выдать растр с разрешением тысячи на тысячи пикселей и, скорее всего, QGIS этого не выдержит.
Вы можете ограничить область интерполяции сами, использовав отдельный полигональный слой — если он есть в проекте, то панель Outline Polygon будет доступна для использования, по умолчанию интерполяция происходит по границам ограничивающего прямоугольника (bounding box).
Помните, что если границы вашего полигона будут далеко от ваших точек, то тогда на определенном удалении от точек начнется экстраполяция, которая будет выражена в повторении среднего значения слоя с вашими точками.
Теперь переходим к вкладке интерполяции (“Interpolation”).
В начале все просто, нажимаем единственную доступную кнопку “Calculate”, чтобы появилась возможность настраивать наши вариограммы.
Теперь у нас появилось много возможностей, о них по порядку. В центре у нас находится блок “Model Adjust”, через нее мы выбираем настройки, чтобы линия вариограммы (это, грубо говоря, смоделированная модель распределения значений в пространстве) максимально пересекала наши точки. Можно выбирать различные модели и настраивать переменные, по которым наша линия будет менять свою форму или расположение относительно осей вариограммы.
На панели справа с названием “Kriagem” можно выбрать радиус и количество соседей от каждой точки, используемых для построения вариограммы, а также нажать на кнопку “Interpolate”, чтобы построить растр интерполированных значений, но перед этим, возможно, есть смысл еще раз проверить модель и провести кросс-валидацию.
Кросс-валидация делается на второй вкладке на панели слева снизу. Она показывает, какие бы значения имели точки с известными значениями, если бы были смоделированы по существующей вариограмме. У меня пару раз было так, что без кросс-валидации QGIS зависал на моменте построения растра, поэтому я просто уже машинально делаю это как ритуал в любом случае.
Нажимаем “Interpolate” и получаем наш растр интерполированных значений, он отобразится как в окне модуля, так и в проекте.

Сравнение методов интерполяции в QGIS

И сейчас я вам предлагаю сравнить результаты интерполяции кригингом, методом обратных взвешенных расстояний и триангуляцией Делоне:
Интерполяция кригингом
Интерполяция методом обратных взвешенных расстояний
Интерполяция методом триангуляции Делоне
Оригинальный растр
А теперь, как говорится, все выводы вы можете сделать сами. Если вам была интересна эта серия статей и вы прочитали их полностью, то ставьте эмодзи с ноготочками к посту в телеграм 💅
Материал подготовил Александр Зуев