Первичная обработка данных SRTM в ГИС SAGA

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/saga-srtm-preprocessing.html

Возможная последовательность шагов по подготовке данных SRTM к тематическому анализу в ГИС SAGA.

Глобальная цифровая модель высот Shuttle Radar Topography Mission (далее – ЦМВ SRTM), находящаяся в открытом доступе с 2003 года, — общедоступный набор геоданных, активно применяющийся в прикладных исследованиях различной направленности. Ее популярность обуславливается простотой получения, практически глобальным охватом и соответствием запросам среднемасштабного картографирования. По разным оценкам, детальность рельефа, представленного данными SRTM, в целом соответствует таковой на топографических картах масштабов 1:100 000 — 1:50 000.

Четвертое поколение данных SRTM прошло несколько стадий обработки, позволивших повысить исходное качество. Основной целью этих улучшений было заполнение пробелов, характерных для территорий со сложным рельефом, поверхностей занятых водными объектами и прочих типов местностей, плохо поддающихся радарной съемке (например, пустынь). Для этого применялись несколько интерполяционных алгоритмов, а в роли вспомогательных источников использовались локальные и национальные ЦМР более высокого разрешения.

Однако, чтобы оценить пригодность ЦМВ для целей геоморфометрического анализа, рекомендуется дать ответы на следующие вопросы (Reuter et al., 2008):

• насколько точно представлены неровности поверхности (микро- и мезорельеф)?
• насколько точно представлена «гидрологическая» форма земной поверхности (вогнутые/ выпуклые формы рельефа, эрозия/ аккумуляция, дивергентность/ конвергентность потока воды)?
• насколько точно могут быть определены реальные тальвеги и водоразделы?
• насколько согласованы измерения высотных отметок по всей территории исследования?

Оценив с таких позиций данные SRTM можно сделать вывод, что их практическое применение все еще усложняется наличием погрешностей, связанных с технологией получения, т.к. обработка не была направлена на их устранение. К таким в первую очередь следует отнести искажения связанные с неоднородностью земного покрова (растительность, застройка), высокочастотный шум (флуктуации отраженного сигнала) и ложные впадины – их совокупное влияние искажает представление о реальном рельефе местности и усложняет моделирование процессов перераспределения вещества и энергии. Поэтому прежде чем приступить к анализу данных SRTM, рекомендуется провести их предварительную обработку, направленную на (Reuter et al., 2008):

• удаление грубых ошибок и артефактов;
• улучшение аппроксимации рельефа;
• улучшение аппроксимации гидрологических/ экологических процессов (таких как перераспределение поверхностного стока, солнечной радиации, отложений и т.д.).

Рассмотрим более детально одну из возможных последовательностей шагов первичной подготовки данных SRTM в ГИС SAGA. В качестве примера используем фрагмент данных SRTM 44_03, полученный из каталога CGIAR-CSI для листа топокарты масштаба 1:100 000 M-37-121, предварительно прошедшего процедуру привязки.

Импорт данных

SAGA использует собственный формат растровых данных SAGA Grid – *.sgrd, поэтому данные GeoTIFF для начала нужно импортировать с помощью модуля Import/ Export – GDAL/ OGR => GDAL: Import Raster. В диалоговом окне укажем путь к основному файлу с расширением *.tif и нажмем Okay.

По окончании работы модуля на вкладке Data появится новый элемент – откройте его в карте двойным щелчком мыши по имени srtm_44_03. Обратите внимание на такие свойства растра как охват по широте/ долготе, количество строк и столбцов, значения NoData, минимальные и максимальные отметки высот.

Через контекстное меню слоя сохраните импортированный файл в формате SAGA Grid (*.sgrd).

Обрезка фрагмента

Поскольку область нашего интереса ограничивается одним листом топокарты, который по охвату намного меньше фрагмента ЦМВ SRTM 5°×5°, для удобства дальнейшей работы обрежем его в соответствии с координатами листа. Для этого воспользуемся модулем Grid – Tools => Cutting [interactive]. Для запуска этого интерактивного модуля в диалоговом окне необходимо указать растр, над которым будут производиться манипуляции.

