Визуализация рельефа по данным SRTM и ASTER GDEM в QGis+SAGA

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

В качестве показательной территории взяты Соликамский и Красновишерский районы Пермского края. В качестве подложки карта OpenStreetMap Mapnik Standart:
Соликамск и Красновишерск

Инструментом для работы послужит QGis 2.18.4 с подключенными алгоритмами SAGA. Все операции, связанные с фильтрацией и созданием изолиний можно выполнять как внутри QGis, так и запуская гис SAGA в качестве автономного приложения.

Рельеф на средне- и крупномасштабных картах в настоящее время в большинстве случаев отображается с помощью данных SRTM или ASTER GDEM, что связано с их глобальным охватом, открытостью и простотой получения. Разрешение этих данных (SRTM 90 м/пикс, ASTER GDEM 30 м/пикс) позволяет, при должной обработке, показывать особенности рельефа примерно до 15 зума. Несмотря на то, что данные ASTER точнее, их использование затруднено необходимостью дополнительной фильтрации для отсеивания значений, не отражающих реальный рельеф (например, высоты леса и жилой застройки). Оптимальных алгоритмов для такой процедуры, которые дают стабильный результат для значительной территории, не разработано, в результате чего, образец визуализации менее точных данных SRTM оказывается обычно более качественным как с геодезической, так и с художественной точек зрения. Однако, севернее 60° с.ш. и южнее 54° ю.ш. данные SRTM отсутствуют, что вынуждает в конечном итоге использовать оба набора данных при визуализации рельефа на территориях, выходящих за границы покрытия SRTM.

Наш случай именно такой (снизу данные SRTM, сверху ASTER GDEM):
ASTER и SRTM

Данные SRTM доступны из различных источников, из которых наиболее удобны сайты cgiar, gis-lab и viewfinderpanoramas. Я предпочитаю использовать последний, поскольку многие сцены там объединены и загружены сразу в geotiff-растры (обычно SRTM представлен в hdt-формате).

ASTER получить немного сложнее: сайты геологической службы США и NASA позволяют скачивать различные данные ДЗЗ, что требует от пользователя определенной подготовки. Кроме того, эти сайты иногда бывают недоступны, либо работают чрезвычайно медленно. В этих случаях можно скачать всю базу через торрент. Дополнительные источники получения данных SRTM и ASTER доступны на странице получения данных.

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

После того как исходные получены, можно запускать QGis:
ASTER и SRTM в QGis

Для начала объединим сцены ASTER в единый растр с помощью меню «растр-прочее-объединение»:
1объединение растров

В диалоговом окне укажем директорию со сценами и название итогового файла.
2объединение растров

Обратите внимание, что в некоторых версиях GDAL отказывается работать, если пути к файлам содержат кириллические символы. В моем случае все прошло успешно:
3объединение растров

Теперь сохраним полученный файл уменьшив его разрешение до разрешения SRTM. Если этого не сделать, в месте соприкосновения сцен из разных источников мы получим вот такую картину:
рельеф SRTM и ASTER

Выделяем в ТОСе наш ASTER, и через правую кнопку мыши вызываем диалоговое окно сохранения растра:
Сохранение растра в QGIS

Здесь обратите внимание на то, что-бы растр сохранялся как данные. Разрешение уменьшаем в три раза, т.е. вместо 18001 столбца вписываем 6000, а вместо 7201 строк вписываем 2400:
Сохранение растра в QGIS2

После сохранения растра он выглядит загрубленным, но все-равно более информативным, чем SRTM:
Сохранение растра в QGIS3

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

По той-же схеме объединим полученный растр с растром SRTM:
DEM-композит

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

image

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

Откроем панель инструментов в меню «анализ данных» и в списке геоалгоритмов SAGA-Grid-Filter выберем алгоритм «Простой фильтр»:

После нескольких минут обработки, алгоритм выдаст сглаженный растр:
отфильтрованный растр

Его мы и будем использовать для отмывки. С помощью меню «растр-морфометрический анализ-теневой рельеф»:

вызовем диалог создания карты теней:

На этом этапе следует знать одно неявное правило. В том случае, когда вы планируете использовать теневую отмывку в качестве подложки саму по себе, можно использовать значения по умолчанию. В том случае, когда ваша отмывка ложится на некую базовую карту (в нашем случае, такой картой служит OpenStreetMap), следует повернуть источник освещения на сто восемьдесят градусов. Дело в том, что стандартная отмывка представляет собой темный растр, который при наложении не только перестает читаться, но и зашумливает перекрываемые слои. Для того, что-бы это избежать, отмывку следует инвертировать, но в этом случае, горы выглядят как впадины, а каньоны напоминают холмы. Учитывая это, мы заранее изменяем источник освещения, что позволит нам при инвертировании сохранить визуальную форму отмывки. По умолчанию, источник освещения расположен в районе трехсот градусов, чего, конечно-же в природе почти никогда не бывает. Еще Салищев указывал на эту особенность — привычная отмывка рельефа обязана своему появлению лампам, которые обычно устанавливали слева от кульмана. Мы поменяем значение «300» на «120» и через несколько секунд алгоритм выдаст нам вот такой результат:

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

Через меню «растр-извлечение-обрезка»

вызовем диалог, в котором укажем исходный и результирующий растр и шейп-маску по которой будет произведена обрезка:

В результате получим вот такую картину:

Двойным кликом по полученному растру в ТОСе вызовем меню свойств, где сменим градиент с «от белого к черному» на «от черного к белому». После применения изменений растр инвертируется. В месте сочленения данных ASTER GDEM и SRTM осталась небольшая белая полоса, однако, после того как будет установлена прозрачность, а сама отмывка наложена на подложку, заметить эту полосу будет практически невозможно:

Для того, что-бы не инвертировать растр при каждой новой загрузке, сохраним его как новый слой, но на этот раз в меню сохранения отметим его не как «данные», а как «изображение». На этом создание отмывки закончено. Установим прозрачность отмывки 95% и подложим под нее OpenStreetMap:

Так выглядит чистая карта OSM:

А так выглядит карта OSM с наложенной на нее картой теней в районе соприкосновения данных SRTM и ASTER:

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

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

С практической точки зрения эта проблема решается так. Создадим временный полигональный слой:
Новый временный слой2

Сделаем его редактируемым (иконка желтого карандаш), после чего установим режим добавления нового объекта
Редактирование в QGis

и обведем один из районов:

Сохраним фрагмент растра для этого участка (не фильтрованного растра, а исходного, мы же усложняем себе задачу!), обрежем его по обведенной области и отфильтруем с теми же параметрами, что и при создании отмывки. Затем с помощью геоалгоритма SAGA-Shapes-GRID-Горизонтали по ЦМР создадим изолинии через каждые 50 метров высоты.

Фильтрация не только убирает излишний шум, упрощая и выравнивая изолинии, но и позволяет «сцепить» наши разнородные данные. Вот пример извлечения изолиний из сырого растра:

Отчетливо просматривается линия сочетания данных ASTER и SRTM. При различных способах фильтрации растра ASTER GDEM эту линию можно делать более или менее заметной о чем я упоминал в начале данной статьи.

Изолинии из отфильтрованного растра на этот район выглядят так:

На границе растра изолинии замыкаются и не несут в себе географического смысла. Такие участки в последующем будут удалены. Именно поэтому обрезка dem-растра производилась нами не по границе района а по внешнему слою.

Аналогичные операции повторяем для второго района. Обратите внимание, что полигоны обрезки растра перекрывают друг друга:

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

После того, как изолинии извлечены, остается только обрезать их по контуру границ, сохранив оригинальные значения атрибутов изолиний. Для этого успешно применяется алгоритм «SAGA-Shapes-Lines-Пересечение линий и полигонов»:

Небольшая настройка стиля и изолинии готовы:

Обычная карта OpenStreetMap:

Тот же фрагмент карты, но с наложенной картой теней и горизонталями:

Отдельно необходимо вынести проблему неоднородности данных ASTER GDEM по качеству. Даже на нашем примере видно, насколько сильно артефакты отсутствующей и ошибочной информации сказываются на качестве визуализации рельефа в целом:

Данная проблема не имеет однозначного механистического решения. Наиболее оптимальные способы ее устранения зависят от требований к визуализации, выбранного региона и доступной вычислительной мощности. В качестве одного из способов решения я предлагаю использовать последовательное применение фильтра DTM (предельный уклон местности 10 градусов, радиус поиска 2 пикселя), заполнение пропусков в образованном в результате DTM-фильтрации растра (порог напряжения 10) и последующая фильтрация простым фильтром (круговой режим поиска, гладкий фильтр, радиус 5px). Этот метод не позволяет полностью избавиться от артефактов, но существенно снижает их число и сглаживает, что определенно положительно сказывается на визуализации рельефа:

Карта OpenStreetMap без отмывки и изолиний:

Карта OpenStreetMap с отмывкой и изолиниями:

Дорога на Морье

Оргазм, вибратор, онанизм, жесткач, котята, индуизм

Ладно-ладно, индуизм добавил только для рифмы, он тут и нахуй не припекся, в отличие от остальной части заголовка. Как вы поняли, сегодня я буду говорить о важном и потому хотел бы донести свою мысль до максимального количества людей.

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

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

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

Да хули пиздеть-то? Давайте лучше прогуляемся с вами на мыс Морье, что в Ленинградской области. Там очевидно какой-то поселок:
мапсми

Сам поселок может и не очень интересен, а вот посмотреть на церковь может быть весьма любопытно:

Да и чего бы туда не сходить, если даже роутинг туда строится без проблем?

Дорога великолепна. Слева слабодренированные кислично-сфагновые березняки, осинники и сероольшатники:
мелиоративная канава

Справа за соснами шумит едва вскрывшееся ото льда Ладожское озеро:
Ладога

Под ногами вполне сносная асфальтовая дорога. Идете вы такой, природе радуетесь. Километр, второй, четвертый, пятый и тут, хуяк!
Стой, стреляют

Что за хуйня? Если кому не видно — надпись на щите гласит: «Стой, стреляют. Проезд и проход запрещен». Бля, да там же поселок по мапсми? А давайте-ка глянем OSM-стандарт, геоданные-то одни и те же:

Мыс Морье

Воттыжблядь-то! Так и есть, я шестой километр иду по вомнскому полигону. А впереди у меня кпп воинской части не то связистов, не то РВСН. Так какого хуя в масми об этом ни слова? Хотя постойте, видите эту тараканью залупу:

Ее просто нужно вовремя увидеть. Она просто появляется не на всех зумах. Просто инженер, который принял решение о таком обозначении запретных и опасных зон — мудак.

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

И мне еще повезло, что я был в курсе куда эта дорога ведет, ибо в 2007 году мы с диаволом на пятнашке там уже приключались, за что едва не получили пиздюлей. Тогда столбики «Стой, стреляют» стояли через каждые пятьдесят метров. За десять лет их просто все выпилили к хуям.

отрисовка тайлов

Весенние поводы

Раз уж случилась весна, это отличный повод что-нибудь замутить. Например улучшить картостиль «Y», призванный представить данные OpenStreetMap в стиле Яндекс-карт. Помните, я писал о нем прошлой осенью? Яндекс проделал за год большую работу, добавил новый зум (теперь их девятнадцать), исправил многие ошибки и неточности стиля. Отличный повод для того, что-бы открыть старый xml-файлик и поработать над ним под музыку из японского мультфильма про тетрадь смерти.

Некст-Гис ввел необходимость регистрации для скачивания шейпов OSM и потребовал за свежие выгрузки триста рублей. Цена небольшая, но это отличный повод переписать стиль целиком на работу с чистыми данными. И поддержка легче и данные обновлять проще. А еще добавить уровни дорог, новые иконки, исправить ошибки и, конечно же, расширить линейку зумов до девятнадцати.

Картостиль Y

Картостиль Яндеск-карт

Сервер не резиновый, так что не обессудьте — пока загрузил только Москву и Московскую область. Если напишете в комментариях про расхождения в стилях — буду рад и постараюсь учесть при следующем обновлении.

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

Рыбалка в Кавголово

Это отличный повод где-нибудь поприключаться. Например в Нижнем Новгороде, где я планирую покартировать деревья вдоль Волги и Оки.

Ну а раз уж я буду в Нижнем, предлагаю устроить какую-нибудь лекцию, картопати или другую движуху. Выберем вечер, засядем в приятном месте, познакомимся и прекрасно проведем время. Я готов ответить на все известные мне вопросы по картографии, геоданным, проекту OpenStreetMap, программе QGis, картографическому дизайну, Mapnik-у, Leaflet-у, геоботанике, географии, полевым измерениям, нечеткому тегированию, фракталам, клеточным автоматам и … ну вы поняли. Главное свяжитесь со мной предварительно по телефону +7-904-614-68-29 или почте schwejk-rpnt@rambler.ru. Выберем удобное время и место.

Для того, что-бы что-то замутить не нужен повод. Но грех чего-нибудь не замутить если повод имеется.