Функции управления цифровыми активами автомобильных дорог. Часть 1 – сегментация

Здравствуйте, уважаемые читатели Хабра!

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

Содержание

Введение
Алгоритм сегментации дороги
Итоговое решение
Примеры сегментации различных по геометрии дорог
Заключение
Ссылки, источники, материалы

Введение

Мы, команда «Инновации-ИТ», занимаемся построением высоконагруженных информационных систем. Помогаем клиентам быстро находить верные ответы благодаря системам поддержки принятия решений, бизнес-аналитике и построению отчётности. Работаем над проектами в сферах Smart City и Smart Transport.

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

**Исходный цельный полигон дороги (увеличенное изображение)**

На рисунке ниже представлен сегментированный полигон. Каждый сегмент имеет уникальный составной идентификатор (layer_id, object_id, zone_id). Также в колонках area и len выведены площадь (ST_Area(geom)) и длинна (ST_MaxDistance(geom)) каждого сегмента.

**Сегментированный исходный полигон дороги**

**Сегментированный исходный полигон дороги (увеличенное изображение)**

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

Зачем нужно автоматизировать данную операцию? Суммарная протяженность дорог в городе составляет несколько тысяч километров. Тем не менее специалистам однократно героически в полном объёме удалось вручную произвести сегментацию полигонов. Однако, прокладываются новые дороги, появляются развязки и мосты. Актуальность данных в системе напрямую зависит от того насколько оперативно специалисты смогут обработать новые полигоны. Заявка на сегментацию может прийти в любой момент. Данная задача требует постоянного мониторинга, но по своей сути она проста и рутинна. К тому же тратится время ценнейших специалистов, у которых есть не менее интересные и важные задачи.

Алгоритм сегментации дороги

Что должна обеспечивать функция:

Длинна сегмента — 300 метров. Такое значение было выбрано эмпирически, в будущем возможно другое. Это оптимальная длина, хорошо отображаемая на используемых масштабах карты. Разумеется, не всегда все участки будут одинаковой длины — общая протяжённость дороги может не быть кратной 300 метрам. Поэтому данное значение является приблизительным. Дорога может не делиться нацело на это число. В результате могут появляться участки с длиной, например, 50 метров. Такие малые сегменты должны присоединяться к соседним, более крупным. Каждому итоговому сегменту присваивается уникальный идентификатор.
Учитывать административные границы округов и районов. Если дорога пересекает границу округа или района, в точке пересечения должен быть принудительный разрыв сегмента.
Быстрое выполнение. Время выполнения функции для дороги протяженностью до 100 км не превышает установленного лимита 5-10 секунд.

Реализуем в виде SQL-функции:

Входные параметры:
- layer_id — идентификатор слоя,
- object_id — идентификатор объекта дороги в слое.

Выходные данные: Запись (INSERT/UPDATE) в целевую таблицу активов с заполненными полями
- layer_id — идентификатор слоя,
- object_id — идентификатор объекта дороги в слое,
- zone_id — идентификатор сегмента,
- geom — геометрия сегмента.

В качестве примера рассмотрим работу алгоритма для объекта дороги как на рисунке выше: layer_id = 179, object_id = ‘10000257’.

**Исходный полигон дороги для разрезки**

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

1) Заносим каждый part_id дороги для последующего разбиения в таблицу source_road_solid_parts. В таблице представлены колонки: part_id; layer_id — идентификатор слоя дороги; object_id – идентификатор объекта в слое; geometry – геометрия дороги. Всё как на рисунке выше.

Если линейные размеры полигона малы для его разбиения (в нашем случае ST_MaxDistance(geom) <= 300 ), то он сразу попадает в результирующую таблицу с уникальным zone_id. Таким образом, все последующие шаги для него пропускаются.

2) Получение скелета дороги. Функция generate_road_skeleton строит скелет дороги и добавляет его в таблицу road_skeleton. В начале упрощаем и производим буферизацию полигона, тем самым сглаживаем границы и избавляемся от вырезов внутри него: клумб, автобусных остановок, люков и прочего ( ST_Buffer, ST_SimplifyPreserveTopology, ST_MakePolygon(ST_ExteriorRing(geometry))). После чего вызываем ST_StraightSkeleton или ST_ApproximateMedialAxis для получения скелета. Если ST_StraightSkeleton не справилась – результат NULL, то вместо неё применяем ST_ApproximateMedialAxis.

**Скелет полигона (увеличенное изображение)**

**Полигон дороги и скелет (увеличенное изображение)**

Функция generate_road_skeleton

-- функция для создания скелета дорогиCREATE OR REPLACE FUNCTION road_processing.generate_road_skeleton(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,    p_tolerance FLOAT default 2.0)RETURNS VOID AS $$begin-- удаляем старые записи из таблицы road_skeletondelete from road_processing.road_skeletonwhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;with…-- Создаем сплошной полигон без дырокfilled_buffered_united_road AS (    SELECT ST_MakePolygon(ST_ExteriorRing(geometry)) as geometry    FROM buffered_united_road),-- Упрощение геометрии для построения скелетаsimplified_filled_buffered_united_road AS (    SELECT ST_SimplifyPreserveTopology(geometry, p_tolerance) as geometry    FROM filled_buffered_united_road),-- Строим скелетskeleton AS (SELECT COALESCE(ST_StraightSkeleton(geometry), ST_ApproximateMedialAxis(geometry)) as geometry--SELECT ST_ApproximateMedialAxis(geometry) as geometry--    SELECT ST_StraightSkeleton(geometry) as geometry-- SELECT CG_StraightSkeletonPartition(geometry, true) as geometry    FROM simplified_filled_buffered_united_road)-- добавляем скелет в таблицуinsert into road_processing.road_skeleton (part_id, layer_id, object_id, geom)select p_part_id,p_layer_id,p_object_id,geometryfrom skeleton;END;$$ LANGUAGE plpgsql;

3) Разделяем скелет на составляющие его линии (LINESTRING). Функция fill_road_segments разделяет скелет полигона на составляющие его линии и заносит их в таблицу road_segments.

Функция fill_road_segments

