Анализ рельефа
1 Цифровая модель рельефа
Цифровая модель рельефа (ЦМР) — это растровое представление непрерывной поверхности, обычно относящееся к поверхности Земли1.
Различают цифровую модель рельефа и цифровую модель местности. Разница между ними состоит в том, что цифровая модель местности включает в себя все объекты, расположенные на поверхности, тогда как цифровая модель рельефа - это только поверхность земли без растительности и строений.
Цифровые модели рельефа, как правило, строят либо по результатам интерполяции, либо на основе открытых данных.
Самым распространенным источником данных о рельефе является цифровая модель рельефа SRTM (shuttle radar topographic mission). Эта цифровая модель была получена в 2000 году на основе спутниковой радарной съемки. Она охватывает планету между 54 градусами южной широты и 60 градусами северной широты. Подробнее на русском можно почитать про нее здесь
Скачать данные SRTM можно на сайте Геологической службы США (https://earthexplorer.usgs.gov/) или с помощью плагина SRTM downloader2.
Для использования модуля и загрузки данных напрямую из QGIS все равно необходимы реквизиты аккаунта на сайте USGS.
Значительным недостатком SRTM является ее ограниченный пространственный охват, из-за которого значительная часть нашей страны не покрыта этой цифровой моделью рельефа.
Для примера и выполнения работы были скачаны данные для окрестностей Эльбруса (файл dem_elbrus.tiff).
Из SRTM выгружается сразу целая 1-градусная трапеция (фрагмент 1 градус по широте и 1 градус по долготе), в этом случае трапеция была обрезана до более маленького фрагмента (вообще не очень рекомендую работать с большими растрами, охватывающими сразу градус и более по долготе и широте).
Более современной и точной цифровой моделью рельефа является модель FABDEM, к тому же она действительно глобальная и охватывает всю территорию земного шара.
Немного почитать о ней на русском можно по ссылке.
Самое главное ее отличие от более ранних моделей в наличии процессинга и удалении из нее объектов, расположенных на поверхности Земли.
Сама модель является открытой и данные можно скачать по ссылке из репозитория - https://data.bris.ac.uk/data/dataset/s5hqmjcdj8yo2ibzi9b4ew3sn
Все данные организованы в архивы, которые включают в себя охват 10 широты на 10 градусов долготы. Растровые данные внутри архивов - это одноградусные трапеции (как и в SRTM).
Для того, чтобы понять, какой именно архив и файл из него нужен вам на исследуемую территорию, можно скачать файл FABDEM_v1-2_tiles.geojson и открыть его либо в QGIS, либо с использованием сервиса geojson.io. Прямоугольники в этом файле показывают одноградусные трапеции, а в их атрибутах указано название архива и файла.
В каждой ячейке исходного растрового изображения содержится значение высоты.
Добавление растрового слоя осуществляется через строку меню Слой \(\longrightarrow\) Добавить слой \(\longrightarrow\) Добавить растровый слой.
После открытия файла его лучше сразу перепроецировать на лету в EPSG:3857.
Перепроецирование на лету не будет изменять систему координат слоя, а только изменит систему координат проекта и то, как ваша карта будет отображаться на экране.
Перепроецирование на лету может быть выполнено с использованием кнопки 
в правом нижнем углу основного окна программы.
В результате будет открыто изображение рельефа в оттенках серого цвета.
Это изображение можно сделать цветным, для этого нужно открыть свойства слоя. В настройках стиля нужно выбрать тип изображения Одноканальное псевдоцветное.
Одноканальное здесь очевидно потому, что изображение не поделено на каналы и существует только в одном.
Градиент можно также создать самостоятельно с помощью редактора градиентов.
Или выбрать из готового каталога градиентов. Для этого при создании нового градиента нужно выбрать каталог cpt-city.
Далее в каталоге слева выбрать тип градиента Topography и подобрать нужный градиент (можно взять elevation или с3t3).
2 Теневой рельеф3
2.1 Расчет теневого рельефа
Для более наглядного отображения рельефа на карте рассчитывается теневой рельеф, который строится с использованием позиции источника направленного света, или теневая отмывка.
В данном случае одним из самых важных параметров является Масштабирование по Z - этот параметр позволяет соотнести между собой единицы измерения плановых координат и высот. В нашем случае высоты в ячейках растра приведены в метрах, а координаты - в градусах., кроме того, стоит учитывать, что в географических системах координат искажение зависит от долготы.
С этим можно справиться двумя способами - перепроецировать растра в прямоугольную систему координат или задать масштабирование по Z (Z-factor). Фактически с помощью этого коэффициента мы переходим от градусов к метрам. Так как у нас окрестности Эльбруса, то можно воспользоваться значением коэффициента для 40 градусов долготы - 0.00001171.4
Азимут и вертикальный угол задают положение источника света.
Этот теневой рельеф можно использовать, чтобы добавить дополнительный объем подложке (предварительно нужно добавить) или исходному рельефу растра. Для этого нужно открыть свойства слоя с теневым рельефом и изменить режим смешивания на Добавление, Осветление или Затемнение (также можно поэкспериментировать с настройками яркости и контрастности).
Также можно воспользоваться функцией Цветной рельеф, которая позволяет создавать классы автоматически и задавать им цвет.
С использованием такого растра может быть построена трехмерная поверхность, отображающая рельеф местности, а также выполнен ряд вычислений для морфометрического анализа рельефа.
Все растровые операции здесь осуществляются за счет разности значений между ячейками.
2.2 Настройка теневого рельефа как стиля растрового слоя
Рассчитывать теневой рельеф и сохранять его в отдельный слой не обязательно, вы можете просто выбрать его как тип стиля для растрового слоя в настройках свойств.
Полученный теневой рельеф после настройки будет отображен на карте.
При создании теневого рельефа с использованием настроек стиля растрового слоя на карте могут появляться артефакты - небольшая “рябь” в картинке.
2.3 Использование модуля Terrain shading
Более художественных результатов при создании теневого рельефа позволяет добиться модуль Terrain shading.
Чуть подробнее про его функции можно почитать в блоге разработчика.
Установка модуля осуществляется стандартно (Модули \(\longrightarrow\) Управление модулями). После установки модуль должен появиться у вас в панели Инструменты анализа.
Модуль содержит сразу три функции, которые позволяют строить теневую отмывку:
- Hillshade (terrain shading) - аналогичен стандартному построению теневого рельефа;
- Shadow depth - инструмент для генерации теней;
- Toposhade - комбинация теневого рельефа и индекса топографического положения (о нем будет еще чуть дальше), который показывает выпуклость и вогнутость рельефа. Такая комбинация позволяет добиться более выраженной карты теневого рельефа.
Полученные новые растры можно комбинировать друг с другом и другими слоями для получения более красивого результата.
3 Расчет параметров рельефа
3.1 Крутизна склонов5
Расчет крутизны склонов или их уклона и экспозиции осуществляется на основе производных первого порядка.
3.1.1 Расчет крутизны склонов
Рассчитаем крутизну склонов - угол наклона каждого склона в градусах.
3.1.2 Выделение территорий пригодных для строительства с учетом уклона
Для того, чтобы выделить области пригодные для строительства, воспользуемся переклассификацией растра.
При переклассификации будем рассматривать территории с уклоном менее 10 %.
Так как уклон был рассчитан в градусах, а параметр пригодности для строительства задан в процентах, то пересчитаем этот параметр в градусы по формуле:
\[ Slope_{degree}=arctg(\dfrac{Slope_{percent}}{100}) \]
Полученное значение пригодности - 6 градусов.
Переклассифицируем растр и выделим области пригодные или непригодные для строительства.
Обратите внимание, что здесь в дополнительных параметрах можно указать, как именно будут задаваться границы интервалов, то есть будут ли включаться в них нижний или верхний предел интервала.
В таблице реклассификации следует ввести границы интервалов и соответствующее им значение (эти значения должны быть целочисленными.
В результате будет получен дискретный растр с двумя значениями (или с тем количеством значений классов, которое вы указали в таблице реклассификации.
3.2 Экспозиция склонов6
Экспозиция поверхности – угол по часовой стрелке между определенным направлением (как правило, на север) и проекцией уклона на горизонтальную плоскость; фиксирует направление (азимут) максимального уклона (градиента) земной поверхности7.
Эта функция позволяет определить ориентацию склонов по сторонам света.
3.3 Расчет неоднородности рельефа8
Также по растру с данными о рельефе может быть рассчитана неоднородность рельефа.
Этот показатель рассчитывается на основе разности значений между пикселями в сетке 3 на 3. Таким образом, каждый полученный пиксель содержит разницу между центральным пикселем и 8 его окружащими.
4 Построение трехмерной поверхности
Также QGIS поддреживает отображение трехмерных поверхностей. Для этого нужно просто в строке меню выбрать Вид - Новая 3D карта.
После чего в окне программы появится окно трехмерной карты (перемещение по карте левой кнопкой мыши, поворот и наклон - колесико, приближение - правая кнопка мыши).
Сразу после открытия нового окна поверхность все еще является плоской, потому что не задано на основе чего строить трехмерную поверхность.
В окне трехмерной карты нужно нажать кнопку Настройки (в панели значков крайняя справа) и в открывшемся окне настроек выбрать тип DEM (raster layer), высота (растровый слой с высотами) и вертикальный масштаб (здесь лучше всзять что-то от 1 до 5, чтобы поверхность была повыразительнее).
Дополнительно можно добавить Terrain shading (фактически теневая отмывка) и подкорректировать tile resolution (чем больше, тем более подробным будет рельеф), и skirt height (если правильно нашла - высота подставки под саму поверхность).
После чего в окне будет отображена трехмерная поверхность.
5 Создание горизонталей
Из растра со значениями высот могут быть извлечены горизонтали
Расстояние между изолиниями - это фактически высота сечения рельефа, то есть через сколько мы будем проводить горизонтали.
Полученные изолинии
Настроим горизонтали таким образом, чтобы каждая пятая была утолщенной, добавим подписи и настроим их так, чтобы они были “головой вверх”.
Начнем с настройки толщины линий. Для этого в свойствах слоя выберем настройку толщины в зависимости от выражения.
Далее пропишем логическое выражение: if(“ELEV”%500=0,0.7,0.3). В этом выражении перед запятой прописано условие - высота горизонтали делится на 500 без остатка, потом толщина линии, если условие выполняется, и толщина, если условие не выполняется.
Результат
Добавим подписи.
Но при таком добавлении подписей, будут подписаны все горизонтали, а мы хотим подписать только утолщенные. Для этого нужно прописать выражение, подобное тому, что мы писали выше: if(“ELEV”%500=0,“ELEV”,““).
Добавим белую обводку подписям, чтобы они не сливались с горизонталями
Сделаем так, чтобы они размещались на линии и изгибались в соответствии с ней.
И настроим подписи так, чтобы все они смотрели головой в сторону увеличения высоты.
В результате получим горизонтали, где каждая пятая будет утолщена и подписана головой в сторону увеличения.
6 Построение профиля местности
В актуальных версиях программы построения профиля местности по цифровой модели рельефа является функцией по умолчанию. Этот инструмент можно открыть из строки меню Вид \(\longrightarrow\) Профиль высот.
Окно профиля высот открывается в нижней части интерфейса. Ваша цифровая модель рельефа (если растр с ней загружен в проект должна автоматически добавиться в левую часть окна построения профиля.
Если ваша цифровая модель рельефа не добавилась в высотный профиль автоматически, то вам нужно проверить ее настройки.
В свойствах слоя (открываются двойным кликом или кликом правой кнопки мыши и выбором пункта Свойства из открывшегося контекстного меню) необходимо найти пункт Высота и поставить галочку напротив настройки “Представляет собой данные высоты”.
После этого ваша цифровая модель будет автоматически добавлена в окно высотного профиля.
Построение профиля может происходить двумя способами:
Захват кривой
- отрисовка линии профиля осуществляется вручную прямо на карте;Захват кривой из объекта
- использование уже существующей линии из векторного слоя.
В первом случае необходимо просто нарисовать нужную вам линию профиля (она может быть отрисована как любая ломаная линия), по которой будет происходить построения.
При перемещении курсора по графику или линии профиля вы увидите перемещающуюся точку, показывающую конкретное положение.
Функция прилипания 
как и при обычной векторизации позволяет использовать другие объекты в качестве исходных.
Полученный профиль вы можете сохранить в виде картинки графика или файла pdf 
или же результаты могут быть сохранены в виде линии профиля (в двумерном и трехмерном видах) для последующего использования в ГИС или CAD программах и в виде таблицы значений высот вдоль профиля.
7 Выявление форм рельефа
Для нижеописанных расчетов используются инструменты программы SAGA. Как и QGIS это бесплатное программное обеспечение.
Ее вы можете установить себе на компьютер как отдельное приложение, скачав, например, по ссылке. Возможно она у вас уже установлена вместе с QGIS (если вы устанавливали его пакетно).
Либо вы можете воспользоваться этими инструментами напрямую из QGIS. Для этого необходимо установить модуль Processing Saga NextGen Provider, после чего скачать программу для своей операционной системы по ссылке и разместить на своем компьютере в папке в надежном месте. Далее в настройках программы QGIS вам нужно прописать, где же хранится SAGA.10
Установку для Windows можно посмотреть в видео по ссылке.
7.1 Расчет структурных линий рельефа
Выделение структурны линий рельефа можно осуществить с использованием инструмента Surface specific points в наборе инструментов SAGA.
Чем большее значение порога вы установите, тем более крупные структурные линии будут выделяться.
7.2 Определение выпуклых и вогнутых форм рельеф
Формы рельефа определяются на основе индекса топографического положения (Topographic position index - TPI). Этот индекс позволяет классифицировать формы рельефа по высоте конкретного пикселя относительно соседей: если он ниже, то значение индекса будет отрицательным и он будет отнесен к вогнутой форме рельефа; если он выше соседей, то индекс будет больше нуля и этот пиксель будет классифицирован как выпуклая часть рельефа.
7.2.1 Расчет стандартными инструментами
Индекс топографического положения может быть рассчитан с использованием стандартного инструмента из строки меню: Растр \(\longrightarrow\) Анализ \(\longrightarrow\) Индекс топографического положения (TPI).
Здесь достаточно указать растр, используемый для расчета и необходимость обработки краев: так как расчеты осуществляются с использованием движущегося окна вокруг каждого из пикселей, то края обычно не обрабатываются, так как эти пиксели не имеют достаточного числа соседей.
Полученный результат вы увидите на растре (на картинке ниже взята цветная шкала, по умолчанию вы получите изображение в оттенках серого).
В примере превалируют либо ровные области, либо небольшие выпуклости и вогнутости рельефа.
Однако такой результат может быть еще связан с небольшим окном для расчетов. Другие варианты расчета дают большую гибкость и позволяют четче выделять формы рельефа.
7.2.2 Расчет с использованием функций SAGA
В наборе инструментов SAGA есть целая группа, предназначенная специально для морфометрического анализа рельефа - Terrain analysis - Morphometry, среди которых в том числе есть и расчет индекса топографического положения - Topographic position index (tpi).
Использование этого инструмента в отличие от стандартного дает больше контроля пользователю, так как здесь можно выбирать размер ячеек для анализа, а также задавать весовые значения.
За счет большей адаптивности этого инструмента получается более легко и наглядно интерпретируемый результат.
Также есть инструмент для классификации форм рельефа на основе индекса топографического положения - Tpi based land form classification.
Но, к сожелению, у меня расчет этим инструментов так и не окончился успехом.
7.2.3 Расчет с использованием модуля Terrain shading
Для расчетов топографического индекса положения необходимо, чтобы данные были в прямоугольной системе координат, поэтому растр с цифровой моделью рельефа необходимо перепроецировать: Растр \(\longrightarrow\) Проекции \(\longrightarrow\) Деформация (перепроецирование).
Наиболее оптимальным является перепроецирование в плоскую прямоугольную систему координат из группы систем координат UTM. Для расчета нужной вам зоны системы координат можно воспользоваться сервисом по ссылке https://www.latlong.net/lat-long-utm.html
Для расчета индекса топографического положения воспользуемся инструментом Topographic position (TPI).
Так же как и для инструментов SAGA вы можете здесь настроить возможность использования весов и радиус обработки.
Чем больший радиус вы установите, тем более крупные формы рельефа будут обнаруживаться.
Недостатком полученного результата будет являться добавление дополнительных пикселей вдоль границ исходного растра со значениями, лежащими за пределами корректного диапазона индекса.
Однако это можно исправить настройкой стиля слоя, в которых можно указать минимальные и максимальные значения растра и исключить все, что не попадает в пределы диапазона.
Полученный результат хоть и будет содержать небольшую рамку, но уже не будет так сильно искажать интерпретацию.
Также полученный растр можно обрезать по маске.
Как показали несколько попыток расчета можно попробовать осуществить его для растра в географической системе координат, так вы не получите дополнительных полей по краям.
Дополнительно вы можете наложить на свой растр с индексом топографического положения горизонтали, чтобы понять, действительно ли рельеф изменяется согласно полученным значениям.
8 Классификация форм рельефа
Алгоритм классификации форм рельефа в наборе инструментов SAGA основан на статье Iwahashi, J. & Pike, R.J. (2007): Automated classifications of topography from DEMs by an unsupervised nested-means algorithm and a three-part geometric signature.
В соответствии с предложенными авторами алгоритмом вся территория автоматически делится на категории по форме поверхности с использованием - градиента склона, локальной выпуклости и текстуре поверхности.
9 Определение открытых и закрытых элементов рельефа
Для расчета открытых и закрытых форм рельефа в комплекте инструментов SAGA есть Topographic openness, которая позволяет рассчитать два растра: с открытыми формами и закрытыми соответственно.
10 Анализ видимости
Анализ видимости здесь рассмотрен на основе модуля Visibility analysis.
10.1 Индекс видимости11
Индекс видимости - это метрика, показывающая поле зрения для любой точки ланшафта.
Индекс видимости рассчитывается как отношение положительных визуальных связей: 1,0 или 100% означает, что точка видна со всех своих соседей.
В качестве исходных данных здесь задается только цифровая модель рельефа, остальные параметры устанавливаются непосредственно при расчете:
радиус анализа - насколько далеко необходимо осуществлять видимость;
высота наблюдателя, то есть высота, с которой осуществляется обзор относительно поверхности земли;
sample - количество линий обзора;
направление - расчет с точки зрения наблюдения за ландшафтом или как точек обзора.
10.2 Расчет бассейнов видимости12
Результатом расчета бассейна видимости, как правило, является бинарный растр со значениями 1 (пиксель видим с точки обзора) или 0 (пиксель не виден с точки обзора).
Если точек обзора много, то их бассейны видимости будут суммироваться при наложении.
На первом этапе необходимо создать точки обзора из имеющегося точечного слоя - Create viewpoints.
В качестве исходных данных в расчет передаются точечный слой, для которого будут осуществляться расчеты, и цифровая модель рельефа.
Важно, что слой с точками обзора и цифровая модель рельефа должны быть в одной системе координат.
Обязательными параметрами являются:
радиус анализа (в метрах) - максимальное расстояние видимости;
высота наблюдателя, то есть высота, с которой ведется обзор (здесь может быть задана как высота среднего роста человека, так и высота здания, если обзор будет осуществляться с него);
целевая высота - значение высоты, которое добавляется ко всем участкам местности, проверяемым на видимость из точки наблюдателя.
В результате вы получите точечный слой из объектов, которые могут использоваться как точки обзора, так как имеют свою собственную высоту.
Далее можно рассчитать бассейны видимости с помощью инструмента Viewshed.
Инструмент позволяет рассчитать три вида бассейнов видимости:
бинарный - видно/не видно (1/0);
расположение ниже горизонта - высота, на которую должно подняться каждое место, чтобы стать видимым;
горизонт - внешний край бассейна видимости.
При использовании двух точек обзора были получены зоны видимости, показанные на рисунке ниже.
Для отображения зон видимости поверх остального рельефа можно настроить стиль полученного слоя как Палитра/Уникальные значения и сделать нулевые значения прозрачными.
В результате на карте у вас будут отображаться только те зоны, которые видны с ваших точек обзора, поверх остальных слоев.
Для построения зон видимости в условиях городской застройки можно воспользоваться библиотекой для python - ObjectNat.
В том случае, если у вас есть цифровая модель поверхности со всеми зданиями для исследуемой территории, можно воспользоваться модулем Visibility analysis и вот этим туториалом.
Для добавления зданий к цифровой модели рельефа можно использовать модуль UMEP.
11 Расчет высоты расположения зданий
В том случае, если необходимо рассчитать значение высот в пределах объектов векторного слоя (административных районов, например) или значение высот, на которых расположены здания, можно воспользоваться стандартным инструментом Зональная статистика.
Этот инструмент позволяет рассчитать значение какой-либо из статистик растрового слоя в пределах векторных объектов.
Этот же инструмент может использоваться, например, для вычисления средних значений вегетационного индекса в пределах кварталов.
В качестве исходных данных была использована цифровая модель рельефа FABDEM для территории Москвы и векторный слой со зданиями, скачанными из OSM.
Доступные для расчета статистические показатели приведены на рисунке ниже. Вы можете использовать один из них или любой набор из нескольких.
По результатам расчета вы получите векторный слой с объектами из исходного слоя, но в их атрибутивную таблицу будет добавлено новое поле с полученными значениями показателя (или несколько полей, если вы рассчитывали несколько показателей).
Footnotes
Что такое цифровые модели рельефа (ЦМР) - https://pro.arcgis.com/ru/pro-app/latest/tool-reference/spatial-analyst/exploring-digital-elevation-models.htm↩︎
Плагин SRTM-Downloader — загрузка SRTM с сервера прямо в QGIS https://cartetika.ru/tpost/98m7ngx971-plagin-srtm-downloader-zagruzka-srtm-s-s↩︎
https://docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/rasterterrainanalysis.html?highlight=relief%20analysis#hillshade↩︎
https://www.esri.com/arcgis-blog/products/product/imagery/setting-the-z-factor-parameter-correctly/↩︎
Документация функции https://docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/rasterterrainanalysis.html?highlight=relief%20analysis#slope↩︎
Экспозиция - https://docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/rasterterrainanalysis.html?highlight=relief%20analysis#aspect↩︎
Основные геоморфометрические параметры: теория - https://gis-lab.info/qa/geomorphometric-parameters-theory.html#%D0%9E%D1%81%D0%BD%D0%BE%D0%B2%D0%BD%D1%8B%D0%B5_%D0%B3%D0%B5%D0%BE%D0%BC%D0%BE%D1%80%D1%84%D0%BE%D0%BC%D0%B5%D1%82%D1%80%D0%B8%D1%87%D0%B5%D1%81%D0%BA%D0%B8%D0%B5_%D0%BF%D0%B0%D1%80%D0%B0%D0%BC%D0%B5%D1%82%D1%80%D1%8B↩︎
https://docs.qgis.org/3.16/en/docs/user_manual/processing_algs/qgis/rasterterrainanalysis.html?highlight=relief%20analysis#ruggedness-index↩︎
https://download.osgeo.org/qgis/doc/reference-docs/Terrain_Ruggedness_Index.pdf↩︎
https://gis.stackexchange.com/questions/468320/installing-saga-next-gen-binaries-in-qgis-on-a-mac↩︎
Visibility index (total viewshed) for QGIS : finally there! - https://landscapearchaeology.org/2020/visibility-index/↩︎
Viewshed analysis in QGIS 3: a tutorial - https://landscapearchaeology.org/2020/viewshed-tutorial/↩︎