После активации модуля в окне сообщений появится уведомление Interactive module execution has been started. Для начала ввода координат инструментом – Action щелкните в любой точке карты. В появившемся окне можно ввести координаты прямоугольника, охватывающего область интереса, – в нашем случае введем координаты углов листа топокарты в десятичных градусах.

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

Нажав Okay, вы увидите, что на вкладке Data появился новый растр – обратите внимание, как отличаются его координаты, число рядов и колонок от исходного.

Теперь остановите работу модуля, убрав галочку в пункте меню Modules => Cutting [interactive]. Согласившись с окончанием работы модуля, вы получите уведомление Interactive module execution has been stopped.

Сохраните новый растр через контекстное меню слоя Save As… под удобным именем, например srtm_m-37-121_gcs.sgrd. Исходный фрагмент теперь можно закрыть, воспользовавшись контекстным меню Close.

Двойным щелчком откроем новый растр: чтобы отрегулировать цветовую шкалу изображения в соответствии с диапазоном значений слоя выберите из контекстного меню пункт Classification => Set Range to Minimum/ Maximum.

Перепроецирование и ресэмплинг

Данные SRTM распространяются в географической СК на основе эллипсоида WGS-84, поэтому для дальнейшего анализа их необходимо перевести в спроецированную СК. Для этого воспользуемся уже знакомым по привязке топокарт модулем Projection – Proj.4 => Proj.4 (Dialog, Grid). В его диалоговом окне сначала введем параметры исходной проекции Source Projection Parameters: выберем географическую СК и датум WGS-84, оставив без изменений прочее.

В качестве исходного растра выберем srtm_m-37-121_gcs и перейдем к диалогу параметров результирующей проекции.

В нем стандартными средствами опишем соответствующую зону проекции UTM по аналогии с тем, как это делалось во время привязки листа топокарты.

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

В результате перед вам появится окно, в котором нужно будет ввести номер зоны UTM, а после — окно с параметрами растра в новой СК (охват, пространственное разрешение, соответствующее количество строк и столбцов матрицы). При этом по умолчанию предлагается размер ячейки ≈90 м, как это и заявлено для данных SRTM. Но для дальнейшего анализа такое пространственное разрешение не очень удобно, поэтому воспользуемся диалогом для его изменения.

Один из простых способов определиться с размером ячейки рассмотрен в статье Hengl, 2006. Согласно предложенному правилу, размер ячейки должен быть равен 0,5 мм аналоговой карты в масштабе исследования. Т.е., если в качестве рабочего масштаба мы выберем масштаб топокарты 1:100 000 (кроме того, именно этому масштабу согласно большинству выводов соответствуют данные SRTM), размер ячейки растра составит 50 м. Обратите внимание, что при вводе числа автоматически пересчитываются и другие значения.

После сообщения Module execution succeeded на вкладке Data появится новый элемент, который через контекстное меню слоя Save As… следует сохранить в рабочую папку проекта под именем srtm_m-37-121_utm.sgrd. Исходный растр srtm_m-37-121_gcs теперь можно закрыть через контекстное меню слоя Close.

Для проверки результатов, загрузите рабочий лист топографической карты М-37-121 в проекции UTM и двойным щелчком откройте его в новом окне, в качестве параметра цветового отображения в свойствах объекта установите Type: RGB. В этой же карте откройте слой srtm_m-37-121_utm. Если все сделано верно, то благодаря тому, что лист топокарты и растр ЦМВ имеют общую проекцию и отвечают одной и той же территории, слои наложатся. На вкладке свойств объекта (справа) в разделе Display установите значение Transparency [%]: 50, нажмите Apply – это сделает слой ЦМВ наполовину прозрачным и вы сможете лучше оценить взаимное соответствие топокарты и данных SRTM.

Для удобства дальнейшей работы можно объединить данные в проект, воспользовавшись меню File => Project => Save Project As….

