Визуализация рельефа по данным 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 с отмывкой и изолиниями:

Метод Бенфорда в оценке достоверности данных

Друзья мои! Вы несомненно знаете больше меня о последних мировых новостях и потому разобщены и тревожны. Но сегодня, у вас будет повод отвлечься. В этот день мы все объединены единым горем утраты. Утрачена флешка, на которой я хранил для вас статью о диссипативной динамике живого напочвенного покрова. Вместе с ней пропало содержимое подарочной бутылки коньяка, мой рукописный реферат на тему «Сатанизм-как социальное явление» и весь тираж осеннего номера «Лабораторного Журнала», отпечатанный в объеме двух с половиной экземпляров. Воистину, в этот день можно посыпать голову пеплом, ибо об этот реферат я в свое время исписал четыре ручки и мне он чертовски дорог, как память о студенческих годах.

Дабы загладить боль утраты, я предлагаю вам статью из пропавшего «Лабораторного Журнала» (а где вы ее теперь прочитаете?), описывающую сущность, принципы применимости и алгоритм метода Бенфорда на примере анализа данных о площадях ООПТ России и площадях, охваченных лесными пожарами в 2009-2013 годах. Сам же я отправляюсь в келью, где буду страдать вплоть до открытия магазина.

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

Свое название бенфорд-последовательность получила в честь Фрэнка Бенфорда Альберта-младшего — американского инженера-электрика, физика и оптика, жившего в штатах в первой половине XX века. Однако, сам «Закон Бенфорда», он же «закон первой цифры» впервые описан за три года до его рождения американским астрономом, математиком и экономистом Саймоном Ньюкомбом. Работая в 1881 году с логарифмическими таблицами в книгах, он обнаружил, что сильнее всего истрепаны страницы на которых содержаться логарифмы чисел, начинающиеся с единицы. На первый взгляд, вероятность оказаться на первом месте в числе одинакова для всех цифр и составляет 1/9. Однако, чем выше по значению было число, состоящее из первой цифры логарифма, тем в большей сохранности находились страницы. Все это наводило на подозрение о неравномерной встречаемости первых цифр в числах.

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

P(d) = logb(d+1)-logb(d) = logb(1+1/d)

где:
b — основание системы счисления, в нашем случае b = 10;
d — первая цифра в числе;

На основе этой формулы была построена бенфорд-последовательность — последовательность вероятности появления различных цифр на первом месте числа. Рассчитанная по формуле, эта последовательность выглядит следующим образом: 30.1, 17.6, 12.5, 9.7, 7.9, 6.7, 5.8, 5.1, 4.6. Вероятность того, что на первом месте в числе окажется единица составляет 30.1%, двойка — 17,6% и так далее до девятки (4.6%).

Долгое время, эта интересная закономерность не находила никакого применения. Однако после 1997 года на нее обратили внимание и стали все активнее использовать для проверки фальсификации данных, например результатов голосования (в том числе и в России). В 1997 году М. Нигрини и Л. Миттермайер в издании «Аудит: Журнал теории и практики» опубликовали шесть разработанных математических тестов, основанных на законе Бенфорда. Тесты были успешно введены в практику аудиторской компанией «Эрнст и Янг» и позволили выявить несоответствие между реальными и заявленными данными клиентов.

Необходимо учитывать, что метод Бенфорда применим не ко всем данным. Он выдает значительные погрешности при работе с выборками для которых заданы максимальные или минимальные значения, с выборками, охватывающими только один или два порядка величин и с малыми по объему выборками.

При решении вопроса применимости метода Бенфорда обычно рекомендуют исходить из «естественности» данных (если данные получены в ходе естественного течения событий, то к ним применим метод Бенфорда). Этот критерий верен, но довольно сложен для использования. В ходе работ с бенфорд-последовательностями я пришел к выводу, что метод бенфорда работает только с данными, топологическое множество которых самоподобно, а элементы могут принимать произвольные значения.

Для проверки применимости метода необходимо аппроксимировать их показательной функцией (чаще всего используется экспонента) и убедиться, что коэффициент аппроксимации составляет 0,9 и выше. Если при этом отсутствуют правила, детерминантно определяющие значение того или иного числа, то метод бенфорда к вашим данным применим.

Алгоритм применения бенфорд-метода в программах LibreOfficeCalc и MS Excel 

1. Исходные данные

