Мультифрактальный анализ данных 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 м (вероятно это связано с квартальной сетью), после чего скорость усложнения падает. Если продолжить это график, он должен опять пойти на спад, что следует хотя-бы из здравого смысла, но ресурсов моего компьютера для таких расчетов уже не хватит. Динамика наблюдается в северной части Москвы, но полагаю, что тенденция будет сходной для всех крупных городов.

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

Взаимосвязь видового разнообразия и проективного покрытия живого напочвенного покрова в результате диссипативного процесса в лесных сообществах Северо-Запада РФ

Тут он мне говорит: «Вы, молодой человек, необразованная охамевшая свинья. Вас пороть надо хорошенько». А я ему отвечаю: «А не пошли бы вы батенька нахуй?». Это выводит его из себя, заставляя морщинистое тело сотрясаться, извергая ругательства: «Ты что себе позволяешь? Ты как со мной разговариваешь? Ах ты ж недоносок, пьянь подзаборная, ты что, совсем охамел!?». А я опять: «Милостивый государь, не ебите мне мозг и потрудитесь проследовать в указанном направлении. То есть — нахуй».

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

Вообще, даже в Шахтах — городе весьма удаленном от мировых столиц, едва ли найдется долбоеб, который не знаком с понятием диссипативной системы. Об этом даже моя собака Лишай знала. Да и хрена ли об этом не знать? Это ведь элементарно. Согласно второму закону термодинамики, любая замкнутая система стремится к полному разпиздосу или, как минимум теоретическому сохранению энтропии. Но в мире, обуреваемом страстями, такую замкнутую систему хуй где найдешь. Даже Новая Азовка с Сидоровкой и то, медленно, но меняются к лучшему. А уж про остальные системы и говорить нечего — какую ни возьми — обязательно открыта.

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

Вот возьмем, к примеру, макароны. Захотели вы заточить пачку спагетти с болгарским кетчупом и жаренной колбасой. Налили воды в кастрюлю и поставили на плиту. Газ пока не зажигаете. Что происходит? Вода, постепенно перестает бултыхать у застаивается без всякого движения как администрация города в Шахтах. Энтропия возрастает, пока не достигает локального максимума. Если кто-то не понимает слово «энтропия», то не расстраивайтесь. «Энтропия» — это «пиздец» по научному. Тишь да гладь, да божья благодать.

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

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

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

А теперь внимание. Ебните свой стакан и слушайте не отвлекаясь.

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

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

Теперь задержите дыхание: я буду кульминировать мою мысль.

Хуй его знает, по какому закону изменяется сложность структуры при изменении потока энергии. Поскольку сие есть тайна пока не раскрытая, возьмем Бритву Оккама и без лишней мозгоебли отрежем все варианты кроме линейного. Что это нам дает? А вот что: при смыкании верхнего яруса изменение структуры напочвенного покрова должно начинаться с резкого уменьшения видового разнообразия и плавного уменьшения проективного покрытия. По мере сгущения крон, скорость уменьшения видового разнообразия снижается, скорость уменьшения проективного покрытия возрастает. Мы получаем две асинхронные волны, по аналогии с синусом и косинусом.

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

Собственно, об этом я уже писал в статье посвященной фрактальному анализу напочвенного покрова несколько лет назад. Но там данных было мало и вообще работа смотрелась несолидно. И вот, наконец-то, я осуществил давнюю мечту и сравнил между собой данные по количеству видов и проективному покрытию травяно-кустарничкового и мохово-лишайникового ярусов, приведенные в монографии В.Н. Федорчук, В.Ю. Нешатаев, М.Л. Кузнецова «Лесные экосистемы северо-западных районов России. Типология, динамика, хозяйственные особенности». Вообще, надо сказать — это охуеннейшая книга. Если где увидите — покупайте обязательно. Даже если наступит полный пиздец и вы не сможете ходить в лес — с помощью одной этой книги вы сможете лет двадцать писать ежегодно новые и актуальные лесоведческие статьи (нахуй они нужны только, но это другой вопрос).

Small_DSCN7736

Итак, у меня было два пакетика травы 75 ампул мескалина 177 геоботанических описаний, каждое из которых сделано с повторностью в 3-10 описаний. Из них я извлек данные по сомкнутости крон, проективному покрытию травяно-кустарничкового, мохово-лишайникового яруса и количество видов в сообщетве. Я сложил покрытия ярусов живого напочвенного покрова. Во-первых, потому что мы не можем учитывать только один произвольно выбранный ярус, а во-вторых, потому что искренне считаю, что всякое необоснованное деление растительности на ярусы есть хуйня на постном масле.

Собственно: данные готовы. Если теория верна, то графики изменения видового разнообразия и проективного покрытия относительно изменения сомкнутости будут асинхронны друг другу. Старик Карл Пирсон вместе с коллегами указывает на наличие обратной связи между суммарным проективным покрытием и видовым разнообразием. И надо сказать весьма заебатой связи: коэффициент корреляции составляет -0,50261 при том, что для преодоления критического значения коэффициента корреляции при уровне значимости 0,01 было достаточно обнаружить тесноту связи |0,15|.

Не будем же тянуть Линнея за яйца. Вот график.

GeoDissipatio1

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

— Нихуя не видно! Че там посредине-то?

— Говно-вопрос. Вот вам центральная вырезка с сомкнутостью от 60 до 80 процентов. Удвоенное количество видов в живом напочвенном покрове показано розовой толстой линией, суммарное проективное покрытие напочвенного покрова показано зеленой.

GeoDissipatio2

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

P.S. Excel-файл с расчетами и графиками.

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

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

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

Дабы загладить боль утраты, я предлагаю вам статью из пропавшего «Лабораторного Журнала» (а где вы ее теперь прочитаете?), описывающую сущность, принципы применимости и алгоритм метода Бенфорда на примере анализа данных о площадях ООПТ России и площадях, охваченных лесными пожарами в 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. Да, я знаю, что качество приведенных картинок отвратительно. Но поверьте, вы встретились с ними в странный момент их жизни.