CREATE OR REPLACE FUNCTION road_processing.fill_road_segments(    p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$BEGIN-- удаляем старые записи из таблицы road_segmentsdelete from road_processing.road_segmentswhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;    INSERT INTO road_processing.road_segments (part_id, layer_id, object_id, id, geom, length)-- получаем линии    WITH dumped AS (        SELECT (ST_Dump(geom)).geom AS geom        FROM road_processing.road_skeleton        WHERE part_id = p_part_id AND  layer_id = p_layer_id AND object_id = p_object_id     ),-- создаём коллекцию (список) из линий    collected AS (        SELECT ST_Collect(geom) AS geom        FROM dumped    ),-- разделяем линии по точкам пересечения,    noded AS (        SELECT ST_Node(geom) AS geom        FROM collected    ),    geometries AS (        SELECT             p_part_id as part_id,             p_layer_id as layer_id,             p_object_id as object_id,             (ST_Dump(geom)).geom AS geom        FROM noded    )    SELECT         part_id,         layer_id,         object_id,        row_number() OVER () as id,        geom,        ST_Length(geom) as length    FROM geometries;END;$$ LANGUAGE plpgsql;

4) Формируем вершины, соединяющие линии скелета. Заполняем таблицу road_vertices, хранящую связи между концами отрезков из таблицы road_segments, применяем функцию pgr_extractVertices. В таблице: id – идентификатор вершины; in_edges, out_edges – списки идентификаторов входных и выходных отрезков (рёбер) из road_segments; geom – геометрия вершины (точка); x, y – координаты вершины.

Функция fill_road_vertices

CREATE OR REPLACE FUNCTION road_processing.fill_road_vertices(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$begin-- удаляем старые записи из таблицы road_verticesdelete from road_processing.road_verticeswhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;INSERT INTO road_processing.road_vertices (part_id, layer_id, object_id, id, in_edges, out_edges, x, y, geom)SELECT    p_part_id,p_layer_id,p_object_id,    id,    in_edges,    out_edges,    x,    y,    geomFROM pgr_extractVertices(format('SELECT id, geom FROM road_processing.road_segments WHERE part_id = %s AND layer_id = %s AND object_id = %L', p_part_id, p_layer_id, p_object_id));END;$$ LANGUAGE plpgsql;

5) Сопоставляем каждому сегменту его вершины. Функция update_road_segments обновляет таблицу road_segments на основе записей в road_vertices. Теперь каждый сегмент (ребро) имеет свой идентификатор (id), начальную (source) и конечную (target) вершины.

**Обновлённые записи в таблице road_segments**

Функция update_road_segments

CREATE OR REPLACE FUNCTION road_processing.update_road_segments(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$beginUPDATE road_processing.road_segments AS r-- находим точку начала сегментаSET source = (    SELECT v.id FROM road_processing.road_vertices v    WHERE part_id = p_part_id AND  layer_id = p_layer_id AND object_id = p_object_id  AND ST_DWithin(ST_StartPoint(r.geom), v.geom, 1)    LIMIT 1),-- находим точку конца сегментаtarget = (    SELECT v.id FROM road_processing.road_vertices v    WHERE part_id = p_part_id AND  layer_id = p_layer_id AND object_id = p_object_id AND ST_DWithin(ST_EndPoint(r.geom), v.geom, 1)    LIMIT 1);END;$$ LANGUAGE plpgsql;

6) Находим конечные вершины (листья). Функция fill_leaf_nodes находит на основе записей в road_segments идентификаторы (vertex_id) конечных вершин (листья) – таких, которые не связаны с другими. Результат заносит в таблицу leaf_nodes.

Функция fill_leaf_nodes

CREATE OR REPLACE FUNCTION road_processing.fill_leaf_nodes(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$begin-- удаляем старые записи из таблицы fill_leaf_nodesdelete from road_processing.leaf_nodeswhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;insert into road_processing.leaf_nodes(part_id, layer_id, object_id, vertex_id)SELECT p_part_id as part_id, p_layer_id as layer_id, p_object_id as object_id, vertex_idFROM (    SELECT source AS vertex_id FROM road_processing.road_segmentsWHERE part_id = p_part_id AND  layer_id = p_layer_id AND object_id = p_object_id     UNION ALL    SELECT target FROM road_processing.road_segmentsWHERE part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id ) AS all_verticesGROUP BY vertex_idHAVING COUNT(*) = 1;END;$$ LANGUAGE plpgsql;

7) Для каждой пары конечных вершин находим длину кратчайшего пути при помощи алгоритма Дейкстры. Функция fill_leaf_cost_matrix находит кратчайшие расстояния между всеми парами конечных вершин из leaf_nodes. Применяем алгоритм Дейкстры (pgr_dijkstraCostMatrix). Результат заносим в таблицу leaf_cost_matrix — матрица кратчайших расстояний между всеми листьями. start_vid, end_vid — идентификаторы начальной и конечной вершины, agg_cost – длинна пути между ними. Путь, идущий от start_vid до end_vid, рассчитывается как сумма всех длин составляющих его сегментов.

**Записи в таблице leaf_cost_matrix. Кратчайшие длинны путей между всеми парами конечных вершин.**

Функция fill_leaf_cost_matrix

CREATE OR REPLACE FUNCTION road_processing.fill_leaf_cost_matrix(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$begin-- удаляем старые записи из таблицы leaf_cost_matrixdelete from road_processing.leaf_cost_matrixwhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;insert into road_processing.leaf_cost_matrix(part_id, layer_id, object_id, start_vid, end_vid, agg_cost)select p_part_id as part_id, p_layer_id as layer_id, p_object_id as object_id, start_vid, end_vid, agg_costFROM pgr_dijkstraCostMatrix(    format('SELECT id, source, target, length AS cost FROM road_processing.road_segments WHERE part_id = %s AND layer_id = %s AND object_id = %L',    p_part_id, p_layer_id, p_object_id),    (    select    array_agg(vertex_id) FROM road_processing.leaf_nodes     WHERE part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id    ),    directed := false);END;$$ LANGUAGE plpgsql;

8) Находим наибольший кратчайший путь среди всех – диаметр графа. Функция fill_diameter_info находит идентификаторы начальной и конечной вершин самого длинного кратчайшего пути среди всех таких в таблице leaf_cost_matrix. Геометрия данного пути в дальнейшем будет выбрана в качестве осевой (центральной) линии полигона, идущей по середине дорожного полотна от его начала до конца.

