Полигоны и лесополосы

Разделение лесов, кустов и полей на старых ландсатах при помощи растра высот

Вот смотрите. Берем седьмой Ландсат за июль 1999 года на территорию Игнатьевского сельского поселения в Адыгее и строим по нему растр вегетационного индекса. Он же — стандартный NDVI, который считается по формуле (nir-red)/(nir+red), где nir — инфракрасный канал, а red — красный. В результате имеем вот такую хренотень:
NDVI-index

А теперь берем снимки Alos, которые улучшенный Астер, который, в свою очередь, отснял эту же территорию несколько месяцев спустя:
Alos

Теперь открываем JOSM и рисуем осевые линии лесополос:
JOSM

Которые затем буферизуем в QGis-e до ширины 100 м (что-бы ухватить разные пикселы Ландсата) и режем на стометровые отрезки поперек. Это в QGis требует особой акробатики, но мы справляемся. Самое главное, что-бы квадраты обязательно выходили за границы лесополосы:
лесополосы

Теперь считаем зональную статистику по двум растрам. С данных Alos снимаем размах, который соответствует разности максимальной и минимальной высоты в каждом полигоне (перепад высот). С растра NDVI снимаем среднее значение индекса в каждом полигоне. Экспортируем все это во внешнюю таблицу и строим примитивный график:
график

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

Остается только повторить процедуру в обратном порядке и получим простейшую карту растительности.

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

Мультифрактальный анализ данных OpenStreetMap в QGis

Бенуа Мандельброт в своей «Геометрии» утверждает, что понятие «фрактал» применимо исключительно к множествам, которые проявляют свойства самоподобия не менее чем на протяжении трех порядков единиц измерения. Манфред Шредер в не менее известной монографии описывает принцип расчета корреляционной размерности, приводя ее к общему виду. Там же он наглядно показывает, как неэвклидовы размерности отражают строение систем в стадии детерминированного хаоса. Аналогичным, но в более популярной форме занимается Хакен, а Федер постулирует чрезвычайно занятную мысль о взаимосвязи самоподобия с теорией перколяции. Однако, несмотря на все это, фракталом называют любой каскад дихотомических ветвлений, а под фрактальным анализом — банальное определение размерности покрытия, забывая даже то, что в оригинальном мануале она идет под названием размерности Хаусдорфа-Безиковича. Еще в одной статье на эту тему нет смысла, но процедура оценки сложности геоданных в QGis настолько проста, что заслуживает краткого упоминания даже среди всеобщей скуки.

Начнем с того, что выберем интересующую область (в моем случае это Москва, район трех вокзалов), которую обводим полигоном. Любой тулзой скачаем OSM данные на этот регион — я использовал плагин OSMDownload, но можно и просто выкачать все через overpass. После того, как данные получены, я рекомендую перевести все в спроецированную систему координат, например в EPSG:2705. Это облегчит дальнейшую работу и позволит вам избежать необходимости перевода градусов в метры для объяснения полученного результата. Я для упрощения работы использовал только точечные данные, но нет препятствий к применению этого метода для линий или даже полигонов.

После этого стандартными средствами (Вектор-выборка-регулярная сетка) строим сетки из квадратов с разной длиной стороны. Чем больше сеток и сильнее разброс площадей у их ячеек тем интереснее результат, но на практике обычно получается не более 20, а если удваивать сторону квадрата для каждой новой сетки, то и того меньше. Можно объединить все сетки в один шейп, это ускорит работу по подсчету заполненных ячеек, но замедлит расчет размерности, так что особого смысла нет. Сам подсчет ведется путем пространственной выборки полигонов сетки по принципу пересечения с точкой OSM:

Я использовал 22 сетки со сторонами квадратов 5, 10, 13, 20, 21, 34, 40, 55, 80, 89, 100, 144, 160, 200, 233, 300, 320, 377, 400, 500, 610 и 640 метров. Не стоит удивляться размерам — эту работу я проводил в рамках изучения встречаемости последовательности Фиббоначчи в геоданных и связи этой последовательности с размерностью покрытия данных. В результате были получены следующие данные:

Сторона
квадрата, м
Полное количество квадратов Количество заполненых
квадратов
Достоверность аппроксимации
(R2)
Размерность
Хаусдорфа-Безиковича
5 241001 4322 1,00 0,00
10 60501 3963 1,00 0,13
13 35574 3755 0,99 0,14
20 15251 3115 0,89 0,22
21 13728 3041 0,91 0,25
34 5251 2151 0,86 0,35
40 3876 1830 0,87 0,41
55 2035 1290 0,87 0,50
80 988 760 0,86 0,62
89 782 650 0,88 0,69
100 651 535 0,89 0,74
144 294 280 0,89 0,83
160 247 237 0,90 0,89
200 176 148 0,90 0,96
233 117 115 0,91 1,02
300 77 69 0,91 1,08
320 70 68 0,92 1,12
377 48 47 0,92 1,16
400 48 40 0,93 1,19
500 35 24 0,93 1,24
610 20 20 0,93 1,27
640 20 19 0,94 1,29


График функции количества заполненных клеток от масштаба клетки как и ожидалось, имеет степенной вид (R2функции аппроксимации = 0,94). Расчеты выполнены в экселе с использованием формул для расчета достоверности авппроксимации:

=ИНДЕКС(ЛИНЕЙН(LN($D$2:$D23);LN($B$2:$B23);;1);3;1)

Для расчета размерности Хаусдорфа-Безиковича:

=-1*ИНДЕКС(ЛИНЕЙН(LN($D$2:$D23);LN($B$2:$B23));1)

На основе полученной таблицы построен график мультифрактального спектра, отражающий изменение сложности данных в зависимости от масштаба (сложность точки = 0, сложность прямой = 1, сложность плоскости = 2):
график мультифрактального спектра

Сложность точечных данных OpenStreetMap возрастает по логарифмическому закону. При увеличении охвата, она стремительно возрастает до тех пор, пока сторона охвата не составит примерно 100 м (вероятно это связано с квартальной сетью), после чего скорость усложнения падает. Если продолжить это график, он должен опять пойти на спад, что следует хотя-бы из здравого смысла, но ресурсов моего компьютера для таких расчетов уже не хватит. Динамика наблюдается в северной части Москвы, но полагаю, что тенденция будет сходной для всех крупных городов.

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

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

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

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

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

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

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

Что было дальше, вы и сами знаете. Потому, что если вы не полный мудак, у вас возникнет только один вопрос: «АSTER или SRTM?». Отвечаю — ASTER:

Я подготовил растр экспозиций через дефолтный функционал QGIS:

кроме того, сделал растр пересеченности рельефа, каждый пиксел которого представляет собой сумму изменений высот в пределах окна 3х3 пикселя (подробнее смотри в статье Riley S.J., DeGloria S.D., Elliot R. A terrain ruggedness index that quantifies topographic heterogeneities // J. Sci. 1999. V. 5. № 1–4. P. 23–27.). Это не входило в изначальное предположение, но преступлением было бы не проверить возможность взаимосвязи частоты обнажений с индексом пересеченности.

Вокруг каждого описания был построен буфер, радиусом в пол-секунды WGS-84 (приблизительно 30х13 метров), не столько, что-бы облегчить вычисление зональной статистики, сколько нивелировать распиздяйство геологов, которые вначале пишут координаты в пикетажку, а после перебивают их в базу. Кроме того, многие обнажения значительны по простиранию и получение точечной статистики для них явно лишено смысла.

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

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

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

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

(кол-во уступов + 1)/(кол-во заглаженных обнажений + 1)

никак не связана с экспозицией склонов:

Более того, распределение точек наблюдения, в которых не обнаружены выходы коренных пород совершенно аналогично предыдущему:

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

На всякий случай проверим это на растре SRTM:

Как я однажды сказал: «Зрение может обмануть, гистограмма — никогда»:

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

А вот так:

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

Ну а что-же с индексом пересеченности? — спросите вы. Да та-же хуйня. Вот график зависимости количества описанных обнажений от величины индекса пересеченности.

А вот гистограмма этого индекса по всему растру:

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

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