Фильтрация

Для начала с помощью модуля Shapes – Grid => Contour Lines from Grid построим изолинии на основе ЦМВ. В диалоговом окне модуля выберем соответствующую систему координат и растр ЦМВ, а также укажем высоту сечения рельефа – 5 м.

После сообщения Module execution succeeded на вкладке Data появится новый элемент – полилинейный шейп-файл: сохраните его в рабочую папку проекта через контекстное меню слоя Save As… под именем srtm_5m_pln, а после этого двойным щелчком добавьте в окно карты. Для смены цвета изолиний на более привычный в свойствах объекта установите Color: Maroon и нажмите Apply.

Рисунок изолиний хорошо передает общие черты рельефа, даже в сравнении с топографической картой, но он содержит много мелких неровностей. Они как раз и являются шумовой компонентой, речь о которой шла ранее. Прежде чем проводить анализ ЦМВ, шум необходимо устранить с помощью фильтрации. Из-за сложности характера распределения полностью избавиться от шума не удастся, поэтому главная задача фильтрации – максимальная элиминация шумовой компоненты без утраты характерных черт рельефа местности.

Воспользуемся простым однородным фильтром Grid – Filter => Simple Filter. Суть работы фильтра состоит в следующем: он получает новые значения ячеек растра в соответствии с некоторым математическим выражением, т.е. пересчитывает значения центральной ячейки на основе значений ее соседей. Результат фильтрации зависит от параметров скользящего окна, представленных в группе Options пунктами Search Mode, Filter и Radius. Search Mode и Radius совместно контролируют число соседних ячеек растра, которые будут учитываться при расчете нового значения центральной ячейки.

Пункт меню группы Options	Что определяет	Варианты	Что означает
Search Mode	окно поиска – форма матрицы для пересчета значений центральной ячейки	Circle	сферическая матрица
		Square	квадратная матрица
Filter	тип фильтра, который будет применяться для пересчета значений	Smooth	сглаживание – усреднит разницу между центральной ячейкой и ее окружением; новое значение рассчитывается по формуле $z'={\bar {z}}$, где ${\bar {z}}$ – среднее арифметическое значение в окне анализа
		Sharpen	заострение – имеет противоположный по сравнению с предыдущим эффект, поскольку усиливает различия в значениях ячеек; новое значение рассчитывается по формуле $z'=2z-{\bar {z}}$
		Edge	усиление кромок – выделение линий с высокой вариативность значений (например, линий перегиба рельефа); новое значение рассчитывается по формуле $z'=z-{\bar {z}}$
Radius	размер матрицы, которая будет использована для пересчета значений	количество ячеек	чем больше значение – тем сильнее выражен эффект выбранного фильтра

Установим параметры фильтрации в окне модуля следующим образом:

После завершения работы модуля и появления сообщения Module execution succeeded на вкладке Data появится новый элемент. Через контекстное меню слоя Save As… сохраните его в рабочую папку проекта под именем srtm_simple_fltr.sgrd и двойным щелчком откройте в окно новой карты.

Для оценки полученного результата уже известным способом с помощью модуля Shapes – Grid => Contour Lines from Grid построим изолинии для растра srtm_simple_fltr. Через меню Save As… сохраним векторный слой изолиний в рабочую папку проекта под именем srtm_simple_fltr_5m_pln.shp и двойным щелчком откроем в окно карты с профильтрованной ЦМВ. Чтобы визуальное сопоставление результатов было более удобным, поместим карты рядом с помощью пункта меню Window – Tile Vertically и синхронизируем их с помощью инструмента панели меню Synchronise Map Extents.

В результате мы можем изменяя масштаб отображения и перемещая карту в одном окне, получать идентичное изображение в другом.

Данные SRTM до и после обработки простым фильтром

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

Кроме того, мы использовали самый простой вариант фильтра, но библиотеки SAGA Grid – Filter и Grid – Filter (Perego 2009) содержат 20 различных модулей фильтров, которые могут использоваться как по отдельности, так и совместно в зависимости от особенностей распределения шумовой компоненты, типов артефактов, характера рельефа и т.д.

