Gis

Хотелось бы поделиться опытом создания систем маршрутизации PostgreSQL/PgRouting на карте OpenStreetMap. Речь пойдет о разработке [коммерческих] решений со сложными требованиями, для более простых проектов, вероятно, достаточно обратиться к документации. Насколько мне известно, такие вещи, как полная поддержка односторонних дорог и направлений движения, быстрый роутинг на тысячах адресов (порядка секунд на обычном лаптопе, к примеру, Macbook Pro 13" 2013 года), создание дорожного графа с заданными свойствами, мета-оптимизация маршрутов вообще нигде и никак не рассматриваются. Как обычно, все данные и результаты доступны в моем GitHub репозитории OSM Routing Tricks, который я буду пополнять по мере публикаций.

Небольшой маршрут из 330 адресов на карте OpenStreetMap (время построения около 5 секунд на вышеупомянутом лаптопе). Можно ли за это же время построить маршрут, скажем, из 5000 точек? Да, можно, и об этом мы тоже поговорим (в следующих частях статьи).

Что такое маршрутизация и зачем она нужна

Пожалуй, легче сказать, где маршрутизация не используется, чем перечислить все ее применения. Например: маршрутизация нужна для доступа к этой статье в сети интернет, для доставки почты или посылок, для путешествий, и даже для спутниковой интерферометрии, о которой я писал ранее, тоже используется маршрутизация (вычисление суммарного сдвига фазы по замкнутым траекториям на растре)! Далее мы будем говорить только о применении расширения PgRouting для СУБД PostgreSQL и карте OpenStreetMap.

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

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

Версии программ и операционные системы

Для новых проектов рекомендую актуальную версию PgRouting 3.0 и СУБД PostgreSQL 12, хотя тот же код обычно работает с PgRouting 2.6 и СУБД PostgreSQL 9 и более давними. Стоит отметить, что на MacOS (любой версии) и Debian/Ubuntu (любой версии) результаты работы функции PgRouting для поиска оптимального маршрута pgr_TSP() сильно отличаются, при этом, как правило, на линуксах результаты при параметрах по умолчанию получаются значительно лучше, при этом используемый алгоритм "отжига" в компиляции на MacOS работает нестабильно, то есть небольшие изменения параметров приводят к непредсказуемым изменениям результатов при недостаточно продуманном дорожном графе, что, кстати сказать, очень помогает при тестировании. Сам я использую и MacOS и Debian/Ubuntu.

Как загрузить OpenStreetMap в базу данных PostgreSQL

Существует много разных утилит, из которых я предпочитаю GDAL, на мой взгляд, наиболее мощную и предсказуемую. Эта утилита позволяет, в частности, преобразовать дамп OpenStreetMap в дамп PostgreSQL и потом загрузить его как обычный SQL скрипт, см. GDAL: PostgreSQL SQL Dump. Можно загрузить и напрямую в базу, см. документацию GDAL. Пример команды загрузки данных из дампа OSM для Германии (germany-latest.osm.pbf) в базу данных PostgreSQL (osmrouting):

ogr2ogr \    -f PGDUMP \    /vsistdout/ "germany-latest.osm.pbf" \    -nln "osm_lines" \    -nlt LINESTRING \    -progress \    --config PG_USE_COPY YES \    --config GEOMETRY_NAME the_geom \    --config OSM_COMPRESS_NODES YES \    --config OSM_CONFIG_FILE "osmconf.ini" \    -sql "select * from lines where highway is not null" \    -lco FID=id \    | psql osmrouting

Дампы OpenStreetMap по странам предоставляет сервис OpenStreetMap Data Extracts.

Как превратить карту дорог OpenStreetMap в дорожную сеть для роутинга

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

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

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

Для пешеходных маршрутов лучше для каждого добавленного узла связи с домами также добавить узел к парному ребру (обратного направления) и определить короткие сегменты между этими узлами, так как во многих случаях пешеходу удобнее пересекать небольшие улицы, нежели двигаться по одной стороне от перекрестка до перекрестка. Для автотранспортных маршрутов такой вариант тоже подходит, только нужно увеличить расстояние между виртуальными узлами. Скажем, 10 м для пешеходного маршрута и 100 м для автомобильного разумная оценка для каждой смены направления вне перекрестка.

Возможно, вас интересует вопрос, почему мы тратим время на создание сложной дорожной сети вместо использования более продвинутых методов построения оптимального маршрута? Все очень просто по хорошей дорожной сети даже простой алгоритм быстро (секунды) построит отличные маршруты, в то время как по посредственной дорожной сети даже лучшие из существующих алгоритмов за разумное время (минуты, часы или дни, в зависимости от задачи) построят плохие маршруты или вообще не смогут завершить вычисления. Это следствие невероятной вычислительной сложности задачи. К примеру, если мы желаем посетить 10 или 100 адресов по обе стороны одной улицы и при этом у нас определены два направления движения и только два перехода между ними, скажем, в начале и в конце улицы задача имеет единственное решение и любой алгоритм его найдет почти мгновенно! В случае же, когда у нас нет заданных направлений движения и разрешено пересечение улицы около каждого заданного адреса задача становится не решаемой вычислительно уже для нескольких десятков адресов и разные алгоритмы и для разного числа адресов вернут разные и весьма не оптимальные маршруты. Таким образом, критически важно именно построить дорожную сеть с такими ограничениями, которые запрещают большинство (неоптимальных) маршрутов, так что пространство перебора становится несравнимо меньше и достаточно оптимальный результат гарантирован любым из методов. Именно ограничения на направление движения и допустимые повороты являются самыми эффективными.

Таблицы базы данных PostgreSQL

Мы будем работать с базой "osmrouting", содержащей несколько таблиц, содержимое которых доступно в виде PostgreSQL SQL дампов в репозитории, см. также скрипт инициализации базы данных и загрузки дампов load.sh (также загружает необходимые расширения PgRouting и PostGIS).

Первая таблица является единственной необходимой и содержит данные маршрутной сети и, по желанию, информацию для визуализации (сам алгоритм роутинга работает с идентификаторами ребер графа start_id и end_id и некоторой стоимостью каждого сегмента, задаваемой, например, его длиной length):

-- таблица с ребрами маршрутного графа (сегментами дорожной сети)-- id - идентификатор для отладки-- the_geom - геометрия дорожного сегмента для визуализации-- oneway - флаг односторонней дороги-- highway - флаг скоростного шоссе (можно запретить на нем повороты и пересечения)-- pedestrian - флаг пешеходной улицы (можно запретить движение транспорта)-- start_id - идетификатор стартового узла сегмента для роутинга-- end_id - идетификатор конечного узла сегмента для роутинга-- length - длина сегмента (может переопределяться в зависимости от типа дороги и проч.)osmrouting=# \d osm_network                       Table "public.osm_network"   Column   |           Type            | Collation | Nullable | Default ------------+---------------------------+-----------+----------+--------- id         | integer                   |           |          |  the_geom   | geometry(LineString,4326) |           |          |  type       | text                      |           |          |  oneway     | boolean                   |           |          |  highway    | boolean                   |           |          |  pedestrian | boolean                   |           |          |  start_id   | bigint                    |           |          |  end_id     | bigint                    |           |          |  length     | double precision          |           |          | 

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

-- таблица с узлами маршрутного графа-- id - идентификатор для отладки-- the_geom - геометрия дорожного узла для визуализацииosmrouting=# \d osm_nodes                     Table "public.osm_nodes"  Column  |         Type         | Collation | Nullable | Default ----------+----------------------+-----------+----------+--------- id       | bigint               |           |          |  the_geom | geometry(Point,4326) |           |          | 

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

osmrouting=# \d osm_buildings                   Table "public.osm_buildings"  Column  |         Type         | Collation | Nullable | Default ----------+----------------------+-----------+----------+--------- id       | character varying    |           |          |  the_geom | geometry(Point,4326) |           |          | ...

Также тестовые скрипты создают дополнительные таблицы.

Поиск оптимального маршрута

Скрипт репозитория route.sql содержит необходимые команды для выбора случайных 330 адресов домов и построения маршрута по ним с помощью функции PgRouting pgr_TSP(). На самом деле, указанная функция работает не с координатами, а с матрицей расстояний (есть функции и для работы с координатами, см. документацию PgRouting). Матрица расстояний может быть создана функцией pgr_dijkstraCostMatrix(). Заметим, что в скрипте для генерации этой матрицы указан флаг directed=false, так как по умолчанию pgr_TSP() не работает с односторонними дорогами (точнее, не работает с несимметричной матрицей, которая получается при наличии односторонних дорог). С помощью моих функций pgr_dijkstraSymmetrizeCostMatrix() и pgr_dijkstraValidateCostMatrix() это ограничение можно обойти, как мы увидим далее. Что интересно, маршрут возвращается в виде списка идентификаторов узлов дорожного графа (в нашем случае все узлы результата соответствуют добавленным нами виртуальным узлам для домов) и для получения маршрута в виде линии нужно полученный список идентификаторов передать в функцию pgr_dijkstraVia() для нахождения всех посещенных ребер дорожной сети в нужном порядке.

Параметр randomize=false обеспечивает выполнение нескольких запусков алгоритма роутинга и возвращение лучшего результата, для целей отладки и ускорения вычислений можно указать randomize=true, но возвращаемый результат в таком случае непредсказуем. В результате выполнения указанного SQL скрипта создается таблица "route" с сегментами маршрута для каждого заданного адреса, которую можно визуализировать с помощью программы QGIS. Файл проекта QGIS также представлен в репозитории, см. route.qgs Полученная карта с маршрутом представлена на картинке до хабраката.

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

Кроме уже означенной проблемы с направлениями движения, местами можно заметить странный порядок нумерации например, смотрите номера 245,246,247 и другие. Вы можете захотеть покрутить доступные параметры для алгоритма отжига но, даже если это удастся сделать, найденные подходящие параметры для конкретного маршрута не помогут с любым другим маршрутом. Вместо попыток подбора параметров и увеличения времени вычисления следует заняться улучшением дорожной сети (и матрицы расстояний), например, за счет жесткого указания направлений движения.

Заключение

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

Продолжаем цикл статей про построение систем роутинга со сложными требованиями на основе Open Source базы данных PostgreSQL и расширения PgRouting на карте OpenStreetMap. Сегодня мы поговорим о том, как добавить поддержку односторонних дорог (направлений движения). Зачастую, именно отсутствие этой поддержки является основной причиной смены PgRouting на другой "движок" маршрутизации. Как обычно, все данные и результаты доступны в моем GitHub репозитории OSM Routing Tricks, который я пополняю по мере публикаций.

Небольшой маршрут из 330 адресов на карте OpenStreetMap.

Что такое односторонние дороги и зачем они нужны

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

Добавляем поддержку односторонних дорог в PgRouting

Официальная документация PgRouting для функции pgr_TSP Using Simulated Annealing approximation algorithm сообщает нам:

If using directed := true, the resulting non symmetric matrix must be converted to symmetric by fixing the non symmetric values according to your application needs.

Отлично, но в документации нет ни слова о том, как это сделать. Нам придется начать с теории и разобраться, возможно ли это и как именно это сделать. Англоязычная страница википедии, посвященная проблеме коммивояжера, содержит нужный нам раздел Travelling salesman problem: Asymmetric:

Solving an asymmetric TSP graph can be somewhat complex. The following is a 33 matrix containing all possible path weights between the nodes A, B and C. One option is to turn an asymmetric matrix of size N into a symmetric matrix of size 2N.

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

Пусть у нас задана асимметричная матрица весов:

Ей соответствует следующая симметричная матрица весов:

где -w это виртуальные соединения для виртуальных узлов, которые рекомендуется задать некоторым отрицательным числом, потому что

w=0 is not always low enough

Теперь мы можем написать соответствующую функцию-обертку для PgRouting. Поскольку PgRouting отбрасывает все отрицательные значения в матрице весов, мы используем нулевое значение (мы это компенсируем заданием очень высоких весов для запрещенных направлений) и добавляем еще [избыточные] ребра с очень высоким весом, поскольку PgRouting требует указания весов между всеми узлами матрицы (в теории этого не требуется, это просто техническое ограничение, которому нам приходится соответствовать).

CREATE OR REPLACE FUNCTION pgr_dijkstraSymmetrizeCostMatrix(matrix_cell_sql text,    OUT start_vid BIGINT, OUT end_vid BIGINT, OUT agg_cost float)RETURNS SETOF RECORD AS$BODY$DECLARE    sql text;    r record;BEGIN    sql := 'with edges as (' || matrix_cell_sql || ')        select 3e9+start_vid as start_vid, end_vid as end_vid, agg_cost from edges        union        select end_vid, 3e9+start_vid, agg_cost from edges        union        select 3e9+start_vid, start_vid, 0 from edges        union        select start_vid, 3e9+start_vid, 0 from edges        union        select start_vid, end_vid, 1e6 from edges        union        select 3e9+start_vid, 3e9+end_vid, 1e6 from edges        ';    FOR r IN EXECUTE sql LOOP        start_vid := r.start_vid;        end_vid   := r.end_vid;        agg_cost  := r.agg_cost;        RETURN NEXT;    END LOOP;END;$BODY$LANGUAGE plpgsql VOLATILE STRICT;

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

An Infinity value was found on the Matrix

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

CREATE OR REPLACE FUNCTION pgr_dijkstraValidateCostMatrix(matrix_cell_sql text,    OUT start_vid BIGINT, OUT end_vid BIGINT, OUT agg_cost float)RETURNS SETOF RECORD AS$BODY$DECLARE    sql text;    r record;BEGIN    sql := 'WITH RECURSIVE src AS (' || matrix_cell_sql || '),    dst AS (        select            *        from src where start_vid =            (select                start_vid            from src            group by start_vid            order by count(*) desc            limit 1)        union        select            src.*        from src, dst        where src.start_vid=dst.end_vid    )    select * from dst';    FOR r IN EXECUTE sql LOOP        start_vid := r.start_vid;        end_vid   := r.end_vid;        agg_cost  := r.agg_cost;        RETURN NEXT;    END LOOP;END;$BODY$LANGUAGE plpgsql VOLATILE STRICT;

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

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

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

Поиск оптимального маршрута для автомобиля

Поскольку для пешехода тротуары доступны в любом направлении, а для автомобиля одностонние дороги доступны лишь в одном направлении, необходимо указывать разные весовые функции для роутинга пешеходного и автомобильного. Как уже упоминалось ранее, наша дорожная сеть оптимизирована для автомобильного роутинга, о нем и поговорим. Итак, запретим (укажем очень большую длину миллион метров) движение в противоположном направлении. Ранее при подготовке дорожной сети мы разделили каждую дорогу на две, одна из которых (type='osm') совпадает с исходной и вторая (type='osmreverse') инвертирована. Соответственно, для односторонней дороги добавленная нами ее инвертированная копия должна быть запрещена для любого движения (вообще говоря, можно ее и вовсе не создавать, но тогда будет немного сложнее объяснить построение дорожной сети). Кроме того, для прямого движения должна быть доступна только исходная дорога (type='osm') и для обратного инвертированная (type='osmreverse'). В таком случае, итоговые прямая и обратная весовые функции имеют вид:

    case when (type='osmreverse' and oneway) then 1000000 else length end as cost,    case when type ilike 'osm%' then 1000000 else length end as reverse_cost,

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

Построим оптимальный маршрут все для тех же случайных 330 адресов домов и с помощью функции PgRouting pgr_TSP() и добавленных выше функций pgr_dijkstraSymmetrizeCostMatrix() и pgr_dijkstraValidateCostMatrix(). Теперь мы можем использовать флаг directed=true, так как теперь мы добавили для pgr_TSP() поддержку односторонних дорог (точнее, симметризацию и заполнение пропущенных значений в несимметричной матрице, которая получается при наличии односторонних дорог). Как и ранее, в результате выполнения немного модифицированного SQL скрипта route.sql создается таблица "route" с сегментами маршрута для каждого заданного адреса, которую можно визуализировать с помощью программы QGIS. Файл проекта QGIS тоже обновлен, см. в репозитории route.qgs Полученная карта с маршрутом представлена на картинке до хабраката.

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

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

Как и указано на карте OpenStreetMap, в построенном маршруте движение выполняется слева направо для участка дороги с пунктами маршрута 329,330,331 на левой стороне дороги:

и в том же направлении (слева направо) для участка дороги с пунктами маршрута 72,73,74 на правой стороне дороги (второй проход по этому участку дороги):

Нумерация домов по сторонам улиц между перекрестками последовательная. Зачастую нумерация последовательная и после перекрестков (см. картинки выше). Таким образом, мы добились заметного улучшения качества маршрута, не меняя дорожную сеть и не занимаясь оптимизацией параметров функции pgr_TSP().

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

Заключение

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

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

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

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

Видите взаимосвязь ортофотоснимка и рельефа? Если да, то вы или геолог или можете им стать: корреляция компонентов (разложения в пространственный спектр) составляет 41% для длины волны 20 м, 58% для 50 м и 99% для 300 м (Jupyter Python ноутбук с вычислениями доступен по ссылкам ниже). Большинство геофизиков поклянутся, что это у вас спектры порченые (записано с натуры), игнорируя и геофизику и прилагаемые вычисления и ссылки на публикации.

Миф первый, рельеф и сила тяжести

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

Важнейшим доказательством изостазии является отсутствие связи между рельефом и силой тяжести.

Увы, именно в отсутствие такой связи верят многие из встреченных мной геофизиков (но не геологи), хотя это абсолютно неверно! Что интересно, в англоязычной версии этой же викистраницы написано прямо противоположное, смотрите пример с айсбергом в разделе Deposition and erosion:

An analogy may be made with an iceberg, which always floats with a certain proportion of its mass below the surface of the water. If snow falls to the top of the iceberg, the iceberg will sink lower in the water. If a layer of ice melts off the top of the iceberg, the remaining iceberg will rise. Similarly, Earth's lithosphere "floats" in the asthenosphere.

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

Посмотрим, что говорит оглавление русскоязычного учебника по геофизике (Викулин, 2004):

Глобальные волны геоида, отсутствие их связи с особенностями строения земной коры

Здесь глобальные волны геоида это изменения (аномалии) поля силы тяжести порядка масштабов планеты, для которых, как и указано, отсутствует связь с рельефом. В то же время, для всех компонентов силы тяжести меньшего масштаба эта связь есть! Могу лишь предположить, что студенты-геофизики зачастую не понимают смысла выражения "глобальные волны геоида" (кстати, для обработки спутниковых данных современные глобальные модели оперируют "изменяющимся во времени квазигеоидом" с явным указанием порядка волны...).

Более того, корреляция компонентов силы тяжести (или вертикальной компоненты гравитационного поля) и рельефа близка к 100%, как показано в работе "Radially symmetric coherence between satellite gravity and multibeam bathymetry grids" (Marks and Smith, 2012). На графике ниже график радиальной когеренции (справа) из указанной статьи дополнен мною двумерной коррелограммой (слева) для глобальной модели гравитационных аномалий Sandwell and Smith Gravity Anomaly и глобальной модели рельефа GEBCO 2019:

Если на графике радиальной когеренции (корреляции) не понятно, почему для первых километров масштаба связь компонентов мала, то на двумерной коррелограмме все становится ясно измеренная сила тяжести относится к поверхности океана, а топография, очевидно, ко дну. Таким образом, здесь мы видим, что до масштабов ~40 км корреляция компонент близка к 100%. Более того, для региона Индонезия высокая корреляция (>75%) наблюдается до масштабов в сотни километров и средняя (>50%) для масштабов в тысячи километров, как показано в статьях по ссылкам ниже.

Миф второй, рельеф и аномалия Буге

Что интересно, миф второй полностью противоречит первому, что не мешает многим геофизикам верить в оба разом. Встречайте Аномалию Буге (Бугера):

Аномалия Бугера расхождение между теоретически вычисленным и реальным значением гравитационного поля Земли в определённой точке референц-эллипсоида.

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

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

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

И следующий шаг, или редукция:

Уточнённая, или полная редукция Бугера позволяет полностью учесть влияние рельефа местности.

Здесь речь о том, что наличие на местности неровностей в виде возвышений или углублений влияет на измеренную в этих местах силу тяжести, так что для сравнения измеренных значений нужно и их компенсировать. Заметим, что такая поправка, то есть необходимость учитывать влияние рельефа местности на силу тяжести, очевидно противоречит первому мифу. И эта поправка Буге проста в использовании нужно было среднюю плотность пород умножить на высоту точки наблюдения (возможно, Буге и до интеграла по глубине бы додумался, только где ему было взять среднюю плотность пород на территории в зависимости от глубины). Для самого Буге, с учетом расстояния между точками наблюдения километров 20 и более, о большой точности речь не шла, и средняя плотность собранных им же геологических образцов по региону его вполне устраивала для вычисления поправки. Увы, для современных детальных измерений этот метод совсем никуда не годится.

Посмотрим практический пример, как компоненты силы тяжести глобальной модели WGM2012 после редукции Буге (слева) и еще одной редукции (справа) коррелируют с компонентами рельефа глобальной модели GEBCO 2019:

Как видим, про "полностью учесть рельеф местности" речь тут явно не идет то есть результат редукции Буге противоречит его определению, а значит, его использование неправомерно. Обратим внимание на правый график тут мы ожидаем максимум корреляции 100%, а наблюдаем всего около 75%, что показывает неточность самой модели WGM2012 (по сравнению с использованной выше Sandwell and Smith Gravity Anomaly, авторы которой вовсе не вычисляют редукцию Буге по уже понятным причинам), таким образом, здесь реальная корреляция редукции Буге с рельефом получится около -100%.

Если обратиться к учебнику "Гравиразведка" (Утёмов, 2009), там в соответствующем разделе есть разъяснение полученного результата:

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

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

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

Миф третий, сферический конь эллипсоид в вакууме

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

В геофизике модель геопотенциала представляет собой теоретический анализ измерения и расчета эффекты Земли гравитационного поля

Да-а, расчета эффекты Земли это, конечно, не совсем русский язык, или даже совсем не русский. Равно как и теоретический анализ измерения тоже тот еще шедевр. Тем не менее, можно понять, что модель геопотенциала это модель (глобального) гравитационного поля Земли. Это определение уже ошибочное на самом деле, рассматривается вертикальная компонента гравитационного поля, то есть сила тяжести (гравитационное поле векторное, а сила тяжести скалярная, и связь между ними весьма сложная). Но, главное, геофизикам такая модель вообще не нужна! А кому нужна, описано дальше в этой же викистатье:

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

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

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

Миф четвертый, о пространственных спектрах

У меня в коллекции есть такие перлы, как спектры у вас порченые (освятить, что ли, хотят?) и другие. Попросту говоря, геофизики не имеют никакого понятия, что это такое.

Давайте посмотрим, как же практически работать с пространственными спектрами. Как всегда в математике, все эти преобразования по сути чрезвычайно просты. Например, результат применения гауссова фильтра масштаба N метров есть низкочастотная фильтрация (low-pass filter), разница между исходным изображением и фильтрованным есть высокочастотная фильтрация (high-pass filter) и разница между двумя фильтрами масштаба N и M метров есть полосовая фильтрация (band-pass filter). Для выделения компонента пространственного спектра шириной l метров для длины волны L метров достаточно выполнить полосовую фильтрацию с фильтрами масштаба L-l/2, L+l/2. Спектр мощности может быть представлен как стандартное отклонение полученных компонент, а через логарифм спектра мощности посчитана и фрактальная размерность. Для примера посмотрите график радиальной когерентности из статьи НАСА и двумерную спектрограмму, полученную указанным методом, в обсуждении первого мифа выше.

Низкочастотная фильтрация это улучшенная замена Буге преобразования:

А если выполнить высокочастотную фильтрацию, то мы сразу получаем все глобальные разломы то есть границы литосферных плит и микроплит:

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

Интересно, что аргентинские геофизики применяют индекс фрактальности для вычисления средней плотности по территории в зависимости от глубины для использования в преобразовании Буге. В самом деле, этот способ работает для территории Латинской Америки (и не только там, но не везде по миру все зависит от характера профиля плотности) я проверил с помощью численного моделирования, по ссылкам ниже вы найдете готовый Jupyter Python ноутбук с вычислениями. А вот глобальные модели гравики для Латинской Америки с помощью сферических разложений получаются очень так себе (впрочем, для космических вычислений это не мешает, там не нужны детальные пространственные компоненты силы тяжести на уровне поверхности) слишком высока фрактальность территории (сильно расчлененные горы с их отрогами), видимо, потому и пришлось придумывать фрактальный метод анализа.

Заключение

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

Бонус

И специальное дополнение для мифоборцев. Не стал добавлять еще один миф, просто "на закуску" приведу выдержку из рецензии, когда-то полученной от журнала "Отечественная геология":

Авторы пытаются нас убедить в том, что математические манипуляции с оптическим сигналом одного из каналов мультиспектральной съемки позволяют нам анализировать глубинное строение Земли. Математические обработки мультиспектральных снимков направлены на улучшение их визуального восприятия и классифицирования представленных на них объектов для создания тематических карт.

Так что учтите если вы ссылаетесь на результаты из публикаций, скажем, НАСА (аналогично корреляции между пространственными спектрами рельефа и гравиметрии посчитав таковую между рельефом и космическими снимками) то будете обвинены в ереси ("математических манипуляциях") и все ваши "порченые спектры" придется немедленно освятить (боюсь, путем окунания компьютера в чан со святой водой и никак иначе). Я предупредил! :)