Запись в таблице leaf_cost_matrix – идентификаторы начальной и конечной вершин самого длинного кратчайшего пути скелета дороги.

Функция fill_diameter_info

CREATE OR REPLACE FUNCTION road_processing.fill_diameter_info(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$begin-- удаляем старые записи из таблицы diameter_infodelete from road_processing.diameter_infowhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;insert into road_processing.diameter_info(part_id, layer_id, object_id, start_vid, end_vid, max_distance)SELECT    p_part_id as part_id, p_layer_id as layer_id, p_object_id as object_id, start_vid, end_vid, agg_cost AS max_distanceFROM road_processing.leaf_cost_matrixWHERE part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_idORDER BY agg_cost DESCLIMIT 1;END;$$ LANGUAGE plpgsql;

9) Получаем геометрию центральной осевой линии (centerline). Функция fill_final_longest_route находит геометрию самого длинного кратчайшего пути на основе записи таблицы leaf_cost_matrix. Результат заносится в таблицу final_longest_route.

**Запись в таблице final_longest_route – центральная (осевая) линия дороги**

10) Строим первичные зоны классификации на основе центральной линии. Функция fill_road_zones разделяет центральную линию (centerline) на одинаковые отрезки с необходимой длинной сегмента (300 метров). Параметры road_segment_len и road_segment_wide отвечают за длину и ширину зоны соответственно. Применяем функцию ST_Buffer для построения зоны (буфера) шириной road_segment_wide вокруг отрезка длиной road_segment_len. Результат заносится в таблицу road_zones.

В нашем случае длина зоны 300 метров, а ширина — 50 метров.

Функция fill_road_zones

CREATE OR REPLACE FUNCTION road_processing.fill_road_zones(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,    road_segment_len DOUBLE precision default 300, -- длинна сегмента дороги    road_segment_wide DOUBLE precision default 50 -- ширина сегмента дороги)RETURNS VOID AS $$begin-- удаляем старые записи из таблицы road_zonesdelete from road_processing.road_zoneswhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;WITH dumped_lines AS (  SELECT     (ST_Dump(geom)).geom AS line_geom,    (ST_Dump(geom)).path[1] as part_order  FROM road_processing.final_longest_route  WHERE part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id),-- Собираем все точки в правильном порядкеordered_points AS (  SELECT     ST_StartPoint(line_geom) as point,    part_order,    1 as point_type  -- start point  FROM dumped_lines  UNION ALL  SELECT     ST_EndPoint(line_geom) as point,    part_order,    2 as point_type  -- end point  FROM dumped_lines),-- Убираем дубликаты (точки соединения)unique_points AS (  SELECT DISTINCT ON (point)     point,    part_order,    point_type  FROM ordered_points  ORDER BY point, part_order, point_type),-- получеам Linestring из Multilinestringcenterline as (-- Создаем линию из точекSELECT   ST_MakeLine(point ORDER BY part_order, point_type) as linestring_geom,  ST_Length(ST_MakeLine(point ORDER BY part_order, point_type)) as total_lengthFROM unique_points),-- Разбиваем на сегменты по 300 метровfinal_segments as (SELECT   generate_series(0, floor(total_length / road_segment_len)::int) as seg_num,  ST_LineSubstring(    linestring_geom,    generate_series(0, floor(total_length / road_segment_len)::int) * road_segment_len / total_length,    LEAST((generate_series(0, floor(total_length / road_segment_len)::int) + 1) * road_segment_len / total_length, 1.0)  ) as segment_geom,  ST_Length(    ST_LineSubstring(      linestring_geom,      generate_series(0, floor(total_length / road_segment_len)::int) * road_segment_len / total_length,      LEAST((generate_series(0, floor(total_length / road_segment_len)::int) + 1) * road_segment_len / total_length, 1.0)    )  ) as segment_lengthFROM centerline)insert into road_processing.road_zones(part_id, layer_id, object_id, zone_id, geom)select p_part_id as part_id, p_layer_id as layer_id, p_object_id as object_id,row_number() over() as zone_id,ST_Buffer(segment_geom, road_segment_wide, 'endcap=flat join=mitre') as geomfrom final_segments;END;$$ LANGUAGE plpgsql;

11) Получаем более качественные зоны. Центральная линия из-за особенности полигона дороги может иметь микроскопические разрывы и зигзаги. Такое чаще всего происходит с криволинейными дорогами или дорогами, имеющими ответвления в виде примыкающих дорог, перекрёстков. Если разрывы ещё можно устранить, то справиться с зигзагами в полной мере не удаётся. Из-за этого зоны дороги не всегда идеально получаются в длине ровно по 300-метров и бывают расположены не вдоль её направления — под углом или поперёк. Выполняем функцию process_again_and_get_better_zones, строящую более качественную осевую линию и соответсвенно зоны. Эта функция объединяет все ранее полученные зоны и повторяет для неё вышерассмотренные пункты 3-11. Таким образом строится centerline по аппроксимированной геометрии дороги, полученной на основе объединения первичных зон.

**Исходный полигон и его улучшенная центральная линия.**

12) Объединяем сильно пересекающиеся по площади зоны. Функция merge_overlapping_zones в цикле объединяет (используем ST_Union) сильно пересекающиеся по площади зоны. За порог процентного значения общей площади отвечает параметр tolerance_percent.

Функция merge_overlapping_zones