Со страницы сайта oopt.aari.ru, разработанного ФГБУ «ААНИИ» и Лабораторией геоинформационных технологий взят перечень особо охраняемых природных территорий России. Список насчитывает 8013 ООПТ, из которых 4410 войдут в нашу обработку. Это действующие или реорганизованные ООПТ, для которых есть данные по площади.

Данные по площади лесных пожаров взяты с сайта федерального агентства лесного хозяйства. Выборка охватывает данные по всем регионам России с первого квартала 2009 года по второй квартал 2013 года. Всего за этот период было охвачено лесным пожаром 949 территорий различной площади.

2. Проверка на распределение

Нам необходимо убедиться, что данные подчиняются экспоненциальному распределению. Сортируем данные по площади и аппроксимируем их экспонентой.

Lj2-24

На рисунках изображены площади ООПТ (верхний рисунок) и площади пожаров (нижний рисунок), отсортированные по значению. Ось ординат показывает площадь в гектарах.   Чем больше площадь особо охраняемой природной территории, тем меньше таких ООПТ в стране. Равно как и значительные площади подвергаются пожарам гораздо реже небольших участков.  Коэффициент аппроксимации обоих наборов данных экспонентой (синяя линия) составил 0,98.

3. Избавление от нулей

Отличительной особенностью фрактальных множеств, к которым относятся и наши данные является их масштабная инвариантность. Распределение не зависит от единиц в которых выражены величины. Будь наши данные выражены в километрах, миллиметрах или ангстремах, мы всегда будем наблюдать одинаковые закономерности.  Масштабная инвариантность позволяет нам избавиться от значений менее единицы простым умножением на 100 (в каждом конкретном случае может быть различный порядок, в зависимости от наименьшего числа в выборке. В нашем случае таким числом было 0,01). Сделать это необходимо, поскольку формула Бенфорда использует логарифмы, а потому не работает с нулевыми числами.

4. Отделение первой цифры и расчет

Методом LEFT() в LibreOfficeCalc или ЛЕВСИМВ() в Excel отделяем первую цифру из каждого числа. Получившийся столбец с первыми цифрами чисел сортируем и подсчитываем количество единиц, двоек, троек и т.д. до девяток. Вероятность встречи каждой цифры рассчитываем как отношение количества чисел, начинающихся с данной цифры к общему количеству чисел. Например, если в выборке по пожарам было 273 числа, начинающихся на единицу, а общий объем выборки 949, то вероятность того, что первой цифрой в числе будет единица составит 100%*273/949=28,8%.   В итоге у вас получится аналог вот таких таблиц (верхняя таблица — данные по площади ООПТ, нижняя таблица — данные по площади пожаров):

Lj2-25

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

Lj2-252

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

Приведенные графики свидетельствуют, что данные по площадям ООПТ России и данные по площади пожаров за 2009-2013 г. достоверны. Наибольшие ошибки приходятся на крайние значения, что связано со сложностью определения массовых (ошибки по единице) и крупных (ошибки по девятке) объектов в натуре, а также с меньшим объемом статистических данных (ошибки по девятке).

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

P.S. Да, я знаю, что качество приведенных картинок отвратительно. Но поверьте, вы встретились с ними в странный момент их жизни.

Мапим по гуглопанорамам — наземное фотограмметрическое картирование в QGIS с помощью плагина stereoSurveys

Read in English

Замечание 1. В данной статье не расматриваются юридические вопросы законности использования описанного метода при работе с данными компании Гугл. Мое дело метод показать, а с юристами сами разбирайтесь.
Замечание 2. Собственно, ничто не мешает использовать любые другие данные, вплоть до своих фотографий.
Замечание 3. Код описываемого ниже модуля от первой до последней строки написан Enrico Ferreguti. Мое значение в этом проекте чисто терапевтическое.

Пару месяцев назад я опубликовал пост о технологии применении снимков Google StreetView для фотограмметрического картирования территории. Этот незатейливый текст вдохновил Enrico Ferreguti (по его словам) на разработку модуля StereoSurveys для QGIS.

Основное назначение модуля — перенос контуров объектов, попавших в объектив камеры гугломобиля в точечный слой QGis. Наличие объекта с известным пеленгом на двух геопривязанных снимках позволяет однозначно установить его местоположение по свойству суммы углов треугольника. Использовать это свойство на практике возможно было и ранее, однако, только после появления StereoSurveys появилась возможность значительно сократить трудоемкость работ и увеличить точноность нанесения данных.