Ссылки

Marks, K.M., Smith, W.H.F., 2012. Radially symmetric coherence between satellite gravity and multibeam bathymetry grids. Mar Geophys Res 33, 223227. https://doi.org/10.1007/s11001-012-9157-1

Викулин А.В., 2004. ВВЕДЕНИЕ В ФИЗИКУ ЗЕМЛИ. УЧЕБНОЕ ПОСОБИЕ ДЛЯ ГЕОФИЗИЧЕСКИХ СПЕЦИАЛЬНОСТЕЙ ВУЗОВ

Утёмов Э.В., 2009. Гравиразведка

Spectral Coherence between Gravity and Bathymetry Grids

WGM2012 spatial components of free-air gravity and topography correlation

Bouguer and Free-Air Gravity Anomalies in terms of spatial spectrum

There is a high correlation between DEM and ortho photos

3D Density Inversion by Circular Hough Transform [Focal Average] and Fractality Index

[Gaussian Filtering on Spheroid [Sandwell & Smith]](

Ниже я картограф Александр Пронин расскажу, как совместно с моим коллегой Никитой Славиным была создана карта возраста домов города Владимира.

Данные

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

Первый мой сайт владимирдом.рф карта застройки. Ее я создал для того, чтобы показать жителям, как застраивался наш город. Для нее данные о годе постройки брал с сайта Мин.ЖКХ, для детских садов и школ с их официальных сайтов. Много данных по году прислали пользователи карты. На карте есть данные только о годе постройки.

Второй кадастровые данные. Меня интересуют объекты капитального строительства. У многих указан год постройки и адрес.

Третий сайт открытых данных Министерства культуры. Нужны объекты культурного наследия. С него собираю данные по названию объекта, адресу и фотографию.

Четвертый сайт wikimapia.org. Есть название объекта, адрес и фотография. После сбора данных приступаем к их обработке.

Геопроцессинг

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

Слой с OpenStreetMap

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

Слой с OpenStreetMap + данные с владимирдом.рф

Слой с OpenStreetMap + данные кадастра

Слой с OpenStreetMap + открытые данные Министерства культуры

Слой с OpenStreetMap + данные с wilimapia.org

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

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

• Название
• Адрес
• Год постройки
• Этажность
• Фото
• Источник фото

Всего на карте получил 21429 зданий, у 8017 год постройки известен.

Красные полигоны год неизвестен, зеленые известен

После того, как построил график, в глаза бросается большое количество домов, построенных в 1917 году.

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

Палитра цветов

Заниматься визуализацией данных буду в QGIS. Палитра цветов, которая зависит от года постройки, будет представлена в непрерывной палитре от горячих цветов к холодным. Беру темную базовую карту от CartoDB и цветовую палитру Spectral, предустановленную в QGIS.

Для красного и синего цветов делаю фотографии здания банка на Соборной площади и панельного дома по улице Василисина. Остальные фотографии беру у фотографа Ивана Медведева. Оранжевый цвет водонапорная башня, зеленый кукольный театр, жёлтый дом 3 по улице Луначарского (бывш. Совнархоз), голубой Городской дворец культуры и серый Политехнический колледж.

Тут же в QGIS делаю базовую карту. На ней будет 3 слоя: сплошной цвет, дороги и водные объекты.

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

Публикация карты

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

Любой пользователь может нажать на кнопку Редактировать и предложить свои правки, а так же загрузить к зданию его фотографию.

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

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


Человейники Гонконга. Интересно, можно ли будет сделать такую фотку в Москве лет через 30?

О проекте

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


Типичный человейник заглядывает в окна кирпичного дома. Автор фотографии Михаил Мельников.

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


Цель проекта ХоумХаб - визуализация качества городской среды, подкрепленная объективными цифрами и понятными расчетами. Нам уже есть, что показать: на карте представлены данные по плотности застройки и экологической ситуации, аналогов которым в открытом доступе найти не удалось.

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


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

Из чего складывается качество жизни в мегаполисе

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

• Уровень шума.
• Социальная инфраструктура.
• Транспортная инфраструктура.

Зато у нас уже готов расклад по двум ключевым факторам:

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

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

Плотность застройки

Плотность застройки характеризует интенсивность использования территории. Слишком плотная застройка таит ряд опасностей:

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


Хватит ли пропускной способности этой дороги, когда достроятся все 45 корпусов? Кстати, как на единственной полосе разъезжаться?

Нехватка солнечного света. Существует специальный ГОСТ о продолжительности инсоляции (естественного освещения) жилых помещений. Застройка окна-в-окна его игнорирует чуть более чем полностью:


Солнечные лучи попадают в этот двор и окна максимум на 40 минут в день. Фотография взята из прекрасного поста про питерский колодец-могильник.

Разрушение природной экосистемы. Сравните просторные зеленые дворы между кирпичными домами советской постройки и заставленные машинами бетонные каньоны новостроек. Березы или яблоневый садик под окнами? Забудьте.

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

Огромная плотность населения безупречные условия для распространения пандемии. COVID-19 непрозрачно намекнул, что ряды высоченных человейников с одним супермаркетом на весь микрорайон так себе идея.

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




Коробки в боевой раскраске. Обе фотографии взяты из блога фотографа, столкнувшегося с человейниками.

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

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

Расчет плотности застройки

Плотность застройки города принято оценивать по двум параметрам:

Процент застройки процентное соотношение площади зданий к общей площади участка, на котором они стоят.

Процент плотности застройки процентное соотношение суммарной поэтажной площади зданий на участке к его площади.

В разных странах и городах исторически закрепились свои нюансы расчетов. Например, в Нью-Йорке площадь здания отмеряется от внутренних стен, а площадь участка берется строго по его границам. В Лондоне площадь здания меряется от наружных стен, а площадь участка захватывает половины улиц, окружающих участок. Во многих странах принято измерять плотность застройки в расчете на гектар земли (1 га = 10 000 м), а в наших мегаполисах прижилось измерение в расчете на квартал. Площади кварталов могут отличаться в разы, и это привносит неразбериху.

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

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

Будем расчитывать процент плотности застройки и процент застройки вот для этого дома:


Москва, Ломоносовский проспект, 14.

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

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


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

Почему конутры участка на скриншоте получились такой странной формы? Они расчитаны с помощью функции PostGIS ST_Buffer(geometry, radius_of_buffer), формирующей полилинию точек, отстоящих от заданной геометрии на расстояние radius_of_buffer. Параметр radius_of_buffer равен 200 м, а сглаживание углов включено по умолчанию. Для наглядности вот примеры результатов ST_Buffer() на этом и других зданиях, но с radius_of_buffer всего в 10 м:


Контуры домов фиолетовые. Контуры участка, полученные с помощью ST_Buffer() синие.

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

Теперь находим все попадающие на участок водоемы. Вычитаем их площадь пересечения с участком из общей площади участка. Мотивация простая: требуется определить площадь участка вокруг дома, по которой можно спокойно перемещаться. Грубо говоря, у дома на скале, вокруг которого океан, застройка примерно 100%.

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

Для некоторых дорог в OSM указан тег lanes, обозначающий количество полос. Для обозначения значимости дороги используется тег highway. По тому, является ли дорога автомагистралью федерального значения, городской трассой или внутриквартальным проездом, можно предположить ее полосность. Ширину полосы определяет ГОСТ: она зависит от типа трассы и максимально допустимой скорости (OSM-тег maxspeed). Например, для highway=motorway ширина полосы может быть 3.75 м, а для дороги во дворе 2.75 м.


Дороги внутри участка отмечены темно-серым, трамвайные пути светло-серым.

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

Процент застройки дома 17%:


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

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

Процент плотности застройки дома 136%:


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

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


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

С рекордсменом по проценту плотности застройки все понятно. А где самый высокий процент застройки? Конечно же, в центре города. Во многих дворах внутри Садового кольца он достигает 85%. Но для исторического центра это нормально: многие европейские города могут похвастать 95% застройки.


Приз в номинации Максимальный процент застройки в Москве вручается дворам возле Тверской улицы.

А вот как выглядит на карте уплотнительная застройка дворов возле метро Домодедовская:


Слева блочки с адекватной плотностью застройки, справа новый ЖК бизнес-класса, застроенный окна-в-окна. Ибо сказано: за слово бизнес в названии и входную группу с искуственным мрамором не стыдно доплатить по тарифу x2.

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

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

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

Поэтому мы реализовали более наглядный способ генерализации:

1) Замостить карту гексагонами с длиной ребра 100 метров для 12 уровня зума, 200 м для 11 зума и так далее вплоть до 6 зума. Для этого была использована функция PostGIS ST_HexagonGrid(hex_edge_len, bounds).

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