Гидрологическая коррекция

Модули гидрологической коррекции объединены библиотекой Terrain Analysis – Preprocessing, которая содержит различные инструменты для заполнения локальных впадин (sinks) для облегчения дальнейшего гидрологического анализа. Воспользуемся модулем Terrain Analysis – Preprocessing => Fill Sinks (Planchon/ Darboux, 2001). Принцип его действия состоит в следующем: вместо постепенного заполнения локальных понижений, вся поверхность сначала «заливается» слоем воды, а после ее излишек удаляется, оставляя после себя заполненные понижения. Гибкости алгоритму дает возможность заполнять впадины как до горизонтальной поверхности, так и с сохранением незначительного уклона (например, 0,01°). Первый вариант удобен в том случае, если нужно оценить объем впадины, второй – для оконтуривания дренажной сети. Если в дальнейшем планируется проводить гидрологический анализ ЦМР, то следует выбрать второй вариант.

После завершения работы модуля и сообщения Module execution succeeded во вкладке Data появится новый элемент srtm_simple_fltr [no sinks], который следует через контекстное меню слоя Save As… следует сохранить в рабочую папку проекта под именем srtm_nosinks.sgrd. Двойным щелчком по имени слоя откройте его в новую карту.

Для оценки полученного результата уже известным способом с помощью модуля Shapes – Grid => Contour Lines from Grid построим изолинии для растра srtm_nosinks. Через меню Save As… сохраним векторный слой изолиний в рабочую папку проекта под именем srtm_nosinks_5m_pln.shp, после чего двойным щелчком откроем в окно карты с гидрологически откорректированной ЦМВ.

Еще один способ оценить и визуализировать результаты гидрологической коррекции – провести количественную оценку отличий между ЦМВ до и после нее. Для этого воспользуемся калькулятором растров Grid – Calculus => Grid Calculator. Для начала в его диалоговом окне из выпадающего списка Grid System выберем систему координат растров для которых будут проводиться расчеты. После этого в поле >>Grids нужно указать те растровые слои, которые войдут в формулу (в нашем случае srtm_simple_fltr и srtm_nosinks).

Вернувшись в основное окно модуля, вводим формулу. Чтобы узнать, где именно расположены заполненные впадины, следует просто отнять от srtm_nosinks (растр g2) слой srtm_simple_fltr (растр g1), т.е. формула будет выглядеть как g2–g1. Для присвоения новому растру удобного содержательного имени (а не формулы, как предлагается по умолчанию), вводим название filled_sinks в поле Name и убираем галочку в поле Take Formula.

После завершения работы модуля и сообщения Module execution succeeded во вкладке Data появится новый элемент filled_sinks, который через контекстное меню слоя Save As… следует сохранить в рабочую папку проекта под именем filled_sinks.sgrd. Двойным щелчком по имени слоя откройте его в новую карту.

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

Сравним между собой результаты каждого этапа. Перейдем на вкладку Data, для каждого из четырех растров ЦМВ на вкладке Settings справа активизируем параметр Show Cell Values, нажмем Apply – в результате при увеличении фрагмента будут отображаться значения в каждом пикселе и нам будет легче оценить насколько в процессе обработки изменились данные.

Разместим все карты в рабочем окне через меню Window – Tile Vertically и синхронизируем экстенты с помощью инструмента – Synchronise Map Extents. Это позволит сопоставить результаты как визуально

так и количественно с отображением значений пикселей при увеличении отдельных участков

Теперь наши данные готовы к дальнейшему полноценному морфометрическому и гидрологическому анализу.

Ссылки по теме

Описание и получение данных SRTM
Открытая настольная ГИС SAGA - общая характеристика
Olaya, V. A gentle introduction to SAGA GIS. 2004.
Cimmery, V. User Guide for SAGA. Vol. 1,2. 2010