CREATE OR REPLACE FUNCTION road_processing.merge_overlapping_zones(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,    tolerance_percent DOUBLE PRECISION DEFAULT 50 -- процентный порог пересечения зон. Чем меньше, тем больше по площади могут получиться зоны.)RETURNS VOID AS $$DECLARE    merged INTEGER;    v_id1 INTEGER;    v_id2 INTEGER;    v_union_geom GEOMETRY(GEOMETRY, 3857);BEGIN    LOOP        -- Находим пару для объединения        SELECT             a.zone_id,             b.zone_id,            ST_Union(a.geom, b.geom)        INTO             v_id1,             v_id2,            v_union_geom        FROM road_processing.road_zones a        JOIN road_processing.road_zones b ON a.zone_id < b.zone_id        WHERE road_processing.safe_st_intersects(a.geom, b.geom)            AND ST_Area(road_processing.safe_st_intersection(a.geom, b.geom)) >                 (tolerance_percent / 100.0) * LEAST(ST_Area(a.geom), ST_Area(b.geom))and a.part_id = p_part_id AND a.layer_id = p_layer_id AND a.object_id = p_object_idand b.part_id = p_part_id AND b.layer_id = p_layer_id AND b.object_id = p_object_id        LIMIT 1;                -- Если ничего не найдено, выходим        IF NOT FOUND THEN            EXIT;        END IF;                -- Обновляем геометрию первого полигона        UPDATE road_processing.road_zones         SET geom = v_union_geom         WHERE zone_id = v_id1 AND part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;                GET DIAGNOSTICS merged = ROW_COUNT;                -- Удаляем второй полигон        DELETE FROM road_processing.road_zones         WHERE zone_id = v_id2 AND part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;            END LOOP;END;$$ LANGUAGE plpgsql;

13) Присваиваем наложениям зон уникальные идентификаторы. Функция split_overlapped_zones выделяет все наложения зон и присваивает им уникальные идентификаторы zone_id.

**Присваиваем пересечениям зон уникальные идентификаторы.**

14) Разбиваем исходный полигон дороги на множество мелких частей (шарды). Функция fill_road_parts разбивает полигон дороги на множество малых частей (шардов) (ST_Subdivide(ST_Segmentize(geom, segment_len))). Результат заносится в таблицу road_parts. id – идентификатор шарда, zone_id – идентификатор зоны, к которой будет в последующем классифицирован каждый шард.

**Первый вариант разбиения полигона на осколки (шарды).**

Также можно использовать функцию ST_TriangulatePolygon для разбиения на треугольные шарды.

**Пример разбиения полигона на треугольные шарды.**

Функция fill_road_parts

CREATE OR REPLACE FUNCTION road_processing.fill_road_parts(    p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,    segment_len double precision default 1.0 -- длинна сегментов )RETURNS VOID AS $$BEGINdelete from road_processing.road_partswhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;WITH original_road_part as (select p_part_id as id, geometry as geomfrom  road_processing.road_solid_partswhere p_part_id = part_id  AND p_layer_id = layer_id  AND p_object_id = object_id),road_shards as(        SELECT ST_Subdivide(ST_Segmentize(geom, segment_len)) as geom         FROM original_road_part    )--    select * from road_shards    INSERT INTO road_processing.road_parts(part_id, layer_id, object_id, zone_id, id, geom)    SELECT         p_part_id as part_id,        p_layer_id as layer_id,        p_object_id as object_id,        NULL as zone_id,        row_number() OVER() as id,         geom     FROM road_shards;END;$$ LANGUAGE plpgsql;

15) Удаляем невалидные шарды. Функция fix_invalid_geometries исправляет невалидные геометрии шардов (используем ST_IsValid, ST_MakeValid) и удаляет те, которые исправить не получилось. Операции с невалидными геометриями могут вызвать ошибки в дальнейшем процессе. Примеры ошибок: самопересечение, незaмкнутость, разрыв в мультиполигоне и другие.

Функция fix_invalid_geometries

CREATE OR REPLACE FUNCTION road_processing.fix_invalid_geometries(    p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT)RETURNS VOID AS $$BEGIN--исправляем невалидные геометрии в road_partsUPDATE road_processing.road_partsSET geom = ST_MakeValid(geom)WHERE NOT ST_IsValid(geom) and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;--удаляем невалидные гемотерии из road_partsDELETE FROM road_processing.road_partsWHERE (NOT ST_IsValid(geom)) and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;END;$$ LANGUAGE plpgsql;

16) Удаляем неподходящие зоны. Функция remove_unsuitable_road_zones удаляет слишком маленькие по площади и длине зоны, в том числе невалидные, пустые и в виде точек.

**Полученные зоны в таблице road_zones.**

17) Проводим классификацию шардов на множестве зон. Функция road_part_classification производит классификацию шардов (road_parts) на множестве зон (road_zones). В результате заполняется колонка zone_id таблицы road_shard.

Классификация производится на основе следующих метрик:

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

Каждая метрика имеет свой вес, определяющий её важность в процессе классификации. Каждому шарду сопоставляется одна наиболее подходящая для него на основе метрик зона. В результате проведения классификации каждому шарду id сопоставляется зона zone_id.

Функция road_part_classification