3) Залить гексагоны соответствующими цветами. В результате генерализация будет выглядеть примерно так:


Гексагональная генерализация процента застройки: средние уровни зума, на которых виден центр города.

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

Обратная сторона ХоумХаба

Настало время поговорить про технологии, которые стоят за нашей картой.

Наши источники данных

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

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

Для карты рельефа мы используем DEM (digital elevation model) SRTM v3. Это данные, полученные радарной топографической съемкой материковых и островных территорий с 54 южной широты по 60 северной широты.

SRTM считается одним из самых качественных открытых датасетов, и его точности (1 измерние высоты примерно на 50 метров поверхности) должно хватить для адекватного использования в моделях экологии и распространения шума.

Пайплайн генерации тайлов

Наша цепочка python-скриптов для скачивания данных, импорта в PostgreSQL, подсчета данных для слоев и нарезки их на тайлы разбита на несколько docker-образов и запускается через docker-compose. Весь процесс выглядит так:

Подготовка растровых тайлов с рельефом местности. Датасет SRTM состоит из набора растровых тайлов 1 x 1, которые мы пропускаем через последовательность GDAL-утилит для генерации прозрачных тайлов с хилшейдингом, тайлов с цветовым рельефом и совмещения этих двух наборов в комбинированные тайлы рельефа местности.