Для иллюстрации работы модуля, нанесем в точечном слое местоположение фонарных столбов на улице Текстильной (ХБК, город Шахты). Со спутникового снимка их не видно, поэтому единственный способ нанести их не от балды — воспользоваться описываемым методом.

0

Улица Текстильная, вдоль которой необходимо отметить фонарные столбы

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

2. Скачанный архив необходимо распаковать в папку хранения модулей QGIS. Для windows XP это обычно C:\Documents and Settings\[username]\.qgis2\python\plugins для windows 7: C:\Users\[username]\.qgis2\python\plugins, для linux: /home/[username]/qgis2/python/plugins. Если в папке home отсутствует папка qgis2 (или .qgis2), то, возможно, у вас не отображаются скрытые файлы. В debian-подобных системах (debian, ubuntu, mint, OSGeoLive и др.) это исправляется так: правая кнопка мыши — отображать скрытые файлы.

3. Если у вас винда — обратите внимание на название папки [username]. Плагин не переносит кириллицу! Если в названии файла, проекта, плагина или пути к ним встретится русская буква, QGIS выдаст следующую ошибку:

Traceback (most recent call last):   File "", line 1, in   File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\pyplugin_installer\installer.py", line 274, in upgradeAllUpgradeable     self.installPlugin(key, quiet=True)   File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\pyplugin_installer\installer.py", line 322, in installPlugin     reloadPlugin(key) # unloadPlugin + loadPlugin + startPlugin   File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 319, in reloadPlugin     loadPlugin(packageName)   File "C:/PROGRA~1/QGISWI~1/apps/qgis/./python\qgis\utils.py", line 200, in loadPlugin     msg = msgTemplate % (packageName, "', '".join(sys.path)) UnicodeDecodeError: 'ascii' codec can't decode byte 0xd0 in position 9: ordinal not in range(128)

Для пользователей Windows 7, причина кроется как правило в том, что папка «Users» в этой системе названа как «Пользователи», кроме того, обычно кириллическое написание имеет папка  [username]. Теоретически, исправить эту беду можно создав новую учетную запись с правами администратора. Через нее следует войти завершив текущую сессию, после чего стандартным способом сменить название директорий на латиницу. Насколько это возможно и действенно, я утверждать не берусь — не проверял.

В Линуксе, как правило, папка [username] всегда может быть названа только латиницей. Поэтому, если у вас Линукс — переходите к следующему пункту.

4. Теперь можно запустить QGIS, либо перезапустить его, если он работал до этого.

В верхнем меню выбираем: модули — управление модулями-с ошибками. Ставим галку в чекбоксе stereo surveys и нажимаем установить

Снимок экрана от 2015-06-20 21_37_45

 

После установки появляется текст «Модуль неисправен invalid syntax». Презираем его и закрывая окно установки модулей.

Снимок экрана от 2015-06-20 21_37_59

 

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

Снимок экрана от 2015-06-20 21_40_00

На панорамах Google отображается площадь Прато-делла-Валле — что находится в Падуе, недалеко от Венеции. Это красивейшие места, но в качестве отправной точки для модуля не самые удачные — нам ведь требуется картировать шахтинскую улицу.

Для точечного слоя, создадим шейп-файл с названием test и системой координат EPSG: 3785 (можно и EPSG: 3857 — все работает, другие датумы не пробовал).

Снимок экрана от 2015-06-20 21_42_06

 

Что-бы вернуть его в рамки окна — перетащим таблицу слоев (менеджер слоев, TOC — все его по-разному называют) в нижнюю часть экрана.

Для того, что-бы картирование было более наглядным, подгрузим слой OpenStreetMap — Mapnik, воспользовавшись плагином OpenLayersPlugin.  Это не обязательно, но очень удобно. Приблизим нужный участок.

Рядом с левым верхним углом левой панорамы находится кнопка с зеленым индейцем. Нажимаем на эту кнопку, после чего выбираем на карте точку, из которой мы хотим видеть панораму. Обычно после этого в левом окне появляется панорама Google StreetView.

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

Первое сообщение об ошибке:

Снимок экрана от 2015-06-20 21_43_57

 

Второе сообщение об ошибке:

Снимок экрана от 2015-06-20 21_44_00

 

После нескольких неудачных попыток модуль срабатывает как надо и пользователь видит следующее:

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

— Непонятные числа, поля ввода и текст not calibrated слева и справа от окон с панорамами (первые два числа слева и справа — это координаты точки, в которой сделана панорама.

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

H — высота точки, взятой на прицел в обоих панорамных окнах в метрах, относительно земли;

+/- — погрешность местоположения точки, взятой на прицел в обоих окнах, в метрах;

Lon, Lat — долгота и широта;

X,Y — неизвестные прямоугольные координаты;

— Под числами, расположена кнопка с текстом, «Digitize on map», по клику на которой, точка, взятая на прицел в обоих панорамах должна переносится в точечный слой.

— В окне отображения данных QGIS, помимо OSM-карты, пользователь видит три точки: зеленую — точка из которой снята панорама в левом окне StereoSurveys, красную — точка, из которой снята панорама в правом окне и желтую — точка, взятая на прицел в обоих окнах просмотра панорам.

Можно начинать работу. Находим на обоих панорамах одну и ту же точку. Кнопка Digitize on map не срабатывает, поэтому, включаем редактирование слоя и наносим ее вручную (кликаем по желтой точке, предварительно выбрав инструмент «добавить объект»). В проекте, который идет как пример к плагину, кнопка Digitize on map работает, но тоже далеко не всегда. Этот вопрос еще необходимо прояснить.

Находим один и тот же объект на двух снимках (в нашем случае это первый фонарный столб со знаком пешеходного перехода).

Снимок экрана от 2015-06-20 21_44_18

 

Прицеливаемся поточнее и стреляем — ставим точечный объект поверх желтой точки.

Снимок экрана от 2015-06-20 21_45_02

 

Точно так-же поступаем для остальных объектов. Картирование фонарного столба на противоположной стороне улицы.
Снимок экрана от 2015-06-20 21_47_33

 

Навигация осуществляется средствами Google StreetView. «Проезжаем» на несколько метров вперед в обоих окнах просмотра панорам и берем на прицел новый столб.

Снимок экрана от 2015-06-20 21_49_59

 

Иногда, без видимых причин возникает сообщение об ошибке:

Снимок экрана от 2015-06-20 21_51_10

 

Игнорируем его и продолжаем работать дальше.

Снимок экрана от 2015-06-20 21_53_17

 

Ради интереса, попытаемся изменить числа в полях StereoSurveys. Числа ввести можно только корректные. Насколько я понял из общения с Enrico, данные поля позволяют корректировать данные о высоте точки. Первое поле отвечает за высоту камеры над поверхностью земли (2.5 метра), второе за высоту поверхности земли (по умолчанию 0). Эти параметры особенно важны в горной местности с большими перепадами высот и при картировании объектов, находящихся на большом удалении от точек съемки панорам.

Не отвлекаемся и стреляем столбы дальше.

Снимок экрана от 2015-06-20 21_58_20

 

Процесс захватывает.

Снимок экрана от 2015-06-20 22_00_02 Снимок экрана от 2015-06-20 22_04_23 Снимок экрана от 2015-06-20 22_08_46 Снимок экрана от 2015-06-20 22_09_57

 

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

Еще точечку.

Снимок экрана от 2015-06-20 22:13:22 Снимок экрана от 2015-06-20 22:23:21 Снимок экрана от 2015-06-20 22:25:19

 

В результате, мы получаем точечный слой, который планировали. После нанесения всех точек, закрываем модуль и сохраняем полученный слой. Иногда, после закрытия модуля красная и зеленая точки (точки съемки панорам) не исчезают, в этом случае, необходимо отключить модуль через меню-модули-управление модулями-снять галку с чекбокса напротив StereoSurveys и нажать «Закрыть».

Любопытный момент: обычно, снимки привязывают к gps-трекам. Я за всю Одессу говорить не буду, но в OpenStreetMap так и происходит. Думаю в Google тоже. При этом трек рассматривается в качестве центра дороги, что частично верно только для проселочных колейных дорог. Игнорирование понятия полосы дороги, приводит к системной ошибке в координатах до 5 метров (обычно +/- 3 м). Сама по себе, это величина небольшая, часто незначительная на фоне погрешности прибора. Однако, в том случае, если координаты объекта используются для расчетов детерминированных показателей, ошибки могут быть колоссальны. Это категорически запрещает использование популярных веб-проектов для расчета площади секторов, лежащих в створе с известными координатами, анализа видимости объектов, привязки объектов по двум точкам…, боюсь, этот список может оказаться длинным.

Обратите внимание на местоположение фонарного столба на улице Ворошилова.

Карта OpenStreetMap — столб стоит на дороге.

Снимок экрана от 2015-06-20 22:31:49

 

Карта Google. Аналогичное местоположение столба. Так-же, обратите внимание на местоположение дороги, относительно столбов освещения.

Снимок экрана от 2015-06-20 22:31:14

 

Загадка. В какую сторону ехал гугломобиль, если дело происходит в стране с правосторонним движением?

Снимок экрана от 2015-06-20 22:30:49

 

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

Вот карта фонарных столбов, полученная с помощью модуля:

Снимок экрана от 2015-06-20 22:27:30

 

А вот, та же карта, полученная два месяца назад вручную. Ошибки нанесения (дублирование столбов) — налицо.

9

 

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

- Не работает с кириллицей;
- Не хватает угла обзора панорам как в модуле go2streetview;
- Возникает ошибка при попытке переместить точку на карте;
- Не работает кнопка Digitize on map;
- Плагин не реагирует, когда на обоих экранах одинаковый снимок (должно высвечиваться предупреждение);
- Не высвечивается предупреждение, если погрешность больше заданной величины;
- При запуске модуля должны отображаться снимки из текущих координат экстента;
- При наличии таблицы слоев, модуль выходит за размеры экрана;
- Не исчезают точки после закрытия плагина;
- Нет явного объяснения значения чисел и полей ввода в плагине;
- Много времени уходит на масштабирование снимков, не хватает кнопки "сбросить увеличение";
- Не хватает кнопки "назад" для каждого из окон панорам;
- Не хватает кнопки "вперед", перемещающей обе точки;

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

Наземное фотограмметрическое картирование объектов с помощью Google StreetView, QGis и модуля go2streetview

Замечание 1. В данной статье не расматриваются юридические вопросы законности использования описанного метода при работе с данными компании Гугл. Мое дело метод показать, а с юристами сами разбирайтесь.
Замечание 2. Собственно, ничто не мешает использовать любые другие данные, вплоть до своих фотографий. Для этого необходимо только переписать модуль. К сожалению, моих познаний в языке для этого недостаточно. Была попытка написать отдельную программу для реализации этого метода (Borland С++ Builder 6) с произвольными фотографиями, но, ввиду отсутствия времени и денег работа была остановлена.

 

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

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

0

 

С помощью меню «Модули» -> «Управление модулями» установим модуль «go2streetview«. После чего откроем в QGis указанный район и запустим модуль. Сделать это можно перетаскиванием желтого человечка на необходимую точку улицы, точно так же как и в Google Maps. В новом окне откроется фотография улицы из указанной точки, а сама точка на карте будет обозначенна синим цветом с двумя лучами, иллюстрирующими рамки окна с фотографией.

1

 

Изменение угла обзора в окне go2streetview, автоматически изменяет положение лучей. Это замечательное свойство модуля, наряду со свойством треугольника, позволяет нам определить местоположение любого объекта, попадающего не менее чем на два снимка StreetView.

В окне панорамы приблизим первый столб и соориентируем вид таким образом, что-бы его пересекала рамка окна. После чего создадим вспомогательный линейный слой, в котором продлим продлим луч go2streetview. В моем случае, таким слоем является слой «водотоки» (по ленному недоразумению). К сожалению, модуль go2streetview не позволяет указать в настройках размеры лучей, поэтому приходится поступать таким методом.

2
Такую же процедуру повторим для столба на противоположной части дороги. Для удобства, можно использовать правую границу окна и правый луч соответственно.

Теперь откроем панораму с этими столбами снятую из другой точки. Сделать это можно кликнув по карте (не забудьте перед этим переключиться из режима редактирования слоя в режим модуля). Еще проще сделать это кликнув по стрелке в окне go2streetview.

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

3
Находим противоположный столб.

4

 

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

5

 

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

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

6

 

Таким методом мы движемся вдоль всей улицы, редактируя точечный слой после каждого определения местоположения столба.

7

 

Получая в итоге следующее:

8

 

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

9

Наземный фотограмметрический WZPr-метод оценки полноты древостоя

«Лесник без спичек, что хуй без яичек»
Народная лесниковская поговорка

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

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

Матчасть в лице учебника таксации говорит, что полнота «представляет собой сумму площадей поперечных сечений всех деревьев на площади на высоте 130 см в пересчете на гектар леса». Естественно, тут и далее в статье пойдет речь про абсолютную полноту, которая выражается обычно в квадратных метрах на гектар. Проще говоря абсолютная полнота это площадь всех пеньков на гектаре, срубленных выше пупка, но ниже головы — одна из важнейших величин при определения запаса древесины, а следовательно всех вытекающих параметров. Для измерения абсолютной полноты существует давно разработанный и опробованный метод угловых проб, который чаще называют методом реласкопических площадок или методом Биттерлиха, по имени изобретателя. В старой литературе, иногда встречается название WZPr-метод, от немецкого die Winkelzahlprobe. Те, кто уже знаком с методом, могут пропустить нудятину с его описанием и сразу перейти к сути вопроса.

Метод очаровательно прост. Вам потребуется изготовить вот такой «прибор»:

Полнотомер

Фактически это любая ровная палка с закрепленной на ее конце пластинкой с прорезью. Можно использовать любые материалы, главное, что-бы отношение прорези к длине палки составляло ровно 1:50. Периодически встречаю у коллег разные китайские поделки, вроде вот такой:

Нормальный полнотомер

Здесь роль палки выполняет натянутая цепочка, длина которой в 50 раз больше ширины прорези. Но большая часть использует то, что есть в карманах: спички, зажигалки, куски картона или «ключи» от пивных банок (у кого что лежит в карманах). Главное не ошибится с размерами, после определенной тренировки, такие «инструменты» зачастую дают лучший результат чем адски дорогая лазерно-оптическая техника.
Вы берете прибор Биттерлиха (полнотомер) и прикладываете конец палки без прорези, либо свободный конец цепочки к глазу и визируете через него на деревья, поворачиваясь по кругу. Первым визируете самое ближнее к вам, что-бы не спутать его с остальными и не пойти на «второй круг». При визировании вы смотрите вдоль палки/цепочки сквозь прорезь и наблюдаете три ситуации:

1. При наведении на дерево, дерево «закрывает» прорезь
2. При наведении на дерево, оно точно вписывается в прорезь (невозможно точно определить «закрывает» дерево прорезь или нет)
3. При наведении на дерево, дерево очевидно не закрывает прорезь

Принцип полнотомера

Деревья из первой категории вы считаете за единицу, деревья из второй категории считаете за пол-единицы. Третью категорию не учитываете. Предположим у вас 10 деревьев вошло в первую категорию, 4 дерева попали во вторую и все остальные попали в третью категорию. Абсолютная полнота в этом случае составляет 10*1+4*0,5=12 кв. метров на 1 гектар. Естественно, одной площадкой не обойдешься — весь мануал доступен в книжке «Методы отвода и таксации лесосек» — рекомендую к прочтению.

Весь фокус метода можно понять из картинки ниже:

Хуевая картинка

если длина полнотомера составляет 100 см, то круг, описанный им будет иметь площадь π кв.м. Если в этом кругу будет расти дерево, которое точно вписывается в разрез пластины полнотомера, значит диаметр этого дерева 2 см, а площадь сечения 0,0001π кв.м., что соответствует абсолютной полноте 1 кв. м на 1 га.

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

Призма Анучина

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

Принцип призмы Анучина

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

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

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

Схема метода

Не правда ли, это похоже на полнотомер Биттерлиха, только как бы «повернутый» прицельным визиром к наблюдателю? Естественно, ни о каком соотношении 1:50 тут не может быть и речи, но тем не менее, это не мешает установив зависимость между линейными размерами пиксела, и углом, который он охватывает в натуре, оценивать значения абсолютных полнот.
Первое что необходимо сделать — это узнать ширину 1 пиксела. Ее можно посмотреть либо в свойствах снимка, либо с помощью фоторедактора, например в GIMP 2.8: меню изображение -> свойства-разрешение:

Гимп скриншот

Разрешение снимка обычно представлено либо одним, либо двумя числами. В первом случае чило означает количество пикселей в дюйме диагонали, во втором случае, количестово пикселей в дюйме ширины и высоты. Соответственно, для вычисления ширины в первом случае нам необходимо воспользоваться теоремой Пифагора (количество пикселей в дюйме ширины равно квадратному корню половины квадрата количества пикселей в дюйме диагонали). В нашем случае в 1 дюйме ширины содержится 180 пикселей изображения. Следовательно, ширина 1 пиксела 0,014 см. Это и будет минимальная ширина прорези в нашем виртуальном полнотомере Биттерлиха.
Теперь необходимо определиться с длиной виртуальной базисной рейки. Для этого необходимо установить какой угол по ширине охватывает объектив вашего фотоаппарата. Учтите, что он изменяется при зуммировании объектива, поэтому «приближения» при снимках использовать не нужно.
Ширину охвата фотоаппарата можно рассчитать исходя из технической документации на ваш объектив, в которой указан угол изображения объектива. Для широкоугольных объективов он превышает 75 градусов, для длиннофокусных составляет менее 30 градусов. Остальные объективы называют нормальными, как правило на недорогих моделях стоят именно они. Угол изображения объектива — это диагональный угол, поэтому для расчета ширины охвата требуется вновь воспользоваться теоремой Пифагора.
Если же документация утеряна много лет назад, как в моем случае, достаточно просто установить фотоаппарат на штатив, сфотографировав удаленную на известное расстояние линейку, после чего измерить получившийся угол. Мой фотоаппарат охватывает ширину в 45 градусов.
Ширина фотографии составляет 3072 пиксела, что видно из свойств растра. Значит на каждый пиксел охватывает ширину в 45/3072=0,0146484375 градуса. Выше мы рассчитали, что каждый пиксел имеет ширину 0,014 см. Все что нам остается, для того, что-бы рассчитать длину виртуальной базисной рейки — это решить задачку для седьмого класса: основание равнобедренного треугольника составляет 0,014 см, вершина имеет угол 0,0146484375 градуса. Всего-то необходимо найти высоту этого треугольника.
Решение 1. Высота делит равнобедренный треугольник на два равных прямоугольных треугольника с длиной противолежащего катета 0,014/2=0,007 см и острым углом в вешине 0,0146484375/2=0,00732421875 градуса. Искомая высота является прилежащим катетом в наших треугольниках и равна отношению катета к тангенсу острого угла: 0,007 / tg(0,00732421875) = 54,7594860417≈55 см.
Решение 2. Будем считать, что этот треугольник не равнобедренный, а прямоугольный. Да, это аморальное математическое допущение, но при таком соотношении сторон, возникающая ошибка незначительно мала. Тогда искомая высота, она же прилежащий катет составляет: 0,014/tg(0,0146484375)=54,7594851469≈55 см.
Размер нашей виртуальной базисной рейки составляет 5,5 м. Визирный прицел составляет 0,00014 м. Следовательно, если на фотографии мы видим, что единственное дерево имеет ширину 1 пиксел, мы можем утверждать, что на площади (π×5,5^2)/(360/45)=11,8791472214≈12 кв.м площадь сечения древостоя составляет π*0,00007^2=1,5393804 × 10^(-8) м, что в пересчете на гектар составит 0,000013 кв.м. Это, конечно, чрезвычайно мало, но ситуация, при которой у вас в кадре будет только один ствол дерева шириной в 1 пиксел может быть разве что в таком случае:

Саванна

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

Таким образом, абсолютная полнота древостоя на 1 га вычисляется по формуле:
P = Σ(10000*(π*(0,5*B)^2)/(360/V*π*((0,5*B)/tg(0,5*0,W*B))^2)) где,
P — полнота, в кв.м
V — угол охвата снимка,
B — ширина дерева на фото, в см (ширина прорези виртуального полнотомера),
W — угол охвата одного пиксела, градусов
Применительно к нашему случаю:
P = Σ(10000*(π*(0,5*B)^2)/(360/45*π*((0,5*B)/tg(0,5*0,014*B))^2))
В материальном виде, использование этого метода можно представить так, будто мы располагаем не одним полнотомером, а бесконечным множеством полнотомеров с разными размерами визиров. Классический метод Биттерлиха совершенно игнорирует все деревья, не совпадающие с размером визира. В случае же фотограмметрического метода, учитываются все деревья. При этом, нет необходимости измерять абсолютно все деревья на снимке, но чем больше их будет измерено, тем точнее окажутся результаты.