CREATE OR REPLACE FUNCTION road_processing.road_part_classification(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,    intersection_area_weight double precision default 0.2, -- вес площади пересечения    boundary_length_weight double precision default 0.2, -- вес длины границы    centroid_distance_weight double precision default 0.2, -- вес расстояния между центроидами    min_distance_weight double precision default 0.2, -- вес минимального расстояния    lies_inside_weight double precision default 1.0 -- вес того лежит ли часть дороги целиком в зоне    )RETURNS VOID AS $$beginWITH metrics AS (    SELECT         rp.id AS shard_id,        rp.geom AS shard_geom,        rz.zone_id AS zone_id,        rz.geom AS zone_geom,        -- Площадь пересечения        ST_Area(road_processing.safe_st_intersection(rp.geom, rz.geom)) AS intersection_area,        -- Длина границы пересечения        ST_Length(road_processing.safe_st_intersection(ST_Boundary(rp.geom), ST_Boundary(rz.geom))) AS boundary_length,        -- Расстояние между центроидами        ST_Distance(ST_Centroid(rp.geom), ST_Centroid(rz.geom)) AS centroid_distance,        -- Расстояние между ближайшими точками        ST_Distance(rp.geom, rz.geom) AS min_distance,-- лежит ли часть дороги целиком в зоне        road_processing.safe_st_within(rp.geom, rz.geom) as lies_inside       FROM road_processing.road_parts rp    CROSS JOIN road_processing.road_zones rzwhere rp.part_id = p_part_id AND rp.layer_id = p_layer_id AND rp.object_id = p_object_id and rz.part_id = p_part_id AND rz.layer_id = p_layer_id AND rz.object_id = p_object_id),--select * from metrics-- Нормализуем метрики для взвешенной оценкиnormalized_metrics AS (    SELECT         shard_id,        zone_id,        intersection_area,        boundary_length,        centroid_distance,        min_distance,        lies_inside,                -- Нормализация площади пересечения (чем больше, тем лучше)        CASE             WHEN MAX(intersection_area) OVER (PARTITION BY shard_id) > 0             THEN intersection_area / MAX(intersection_area) OVER (PARTITION BY shard_id)            ELSE 0         END AS norm_intersection_area,                -- Нормализация длины границы (чем больше, тем лучше)        CASE             WHEN MAX(boundary_length) OVER (PARTITION BY shard_id) > 0             THEN boundary_length / MAX(boundary_length) OVER (PARTITION BY shard_id)            ELSE 0         END AS norm_boundary_length,                -- Нормализация расстояния между центроидами (чем меньше, тем лучше)        CASE             WHEN MAX(centroid_distance) OVER (PARTITION BY shard_id) > 0             THEN 1 - (centroid_distance / MAX(centroid_distance) OVER (PARTITION BY shard_id))            ELSE 1         END AS norm_centroid_distance,                -- Нормализация минимального расстояния (чем меньше, тем лучше)        CASE             WHEN MAX(min_distance) OVER (PARTITION BY shard_id) > 0             THEN 1 - (min_distance / MAX(min_distance) OVER (PARTITION BY shard_id))            ELSE 1         END AS norm_min_distance            FROM metrics),--select * from normalized_metrics-- Рассчитываем взвешенную оценку для каждой парыscored_matches AS (    SELECT         shard_id,        zone_id,        intersection_area,        boundary_length,        centroid_distance,        min_distance,        norm_intersection_area,        norm_boundary_length,        norm_centroid_distance,        norm_min_distance,        lies_inside,                -- Взвешенная сумма метрик (веса можно настроить в зависимости от важности)        (            COALESCE(norm_intersection_area, 0) * intersection_area_weight +      -- вес площади пересечения 40%            COALESCE(norm_boundary_length, 0) * boundary_length_weight +        -- вес длины границы 30%            COALESCE(norm_centroid_distance, 0) * centroid_distance_weight +      -- вес расстояния между центроидами 20%            COALESCE(norm_min_distance, 0) * min_distance_weight +            -- вес минимального расстояния 10%            COALESCE(lies_inside, 0) * lies_inside_weight             -- лежит ли часть дороги целиком в зоне        ) AS total_score            FROM normalized_metrics),--select * from scored_matches-- Выбираем лучшую зону для каждой части дорогиbest_matches AS (    SELECT DISTINCT ON (shard_id)        shard_id,        zone_id,        total_score,        intersection_area,        boundary_length,        centroid_distance,        min_distance,        norm_intersection_area,        norm_boundary_length,        norm_centroid_distance,        norm_min_distance    FROM scored_matches    ORDER BY shard_id, total_score DESC)--select * from best_matches-- Обновляем таблицу road_parts, присваивая zone_idUPDATE road_processing.road_parts rpSET zone_id = bm.zone_idFROM best_matches bmWHERE rp.id = bm.shard_id AND rp.part_id = p_part_id AND rp.layer_id = p_layer_id AND rp.object_id = p_object_id;END;$$ LANGUAGE plpgsql;

18) Объединяем шарды в сегменты. Функция group_zones_for_road_partition_improved производит группировку и объединение шардов по zone_id. Результат записывается в таблицу road_partition.

19) Сливаем маленькие зоны. Функция merge_small_intersected_zones циклически устраняет маленькие зоны разбиения, соединяя их с соседними. Обновляет записи в таблице road_partition.

Функция merge_small_intersected_zones

CREATE OR REPLACE FUNCTION road_processing.merge_small_intersected_zones(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,    max_area DOUBLE PRECISION DEFAULT 20000,-- максимальная площадь любого полигона    max_length DOUBLE PRECISION DEFAULT 350,-- максимальная длинна любого полигона    min_area DOUBLE PRECISION DEFAULT 5000,      -- полигоны меньше этой площади считаются "малыми"    min_length DOUBLE PRECISION DEFAULT 150      -- полигоны меньше этого размера также считаются "малыми")--RETURNS INTEGER AS $$RETURNS VOID AS $$DECLAREmerged INTEGER;    v_small_id INTEGER;    v_target_id INTEGER;    v_new_geom GEOMETRY(GEOMETRY, 3857);    v_small_geom GEOMETRY(GEOMETRY, 3857);    v_small_area DOUBLE PRECISION;    v_small_length DOUBLE PRECISION;BEGIN    LOOP        -- Находим один малый полигон (по площади или линейному размеру)        SELECT             zone_id,             geom,            ST_Area(geom),            ST_MaxDistance(geom, geom)  -- приблизительный линейный размер        INTO             v_small_id,            v_small_geom,            v_small_area,            v_small_length        FROM road_processing.road_partition        WHERE (--ST_Area(geom) <= min_area OR ST_MaxDistance(geom, geom) <= min_length)and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id  AND ST_IsValid(geom)         LIMIT 1;        --         Если малых полигонов нет - выходим        IF NOT FOUND THEN            EXIT;        END IF;        --         Ищем пересекающийся полигон с НАИМЕНЬШЕЙ бы полученной площадью и длине        SELECT             b.zone_id,            road_processing.safe_st_union(v_small_geom, b.geom)        INTO             v_target_id,            v_new_geom        FROM road_processing.road_partition b        WHERE b.zone_id != v_small_idand road_processing.safe_st_intersects(v_small_geom, b.geom)--          AND ST_Intersects(v_small_geom, b.geom)   and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id-- выбираем наименьший по площади и длине полигон        ORDER BY  (COALESCE(ST_MaxDistance(v_small_geom, v_small_geom), 0) +              COALESCE(ST_MaxDistance(b.geom, b.geom), 0)) * 0.5 +             (COALESCE(ST_Area(v_small_geom), 0) + COALESCE(ST_Area(b.geom), 0)) * 0.5 ASC        LIMIT 1;                -- Если ничего не найдено, выходим        IF NOT FOUND THEN            EXIT;        END IF;                -- Обновляем геометрию первого полигона        UPDATE road_processing.road_partition         SET geom = v_new_geom         WHERE zone_id = v_target_id and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;                GET DIAGNOSTICS merged = ROW_COUNT;                -- Удаляем малый полигон        DELETE FROM road_processing.road_partition WHERE zone_id = v_small_id and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;    END LOOP;END;$$ LANGUAGE plpgsql;