Подготовка векторных тайлов с данными из OSM, слоями застройки и экологии. С ресурса Geofabrik скачивается фрагмент OSM, а конкретно Центральный федеральный округ России. Из него с помощью osmium-extract извлекается окружность, включающая Москву и некоторую часть МО.

Для заливки данных в PostgreSQL запускается импортер Imposm3. Он хорош тем, что позволяет удобно сконфигурировать параметры генерализации для лесов, озер, дорог и произвольных объектов.

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

Финальный аккорд запуск генерации векторных тайлов. Для этого используется самописный скрипт с конфигом, конкретизирующим, какие зумы и с какими данными генерировать. В PostGIS есть все необходимые для этого функции, такие как ST_MakeEnvelope(), ST_Intersects(), ST_AsMVT(), ST_AsMvtGeom().

Вот карта и готова. Для отрисовки в браузере используется опенсорсная библиотека Mapbox GL JS. Стиль карты был запилен по-простому, безо всяких WYSIWYG-примочек, в виде гигантского json в формате Mapbox.

О планах и о команде

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

Архитектор проекта, фулстек-разработчик, девопсер, дизайнер, веб-мастер и копирайтер-на-пол-ставки это все я. В душе я бэкендер, поэтому вместо мобильной версии сайта на js у меня получилась консольная тулза на C++. Пришлось переделывать. Вышел микросервис. Разработка забуксовала. К счастью, в проект ворвался Rostixman: он за одно воскресенье запилил вменяемую мобильную версию и пофиксил косяки интерфейса, за что ему отдельное спасибо.

