Здравствуйте, уважаемые читатели Хабра!
В первой части мы рассказали, как решили задачу сегментации полигона дороги в PostGIS. А теперь рассмотрим сопоставление сегментов двух разных версий дороги для сохранения учёта историчности привязанным к ним событий.
Содержание
Введение
Мы в команде занимаемся построением высоконагруженных информационных систем. Помогаем клиентам быстро находить верные ответы благодаря системам поддержки принятия решений, бизнес-аналитике и построению отчётности. Работаем над проектами в сферах Smart City и Smart Transport.
В рамках проекта по созданию информационной системы мониторинга событий в городе перед нами стала задача создания функций для управления цифровыми активами автомобильных дорог. В предыдущей статье мы рассмотрели реализацию алгоритма сегментации полигона дороги в PostGIS. Теперь же обратим внимание на не менее важную задачу – сопоставление сегментов полигона двух разных версий дороги.
Что мы хотим получить и зачем нам это нужно? На рисунке ниже представлены две сегментированные версии дороги. Синим цветом обозначена актуальная часть дороги, красным – старая. Актуальная дорога находится на картинке под предыдущей. Новая дорога получена путём продления старой. Обе версии полигона были сегментированы и их общие части полностью совпадают. Это тот идеальный результат, который в итоге мы хотим получить.
Напомню, что каждое связанное с дорогой событие принадлежит сегменту дороги, а не всему многокилометровому полигону. Это позволяет отслеживать с течением времени обстановку на участках, выявлять среди них проблемные и опасные.
Чтобы было совсем понятно, приведём следующий пример. Изначально была 1-я версия дороги, её сегментировали и к её 5-му участку привязали инцидент «отсутствие разметки». Затем спустя некоторое время всю дорогу решили удлинить и расширить. Это означает, что у нас появляется новая 2-ая версия цельного полигона дороги. Затем если мы её просто разрежем на сегменты, то с высокой долей вероятности наш 5-й участок будет иметь совсем другой идентификатор, и мы потеряем привязанный к нему инцидент. Но если есть механизм маппинга сегментов, то 5-й участок и привязанное к нему событие останутся на своих местах. Однако стоит отметить, что если произошли существенные изменения геометрии участка, при которых его новая и старая версия несопоставимы, то нам следует участку в новой версии дороги присвоить ранее неиспользовавшийся уникальный идентификатор.
Таким образом, функция маппинга сегментов должна сопоставлять участки двух разных дорог, а при их существенном различии присваивать новому участку ранее неиспользовавшийся уникальный идентификатор. В итоге при вставке новой версии полигона в таблицу (срабатывает по триггеру) должна появляться новая запись (INSERT) в целевой таблице активов со следующими заполненными полями:
-
layer_id — идентификатор слоя,
-
object_id — идентификатор объекта дороги в слое,
-
zone_id — идентификатор сегмента,
-
geom — геометрия сегмента.
Алгоритм сопоставления сегментов дороги
Для демонстрации качества работы каждого этапа алгоритма сначала применим его для не отличающихся по геометрии полигонов двух версий дороги. То есть представим, что к нам пришло крупное обновление дорожной сети, в которой данная дорога осталась прежней.
1) Заносим старое разбиение дороги. Функция fill_old_road_partition заносит старое разбиение дороги в таблицу old_road_partition.
2) Заносим новую цельную дорогу. Функция fill_new_solid_road добавляет геометрию новой цельной дороги в таблицу new_solid_road.
3) Разбиваем новый полигон дороги на шарды. Функция fill_new_road_shards разбивает новый полигон дороги на шарды и заполняет таблицу new_road_shards (применяем ST_Subdivide(ST_Segmentize(geom, segment_len))) аналогично как это делалось на 14-ом этапе в функции fill_road_parts в задаче сегментации.
4) Классифицируем шарды. Функция new_road_shards_classification проводит классификацию каждого шарда на множестве сегментов старой дороги аналогично как это производилось на 17-ом этапе в функции road_part_classification при разрезке дороги с небольшим отличием в том, что здесь в качестве зон классификации выступают сегменты старого полигона. По итогу обновляется колонка zone_id в таблице new_road_shards.
5) Объединяем шарды в сегменты. Функция group_zones_for_new_road_partition группирует шарды по одинаковым значениям zone_id, сливает их геометрии и тем самым формирует сегменты. Результат заносится в таблицу new_road_partition.
Функция group_zones_for_new_road_partition
-- группируем зоны и заполнеям таблицу new_road_partitionCREATE OR REPLACE FUNCTION road_processing.group_zones_for_new_road_partition(p_part_id INTEGER, p_layer_id INTEGER, p_object_id TEXT)RETURNS VOID AS $$begin-- удаляем старый результатdelete from road_processing.new_road_partitionwhere part_id = p_part_id and layer_id = p_layer_id and object_id = p_object_id;insert into road_processing.new_road_partition(part_id, layer_id, object_id, zone_id, geom)SELECTp_part_id as part_id, p_layer_id as layer_id, p_object_id as object_id, zone_id,ST_Union(geom) as geomFROM road_processing.new_road_shardswhere part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_idgroup by zone_id;END;$$ LANGUAGE plpgsql;
В целом сопоставление сегментов прошло успешно. Хотя и можно выявить небольшие различия в границах участках и неперпендикулярность разрезов.
При маппинге могут быть следующие варианты:
-
Геометрия новой дороги не изменились или изменились незначительно. Этот вариант мы только что разобрали.
-
Геометрия новой дороги уменьшилась в размерах.
-
Геометрия новой дороги увеличилась в размерах.
4) Обрабатываем большие сегменты в случае увеличения дороги.
Функция fill_diff_zones определяет различия между старыми и получившимися новыми зонами с одинаковыми идентификаторами (zone_id).
Учитываем во внимание относительную разницу между площадями и периметрами. Результат заносится в таблицу diff_zones.
Функция fill_diff_zones
-- функция для определения различий старых и новых зон при маппингеCREATE OR REPLACE FUNCTION road_processing.fill_diff_zones(p_part_id INTEGER, p_layer_id INTEGER, p_object_id TEXT)RETURNS VOID AS $$begin-- удаляем старые записиdelete from road_processing.diff_zoneswhere part_id = p_part_id and p_layer_id = layer_id and p_object_id = object_id;insert into road_processing.diff_zones(part_id, layer_id, object_id, zone_id, diff_area, diff_len)select p_part_id as part_id,p_layer_id as layer_id,p_object_id as object_id,old_road.zone_id as zone_id,(ST_Area(old_road.geom) - ST_Area(new_road.geom)) / NULLIF(ST_Area(old_road.geom), 0) * 100 as diff_area,(ST_Length(ST_Boundary(old_road.geom))-ST_Length(ST_Boundary(new_road.geom))) / NULLIF(ST_Length(ST_Boundary(new_road.geom)), 0) * 100 as diff_lenfrom road_processing.old_road_partition as old_roadjoin road_processing.new_road_partition as new_road on old_road.zone_id = new_road.zone_id where old_road.part_id = p_part_id AND old_road.layer_id = p_layer_id AND old_road.object_id = p_object_idand new_road.part_id = p_part_id AND new_road.layer_id = p_layer_id AND new_road.object_id = p_object_id;END;$$ LANGUAGE plpgsql;
Затем функция fill_large_zones выбирает увеличенные зоны из таблицы diff_zones на основе относительных пороговых значений отличий площадей и периметров, добавляет их в таблицу large_zones. То есть на данном этапе отбираем слишком большие зоны.Так на рисунке ниже крайний сверху сегмент старой дороги (красный) выступает в качестве наиболее подходящей зоны для классификации достроенного участка (синий).
Функция fill_large_zones
CREATE OR REPLACE FUNCTION road_processing.fill_large_zones(p_part_id INTEGER, p_layer_id INTEGER, p_object_id TEXT, area_threshold double precision default 50, len_threshold double precision default 50, standart_segment_len double precision default 300, diff_len_coefficient double precision default 2)RETURNS VOID AS $$begin-- удаляем старые записиdelete from road_processing.large_zoneswhere part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id;-- находим увеличенные зоныwith increased_zones as (select zone_idfrom road_processing.diff_zoneswhere (diff_area < 0 and abs(diff_area) >= area_threshold ) or (diff_len < 0 and abs(diff_len) >= len_threshold ) and part_id = p_part_id AND layer_id = p_layer_id AND object_id = p_object_id),-- отбираем большие зоны для разбиенияzones_for_partitions as (select iz.zone_id as zone_id, nrp.geom as geomfrom increased_zones as iz join road_processing.new_road_partition as nrp on iz.zone_id = nrp.zone_idwhere ST_Length(ST_Boundary(nrp.geom))>= diff_len_coefficient*standart_segment_len AND nrp.part_id = p_part_id AND nrp.layer_id = p_layer_id AND nrp.object_id = p_object_id)insert into road_processing.large_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,zone_id as zone_id,geomfrom zones_for_partitions;END;$$ LANGUAGE plpgsql;
Функция make_large_zones_partition производит разбиение больших зон на сегменты из таблицы large_zones. Для этого воспользуемся функцией make_part_of_road_segmentation из задачи сегментации. Результат заносится в таблицу new_large_zones_partition.
После чего функция new_large_zones_classification заново производит классификацию для зон из new_large_zones_partition.
И уже функция manage_new_large_zones добавляет зоны из таблицы new_large_zones_partition в new_road_partition.
5) Обрабатываем малые сегменты. Функция merge_small_intersected_mapped_zones итеративно находит и сливает смежные малые зоны. Обновляет записи в финальной таблице new_road_partition.
Функция merge_small_intersected_mapped_zones
CREATE OR REPLACE FUNCTION road_processing.merge_small_intersected_mapped_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.new_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 LIMIT 1; -- Если малых полигонов нет - выходим IF NOT FOUND THEN EXIT; END IF; -- Ищем пересекающийся полигон с НАИМЕНЬШЕЙ бы полученной площадью и длинне SELECT b.zone_id, ST_Union(v_small_geom, b.geom) INTO v_target_id, v_new_geom FROM road_processing.new_road_partition b WHERE b.zone_id != v_small_id 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.new_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.new_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;
Итоговое решение
Все вышеперечисленные этапы были включены в функцию make_part_of_road_mapping, она обрабатывает один полигон дороги. Как и в задаче сегментации для обработки мультиполигона она вызывается в цикле функции main_road_mapping для каждого находящегося в нём полигона (part_id).
Функция make_part_of_road_mapping
CREATE OR REPLACE FUNCTION road_processing.make_part_of_road_mapping(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, -- величина буфера в функции группировки классифицированных шардов unsuitable_road_zones_area_threshold double precision default 100, -- избавляемся от зон с меньшей площадью unsuitable_road_zones_len_threshold double precision default 100, -- избавляемся от зон с меньшей длинной area_threshold double precision default 50, -- процентный порог разницы площадей зон для функции fill_large_zones len_threshold DOUBLE PRECISION DEFAULT 50, -- процентный порог разницы длин зон для функции fill_large_zones diff_len_coefficient DOUBLE PRECISION DEFAULT 2.0 -- коэффициент разницы длины зоны функции fill_large_zones)RETURNS VOID AS $$begin-- разбиваем дорогу на шардыperform road_processing.fill_new_road_shards(p_part_id => p_part_id, p_object_id => p_object_id, p_layer_id => p_layer_id, segment_len => shard_segment_len);-- проводим классификацию каждого шардаperform road_processing.new_road_shards_classification(p_part_id => p_part_id, p_object_id => p_object_id, p_layer_id => p_layer_id,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);-- группируем зоныperform road_processing.group_zones_for_new_road_partition(p_part_id => p_part_id, p_object_id => p_object_id, p_layer_id => p_layer_id);-- определяем различия зонperform road_processing.fill_diff_zones(p_part_id => p_part_id, p_layer_id => p_layer_id, p_object_id => p_object_id);-- определяем наибольшие зоныperform road_processing.fill_large_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,area_threshold => area_threshold, len_threshold => len_threshold,diff_len_coefficient=> diff_len_coefficient);-- разбиваем большие зоныperform road_processing.make_large_zones_partition(p_part_id => p_part_id, p_layer_id => p_layer_id, p_object_id => p_object_id,process_for_better_zones => process_for_better_zones, process_for_split_overlapped_zones => process_for_split_overlapped_zones,process_for_fill_road_parts_considering_zones => process_for_fill_road_parts_considering_zones, standart_segment_len => standart_segment_len, standart_segment_wide => standart_segment_wide, max_segment_area => max_segment_area, min_segment_area => min_segment_area, max_segment_length => max_segment_length, 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,shard_segment_len => shard_segment_len,buffer_distance => buffer_distance,unsuitable_road_zones_area_threshold => unsuitable_road_zones_area_threshold,unsuitable_road_zones_len_threshold => unsuitable_road_zones_len_threshold);-- производим маппинг для разбитых сегментов из больших зон (сопоставляем zone_id)perform road_processing.new_large_zones_classification(p_layer_id => p_layer_id, p_object_id => p_object_id,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);-- добавляем классифицированные сегменты из больших зон в итоговое разбиениеperform road_processing.manage_new_large_zones(p_layer_id => p_layer_id,p_object_id => p_object_id);-- объединяем малые зоныperform road_processing.merge_small_intersected_mapped_zones(p_part_id => p_part_id,p_layer_id => p_layer_id,p_object_id => p_object_id,max_area => max_segment_area, max_length => max_segment_length,min_area => min_segment_area, min_length => min_segment_length);END;$$ LANGUAGE plpgsql;
Параметры функции main_road_mapping
|
№ |
Имя параметра |
Смысловое значение |
Тип |
Значение по умолчанию |
|
1 |
p_layer_id |
Идентификатор слоя |
INTEGER |
— |
|
2 |
p_object_id |
Идентификатор объекта |
TEXT |
— |
|
3 |
standart_segment_len |
Желаемая длинна сегмента дороги |
DOUBLE PRECISION |
300 |
|
4 |
standart_segment_wide |
Предполагаемая ширина зоны дороги |
DOUBLE PRECISION |
50 |
|
5 |
max_segment_area |
Максимальная площадь сегмента |
DOUBLE PRECISION |
20 000 |
|
6 |
max_segment_length |
Максимальная длинна сегмента |
DOUBLE PRECISION |
350 |
|
7 |
min_segment_area |
Сегменты меньше этой площади считаются «малыми» |
DOUBLE PRECISION |
5 000 |
|
8 |
min_segment_length |
Сегменты меньше этого линейного размера считаются «малыми» |
DOUBLE PRECISION |
150 |
|
9 |
merge_area_tolerance_percent |
Процентный порог для слияния зон по общей площади их пересечения |
DOUBLE PRECISION |
50 |
|
10 |
intersection_area_weight |
Вес площади пересечения для функции классификации шарда road_part_classification |
DOUBLE PRECISION |
0.2 |
|
11 |
boundary_length_weight |
Вес длины границы для функции классификации шарда road_part_classification |
DOUBLE PRECISION |
0.2 |
|
12 |
centroid_distance_weight |
Вес расстояния между центроидами для функции классификации шарда road_part_classification |
DOUBLE PRECISION |
0.2 |
|
13 |
min_distance_weight |
Вес минимального расстояния для функции классификации шарда road_part_classification |
DOUBLE PRECISION |
0.2 |
|
14 |
lies_inside_weight |
Вес того лежит ли часть дороги целиком в зоне для функции классификации шарда road_part_classification |
DOUBLE PRECISION |
1.0 |
|
15 |
topology_distance_tolerance |
Параметр для функции generate_road_skeleton – дистанция в ST_SimplifyPreserveTopology |
DOUBLE PRECISION |
2.0 |
|
16 |
shard_segment_len |
Длинна «осколка» дороги. Функция fill_road_parts, параметр для ST_Segmentize |
DOUBLE PRECISION |
1.0 |
|
17 |
area_threshold |
Процентный порог разницы площадей зон для функции fill_large_zones |
DOUBLE PRECISION |
50 |
|
18 |
len_threshold |
Процентный порог разницы длин зон для функции fill_large_zones |
DOUBLE PRECISION |
50 |
|
19 |
diff_len_coefficient |
Коэффициент разницы длины зон для функции fill_large_zones |
DOUBLE PRECISION |
2.0 |
Основная функция — main_road_mapping.
select road_processing.main_road_mapping(p_layer_id => 179, p_object_id => '10000257');
Примеры маппинга
№1
№2
Маппинг мультиполигона – каждый полигон имеет свой идентификатор part_id.
№3
№4
№5
№6
№7
№8
№9
№10
select road_processing.main_road_mapping(p_layer_id => 179, p_object_id => '10000502', shard_segment_len => 1);
№11
№12
select road_processing.main_road_mapping(p_layer_id => 179, p_object_id => '10001480', shard_segment_len => 1);
№13
Заключение
Разработанный алгоритм маппинга сохраняет историчность событий при изменениях связанных с ними территорий.
Время его выполнения зависит от длины полигона и в среднем составляет от 5 до 30 секунд: чем крупнее полигон, тем дольше обработка.
Основные недостатки алгоритма — небольшие расхождения в границах участков и отсутствие строгой перпендикулярности разрезов.
Алгоритмы сегментации и маппинга интегрированы в систему управления цифровыми активами дорог. В некоторых случаях всё ещё требуется участие специалиста, однако, по нашим оценкам, около 90% дорог могут обрабатываться полностью автоматически. Это обеспечивает актуальность данных в системе и значительно сокращает время сотрудников на рутинные операции.
Как бы вы решили эту задачу? Делитесь идеями в комментариях. Задавайте вопросы.
Спасибо за внимание!
Ссылки, источники, материалы
-
Stephen Vincent Mather, Pedro Wightma, Bborie Park, Thomas Kraft. PostGIS Cookbook: Store, organize, manipulate, and analyze spatial data, Second Edition, 2018, стр. 576
ссылка на оригинал статьи https://habr.com/ru/articles/1042618/