20) Сохраняем результат разбиения в таблицу final_road_partition. Функция fill_final_road_partition добавляет результат разбиения в таблицу final_road_partition.

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

Итоговое решение

Все разобранные выше этапы были объединены в одну функцию — make_part_of_road_segmentation. И достаточно передать в неё два обязательных параметра: p_layer_id – идентификатор слоя (layer_id), p_object_id – идентификатор объекта в слое (object_id). Простейшие вызов данной функции:

select road_processing.make_part_of_road_segmentation (p_layer_id => 179, p_object_id => '10000257');

Функция make_part_of_road_segmentation

CREATE OR REPLACE FUNCTION road_processing.make_part_of_road_segmentation(p_part_id INTEGER,    p_layer_id INTEGER,    p_object_id TEXT,        standart_segment_len DOUBLE PRECISION DEFAULT 300, -- желаемая длинна сегмента дороги    standart_segment_wide DOUBLE PRECISION DEFAULT 50, -- предполагаемая ширина дороги    max_segment_area DOUBLE PRECISION DEFAULT 20000,-- максимальная площадь любого полигона    max_segment_length DOUBLE PRECISION DEFAULT 350,-- максимальная длинна любого полигона    min_segment_area DOUBLE PRECISION DEFAULT 5000,     -- полигоны меньше этой площади считаются "малыми"    min_segment_length DOUBLE PRECISION DEFAULT 150,     -- полигоны меньше этого размера также считаются "малыми"    merge_area_tolerance_percent DOUBLE PRECISION DEFAULT 50, -- процентный порог зон для их слияния по общей площади их пересеченияintersection_area_weight double precision default 0.2, -- вес площади пересечения для функции классфикации шарда road_part_classification    boundary_length_weight double precision default 0.2, -- вес длины границы для функции классфикации шарда road_part_classification    centroid_distance_weight double precision default 0.2, -- вес расстояния между центроидами для функции классфикации шарда road_part_classification    min_distance_weight double precision default 0.2, -- вес минимального расстояния для функции классфикации шарда road_part_classification    lies_inside_weight double precision default 1.0, -- вес того лежит ли часть дороги целиком в зоне для функции классфикации шарда road_part_classification    topology_distance_tolerance DOUBLE PRECISION DEFAULT 2.0, --  параметр для функции generate_road_skeleton - дистанция в ST_SimplifyPreserveTopology    shard_segment_len DOUBLE PRECISION DEFAULT 1.0, -- длинна "осколка" дороги. Функция fill_road_parts, параметр для ST_Segmentize    buffer_distance DOUBLE PRECISION DEFAULT 1.0, -- величина буффер для функции ST_Buffer в process_again_and_get_better_zones    buffer_distance_for_group_zones DOUBLE PRECISION DEFAULT 0.01, -- величина буфера в функции группировки классифицированных шардов    topology_distance_tolerance_for_group_zones DOUBLE PRECISION default 0.01, -- величина допуска в функции группировки классифицированных шардов    grid_size_for_group_zones DOUBLE PRECISION default 0.1, --величина сетки в функции группировки классифицированных шардов    unsuitable_road_zones_area_threshold double precision default 100, -- избавляемся от зон с меньшей площадью    unsuitable_road_zones_len_threshold double precision default 100  -- избавляемся от зон с меньшей длинной)RETURNS VOID AS $$DECLAREv_should_skip boolean;begin--  стоит ли пропустить разбивку дороги и сразу добавить в результирующую таблицуselect road_processing.should_skip_segmentation(p_part_id, p_layer_id, p_object_id, max_segment_length) into v_should_skip;if v_should_skipthen perform road_processing.skip_segmentation_and_fill_final_road_partition(p_part_id, p_layer_id, p_object_id)return;end if;PERFORM  road_processing.generate_road_skeleton(p_part_id, p_layer_id, p_object_id, topology_distance_tolerance);--извлекаем сегменты(отдельные линии/отрезки) из "скелета" дорогиPERFORM  road_processing.fill_road_segments(p_part_id, p_layer_id, p_object_id);-- находим все соеденительные точки "скелета" дорогиPERFORM  road_processing.fill_road_vertices(p_part_id, p_layer_id, p_object_id);-- обновляем колонки source и target таблицы road_segmentsPERFORM  road_processing.update_road_segments(p_part_id, p_layer_id, p_object_id);-- находим вершины(листья, конечные узлы) "скелета" дорогиPERFORM  road_processing.fill_leaf_nodes(p_part_id, p_layer_id, p_object_id);-- вычисляем матрицу расстояний между всеми конечными узалми leaf_cost_matrixPERFORM  road_processing.fill_leaf_cost_matrix(p_part_id, p_layer_id, p_object_id);-- вычисляем диаметр дерева - самое длинное кратчайшее расстояниеPERFORM  road_processing.fill_diameter_info(p_part_id, p_layer_id, p_object_id);-- получаем геометрию диаметра - центральную осевую линию дорогиPERFORM  road_processing.fill_final_longest_route(p_part_id, p_layer_id, p_object_id);-- получаем зоны разбиения дорогиPERFORM  road_processing.fill_road_zones(p_part_id, p_layer_id, p_object_id, standart_segment_len, standart_segment_wide);-- заносим полученные зоны, заново обрабатываем их как полигон для лучшего получения зон   PERFORM road_processing.process_again_and_get_better_zones(p_part_id => p_part_id, p_layer_id => p_layer_id, p_object_id => p_object_id, standart_segment_len => standart_segment_len, standart_segment_wide => standart_segment_wide,merge_area_tolerance_percent => merge_area_tolerance_percent, buffer_distance => buffer_distance, topology_distance_tolerance => topology_distance_tolerance);-- сливаем сильно перескающиеся по площади зоныPERFORM  road_processing.merge_overlapping_zones(p_part_id, p_layer_id, p_object_id, merge_area_tolerance_percent);---- выделяем пересечения/наложения зон также в качестве обычных зонroad_processing.split_overlapped_zones(p_part_id, p_layer_id, p_object_id);---- разбиваем исходную дорогу на части    PERFORM road_processing.fill_road_parts(p_part_id, p_layer_id, p_object_id, shard_segment_len);-- исправляем/удаляем невалидные геометрииPERFORM  road_processing.fix_invalid_geometries(p_part_id, p_layer_id, p_object_id);---- удаляем неподхядящие геометрии из road_partsPERFORM road_processing.remove_unsuitable_road_parts(p_part_id, p_layer_id, p_object_id);---- удаляем неподходящие зоныPERFORM road_processing.remove_unsuitable_road_zones(p_part_id, p_layer_id, p_object_id, unsuitable_road_zones_area_threshold => unsuitable_road_zones_area_threshold,unsuitable_road_zones_len_threshold => unsuitable_road_zones_len_threshold);----проводим классфикацию кусочков дорогиPERFORM  road_processing.road_part_classification(p_part_id, p_layer_id, p_object_id,intersection_area_weight, boundary_length_weight, centroid_distance_weight,min_distance_weight, lies_inside_weight);perform road_processing.group_zones_for_road_partition (p_part_id => p_part_id, p_layer_id => p_layer_id, p_object_id => p_object_id,buffer_distance_for_group_zones => buffer_distance_for_group_zones,topology_distance_tolerance_for_group_zones => topology_distance_tolerance_for_group_zones,grid_size_for_group_zones => grid_size_for_group_zones);---- сливаем малые зоны PERFORM  road_processing.merge_small_intersected_zones(p_part_id, p_layer_id, p_object_id,max_segment_area, max_segment_length, min_segment_area, min_segment_length);---- записываем результат в таблицу final_road_partitionPERFORM  road_processing.fill_final_road_partition(p_part_id, p_layer_id, p_object_id);END;$$ LANGUAGE plpgsql;