ХоумХаб развивается, задач становится все больше и они все заковыристее. Поэтому наша команда упоротых исследователей мегаполиса ищет новых разработчиков, особенно фронтендеров! Будем рады разработчикам, которым близка идея проекта и которые готовы посвящать ему свое время.

Мы планируем добавить на карту новые важные факторы, влияющие на качество городской среды:

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

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

Транспортная инфраструктура. Боль многих городов-сателлитов в вечной пробке на единственной дороге между сателлитом и мегаполисом. А у многих микрорайонов проблемы с доступностью на общественном транспорте. Поэтому мы планируем сделать два слоя: для автомобилистов и для приверженцев общественного транспорта.

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

Кроме того, мы собираемся охватить как можно больше городов, в идеале всю Россию. Было бы интересно посравнивать наши мегаполисы с Гонконгом, Нью-Йорком и прочими знаменитыми рассадниками человейников. А исторический центр сравить например с центром Праги или Рима.

А что за карта без гео-поиска? Мы сделаем и его. Но конечно же приоритет 1 улучшать, валидировать и снова улучшать наши модели.

Заключение

ХоумХаб это карта со слоями, которые помогают находить ответы на наболевшие вопросы. Какая плотность застройки в новом микрорайоне? Обнаружены ли в нем зачатки инфраструктуры и можно ли из него выбраться своим ходом? Сколько заводов и промзон в Капотне? Стоит ли ценителю тишины покупать квартиру с окнами на Ленинградское шоссе? Как долго идти от ЖК Грин Парк до ближайшего парка?

А на какие вопросы ищете ответа вы? Делитесь мнением, реквестуйте фичи. Критикуйте, это же Хабр в конце концов! А мы в предстоящих релизах будем приоткрывать новые грани жизни в мегаполисе.


Под занавес: это не Гонконг, это Марфино. А вы бы хотели, чтобы наши мегаполисы превратились в сплошной человейник?