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

Эта страница является черновиком статьи.

Последовательность шагов по подготовке данных SRTM к анализу

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

Импорт данных SRTM в среду ГИС SAGA

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. Это позволит сопоставить результаты как визуально

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

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

Ссылки

Больше про работу с ГИС SAGA