Но функция make_part_of_road_segmentation обрабатывет лишь один part_id. Для обработки всех частей была написана функция main_road_segmentation, которая в цикле вызывает make_part_of_road_segmentation для каждого part_id.

select road_processing.main_road_segmentation(p_layer_id => 179, p_object_id => '10000257');

Функция main_road_segmentation

CREATE OR REPLACE FUNCTION road_processing.main_road_segmentation(    p_layer_id INTEGER,    p_object_id TEXT,    standart_segment_len DOUBLE PRECISION DEFAULT 300,    standart_segment_wide DOUBLE PRECISION DEFAULT 50,    max_segment_area DOUBLE PRECISION DEFAULT 20000,    max_segment_length DOUBLE PRECISION DEFAULT 350,    min_segment_area DOUBLE PRECISION DEFAULT 5000,    min_segment_length DOUBLE PRECISION DEFAULT 150,    merge_area_tolerance_percent DOUBLE PRECISION DEFAULT 50,    intersection_area_weight DOUBLE PRECISION DEFAULT 0.2,    boundary_length_weight DOUBLE PRECISION DEFAULT 0.2,    centroid_distance_weight DOUBLE PRECISION DEFAULT 0.2,    min_distance_weight DOUBLE PRECISION DEFAULT 0.2,    lies_inside_weight DOUBLE PRECISION DEFAULT 1.0,    topology_distance_tolerance DOUBLE PRECISION DEFAULT 2.0,    shard_segment_len DOUBLE PRECISION DEFAULT 1.0,    unsuitable_road_zones_area_threshold double precision default 100,     unsuitable_road_zones_len_threshold double precision default 100  )RETURNS VOID AS $$DECLARE    road_part_record RECORD;BEGIN-- Заносим части Multipolygon (мультиполигона) как отдельные part_id в таблицу source_road_solid_partsPERFORM road_processing.fill_source_road_solid_parts_multipolygon(p_layer_id, p_object_id);        -- Заносим все части дороги из таблицы source_road_solid_parts в road_solid_parts    PERFORM road_processing.fill_road_solid_parts(p_layer_id, p_object_id);        -- Удаляем старое разбиение дороги    PERFORM road_processing.delete_from_final_road_partition(p_layer_id, p_object_id);        -- Выполняем основное разбиение дороги по каждой части    FOR road_part_record IN         SELECT DISTINCT part_id         FROM road_processing.road_solid_parts        WHERE layer_id = p_layer_id AND object_id = p_object_id    LOOP        PERFORM road_processing.make_part_of_road_segmentation(            p_part_id => road_part_record.part_id,            p_layer_id => p_layer_id,            p_object_id => p_object_id,            standart_segment_len => standart_segment_len,            standart_segment_wide => standart_segment_wide,            max_segment_area => max_segment_area,            max_segment_length => max_segment_length,            min_segment_area => min_segment_area,            min_segment_length => min_segment_length,            merge_area_tolerance_percent => merge_area_tolerance_percent,            intersection_area_weight => intersection_area_weight,            boundary_length_weight => boundary_length_weight,            centroid_distance_weight => centroid_distance_weight,            min_distance_weight => min_distance_weight,            lies_inside_weight => lies_inside_weight,            topology_distance_tolerance => topology_distance_tolerance,            shard_segment_len => shard_segment_len        );    END LOOP;END;$$ LANGUAGE plpgsql;

Параметры для функций make_part_of_road_segmentation и main_road_segmentation.

