GIS-LAB

Географические информационные системы и дистанционное зондирование

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

Обсудить в форуме Комментариев — 37Редактировать в вики

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


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

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

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

Однако, чтобы оценить пригодность ЦМВ для целей геоморфометрического анализа, рекомендуется дать ответы на следующие вопросы (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.

Saga dem prep 01.png

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

Saga dem prep 02.png

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

Saga dem prep 03.png
Saga dem prep 04.png

Обрезка фрагмента[править]

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

Saga dem prep 05.png

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

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

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

Saga dem prep 06.png
Saga dem prep 07.png

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

Saga dem prep 08.png

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

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

Saga dem prep 09.png

Перепроецирование и ресэмплинг[править]

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

Saga dem prep 10.png

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

Saga dem prep 11.png

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

Saga dem prep 12.png

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

Saga dem prep 13.png

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

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

Saga dem prep 14.png
Saga dem prep 15.png

После сообщения 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.

Saga dem prep 16.png

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

Фильтрация[править]

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

Saga dem prep 17.png

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

Saga dem prep 18.png

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

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

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

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

Saga dem prep 19.png

После завершения работы модуля и появления сообщения 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.

Saga dem prep 20.png
Saga dem prep 21.png

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

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

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

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

Гидрологическая коррекция[править]

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

Saga dem prep 23.png

После завершения работы модуля и сообщения 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).

Saga dem prep 24.png

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

Saga dem prep 25.png

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

Saga dem prep 26.png

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

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

Saga dem prep 27.png
Saga dem prep 28.png

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

Saga dem prep 29.png

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

Saga dem prep 30.png

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

Ссылки по теме[править]

  1. Соотношение между точностью данных SRTM и различными картографическими масштабами детально рассмотрено в публикациях:
  2. Jarvis A., Reuter H., Nelson A., Guevara E. Hole-filled seamless SRTM data V.4. International Centre for Tropical Agriculture (CIAT). – 2008. – http://srtm.csi.cgiar.org
  3. Подробнее о работе некоторых фильтров:
    1. Grid – Filter (Perego 2009)
    2. Grid – Filter => DTM-Filter (slope-based), пример
    3. Grid-Filter => Multi Direction Lee Filter
    4. Grid-Filter => Mesh Denoise, пример
  4. Ознакомиться с использованным алгоритмом гидрологической коррекции можно в оригинальной публикации Planchon O., Darboux F. A fast, simple and versatile algorithm to fill the depressions of digital elevation models // Catena, 2002, № 46(2-3), р. 159-176

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

Обсудить в форуме Комментариев — 37Редактировать в вики

Последнее обновление: 2014-05-14 23:47

Дата создания: 13.06.2013
Автор(ы): Дарья Свидзинская