№	Имя параметра	Смысловое значение	Тип	Значение по умолчанию
1	p_layer_id	Идентификатор слоя layer_id в таблице слоя features	INTEGER	—
2	p_object_id	Идентификатор объекта semantics->>’objectId’ в таблице слоя features	TEXT	—
7	standart_segment_len	Желаемая длинна сегмента дороги	DOUBLE PRECISION	300
8	standart_segment_wide	Предполагаемая ширина зоны дороги	DOUBLE PRECISION	50
9	max_segment_area	Максимальная площадь сегмента	DOUBLE PRECISION	20 000
10	max_segment_length	Максимальная длинна сегмента	DOUBLE PRECISION	350
11	min_segment_area	Сегменты меньше этой площади считаются «малыми»	DOUBLE PRECISION	5 000
12	min_segment_length	Сегменты меньше этого линейного размера считаются «малыми»	DOUBLE PRECISION	150
13	merge_area_tolerance_percent	Процентный порог для слияния зон по общей площади их пересечения	DOUBLE PRECISION	50
14	intersection_area_weight	Вес площади пересечения для функции классфикации шарда road_part_classification	DOUBLE PRECISION	0.2
15	boundary_length_weight	Вес длины границы для функции классфикации шарда road_part_classification	DOUBLE PRECISIO DOUBLE PRECISION N	0.2
16	centroid_distance_weight	Вес расстояния между центроидами для функции классфикации шарда road_part_classification	DOUBLE PRECISION	0.2
17	min_distance_weight	Вес минимального расстояния для функции классфикации шарда road_part_classification	DOUBLE PRECISION	0.2
18	lies_inside_weight	Вес того лежит ли часть дороги целиком в зоне для функции классфикации шарда road_part_classification	DOUBLE PRECISION	1.0
19	topology_distance_tolerance	Параметр для функции generate_road_skeleton — дистанция в ST_SimplifyPreserveTopology	DOUBLE PRECISION	2.0
20	shard_segment_len	Длинна «осколка» дороги. Функция fill_road_parts, параметр для ST_Segmentize	DOUBLE PRECISION	1.0
25	unsuitable_road_zones_area_threshold	Избавляемся от зон с меньшей площадью в функции remove_unsuitable_road_zones	DOUBLE PRECISION	100
26	unsuitable_road_zones_len_threshold	Избавляемся от зон с меньшей длиной в функции remove_unsuitable_road_zones	DOUBLE PRECISION	100

Пример создания таблицы:

create table if not exists road_processing.road_parts(part_id INTEGER,    layer_id INTEGER,    object_id TEXT,    zone_id INTEGER, id INTEGER,    geom GEOMETRY(GEOMETRY, 3857) );

В схеме road_processing были созданы таблицы, индексы и функции, необходимые для работы функций. Тем самым мы вынесли в отдельную область наше решение, сделали своеобразное расширение (extension).

Для ускорения поиска практически для каждой таблицы создаётся составной индекс по полям object_id, layer_id, part_id.

CREATE INDEX idx_road_partsON road_processing.road_parts (object_id, layer_id, part_id);

А также индекс типа GIST по геометрическому полю.

CREATE INDEX idx_road_parts_geom ON road_processing.road_parts USING GIST (geom);

Примеры сегментации различных по геометрии дорог

№1

-- выполняем основную функциюselect road_processing.main_road_segmentation(p_layer_id => 179, p_object_id => '10000005',merge_area_tolerance_percent => 50, min_segment_length => 250, max_segment_length => 350, shard_segment_len => 1);-- смотрим на финальный результатselect *, ST_Area(geom), ST_MaxDistance(geom, geom) as len from road_processing.final_road_partitionwhere layer_id = 179 and object_id = '10000005';

Визуально алгоритм хорошо справился с поставленной задачей на данном полигоне, разрезав его на 12 сегментов.

Но есть некоторые недостатки:

В колонке len отображается максимальная длина сегмента. Для упрощения используем функцию ST_MaxDistance, которая возвращает расстояние между двумя максимально удалёнными точками полигона. Поэтому данные значения являются приближёнными. Тем не менее, существует разброс длины сегментов: min(len) – 276, max(len) – 652. Хотелось бы минимизировать данный разлёт, чтобы все длины сегментов были бы в диапазоне 300-350 метров.
Малые и большие сегменты встречаются в середине полигона.
Линии разрезов неидеально перпендикулярны направлению дороги.

№2

1) Первый вариант параметров:

-- выполняем основную функциюselect road_processing.main_road_segmentation(p_layer_id => 179, p_object_id => '10000011',merge_area_tolerance_percent => 50, min_segment_length => 250, max_segment_length => 350,shard_segment_len => 1);

В итоге полигон разрезан на две части.

**Сегментация при первом варианте параметров**

2) Второй вариант параметров.

Изменили min_segment_length: вместо 250 поставили 150. Получили три части вместо двух.

**Разрезка дороги с уменьшенным значением допустимой длины сегмента**

№3

**Сегментация разделяющейся на две полосы дороги**

№4

**Замкнутая территория разбилась на две части**

№5

select road_processing.main_road_segmentation(p_layer_id => 179, p_object_id => '10000991',merge_area_tolerance_percent => 50, min_segment_length => 250, max_segment_length => 350,shard_segment_len => 1);

Полигон дороги считается целиком отдельным сегментом так как имеет подходящие размеры.

**Исходный полигон как отдельный сегмент**

№6

№7

**Сегментация длинной дороги с прямоугольным поворотом**

№8

Алгоритм не предназначен для сегментации сети дорог. Это связано с тем, что в его основе лежит нахождение единственной осевой (центральной) линии полигона. И разрезка происходит вдоль неё.

**Полигон дороги с большим ответвлением**

**Центральная (осевая) линия полигона не имеет развязок**

**Центральная (осевая) линия и исходный полигон**

№9

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

select road_processing.main_road_segmentation(p_layer_id => 179,p_object_id => '10000188', process_for_better_zones => true, process_for_fill_road_parts_considering_zones => true,merge_area_tolerance_percent => 50, shard_segment_len => 1);-- Смотрим на результат:select *, ST_Area(geom), ST_MaxDistance(geom, geom) as lenfrom road_processing.final_road_partitionwhere layer_id = 179 and object_id = '10000188';

Обратите внимание, что исходная дорога имеет два разрыва и каждая часть имеет свой part_id.

Заключение

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

Он не предназначен для обработки сети дорог — дорог с большими широкими перекрёстками, с крупными ответвлениями
Для некоторых полигонов длины получающихся сегментов превышают установленное стандартное значение
В середине полигона встречаются большие и маленькие сегменты
Линии разреза сегментов неидеально перпендикулярны направлению дороги

Подбор параметров помогает улучшить результаты.

Время выполнения алгоритма зависит от длинны полигона и обычно занимает 5-30 секунд. Большие полигоны обрабатываются дольше малых.

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

А как бы Вы решали данную задачу? Пишите свои варианты в комментариях. Задавайте вопросы, мы с радостью на них ответим.

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

Спасибо за внимание!

Ссылки, источники, материалы

Stephen Vincent Mather, Pedro Wightma, Bborie Park, Thomas Kraft. PostGIS Cookbook: Store, organize, manipulate, and analyze spatial data, Second Edition, 2018, стр. 576
Документация Postgres
Документация по функциям PostGIS

ссылка на оригинал статьи https://habr.com/ru/articles/1040940/