Здравствуйте, уважаемые читатели Хабра!
В серии статей хочу рассказать о создании основного функционала MVP (Minimum Value Product) системы по управлению цифровыми активами для базы данных PostGIS. Полный перечень возможностей разрабатываемого проекта представлен на картинке ниже.

В этой публикации будут рассмотрены связанные с геометрией две задачи:
Кластеризация всех объектов в базе данных на основе геометрического поля
Поиск одинаковых и похожих объектов на основе исходной геометрии
Содержание:
Введение
Для начала несколько слов о том, чем занимаюсь и откуда вообще возникла необходимость построения такой системы. В данный момент работаю на позиции старшего разработчика Отдела информационных систем Управления автоматизации. В неполный перечень моих задач входит разработка плагинов по созданию автоматизированных отчётов для QGIS на Python; создание микросервисов, оперирующих геометрическими слоями.

Одним из основных источником данных в моей работе и для сотрудников является база данных PostGIS. PostGIS представляет собой расширение для PostgreSQL, поддерживающее работу с географическими объектами. У нас в ней насчитывается более чем 1200 несвязанных между собой таблиц, в которых суммарно содержится сотни миллионов объектов (304’559’489 на данный момент, если быть точнее) и этот объём продолжает расти. Как вы поняли у нас очень много данных.

Какие задачи приходится решать коллегам и как DAM-система сможет им в этом помочь:
Нахождение одинаковых объектов. Представьте у вас есть запись в таблице "schema_name_1"."table_name_1" (полное название таблицы это "имя_схемы"."имя_таблицы") о каком-то объекте. Этим объектом может быть строение, жилое здание, школа, поликлиника, детский сад и т.д. В этой записи есть геометрия объекта и подходящие для вашей задачи некоторые колонки, но не все. Вы помните, что где-то в этом огромном океане данных есть еще одна или несколько таблиц, в которых находится информация об этом объекте с другими так нужными вам колонками. Вот только не знаете точно где...И тут либо сидеть вспоминать, либо перебирать имена схем и таблиц; либо тревожить товарищей и аналитиков в надежде, что они уж точно помнят и быстро ответят. А как здорово было бы иметь функционал, за мгновения выдающий имена таблиц и идентификаторы одинаковых по геометрии объектов.
Нахождение похожих объектов. Аналогично предыдущему пункту, только теперь объекты неравны. Геометрия объектов со временем может меняться: например, могут добавляться или сноситься пр��стройки. К тому же в двух таблицах геометрии одного объекта могут немного отличаться друг от друга. Это представляет собой отдельный интерес и может давать инсайты.
Кластеризация. Объединение объектов в группы на основе их местоположения позволяет оценить как распределены объекты и помочь в сборе информации о территории. Этот функционал по большей части представляет интерес для аналитиков. Например, это может служить для обнаружения аномально плотных или разреженных зон.
Исходные данные
Давайте посмотрим на то, как выглядят наши данные. На картинке ниже представлены несколько записей. Видим несколько колонок, крайнее представляет собой геометрическое поле с именем wkb_geometry в формате WKB (Well-Known Binary).

В pgAdmin можно посмотреть на геометрии объектов.


Теперь, когда задачи ясны и понятны имеющиеся данные, можно приступать к само��у интересному.
Кластеризация объектов по геометрии
Предварительные сведения
Всего существуют 4 оконные функции PostGIS, возвращающие номер кластера для каждой входной геометрии:
integer ST_ClusterDBSCAN(geometry winset geom, float8 eps, integer minpoints);
Применяет алгоритм DBSCAN. Не требует указания количества кластеров. На вход передаём геометрию, а так же параметры расстояния (eps) и плотности (minpoints). Кластеры могут иметь произвольную форму.
Входная геометрия добавляется в кластер, если она: «Основная» - находящаяся в пределах расстояния eps не менее чем от minpoints входных геометрий (включая себя); или «Граничная», находящаяся в пределах расстояния eps от основной геометрии. А также есть «Шумовые» объекты - те, которые не попали в первые две категории. Их значение в cluster_id - null. Могут представлять дополнительный интерес.
Более подробно с описанием данной функции можно ознакомиться здесь.
integer ST_ClusterKMeans(geometry winset geom, integer k, float8 max_radius);
Применяет алгоритм KMeans. За расстояние между двумя объектами принимается расстояние между их центроидами. Параметр k отвечает за кол-во кластеров. Если задан параметр max_radius, то сгенерирует больше кластеров, чем k, гарантируя, что ни один кластер в выходных данных не будет иметь радиус больше чем max_radius. Кластеры имеют сферическую форму.
Детальнее с ней можно ознакомиться здесь.
integer ST_ClusterIntersectingWin(geometry winset geom);
Строит связанные кластеры пересекающихся геометрических объектов. Работает по принципу транзитивности: если A пересекается с B, а B с C, то все три попадут в один кластер.
Более подробно она описана тут.
integer ST_ClusterWithinWin(geometry winset geom, float8 distance);
Объекты попадают в один кластер, если расстояние между ними не более чем distance. Вызов ST_ClusterWithinWin эквивалентен запуску ST_ClusterDBSCAN с minpoints = 0.
Более подробно прочитать о ней можно тут.
Кластеризация двух таблиц
Перед нами стоит задача кластеризации всех объектов в БД. Но давайте начнём с малого. Ниже представлен SQL-запрос для кластеризации объектов двух таблиц: 'osnova_marketing_2024'.'building' и 'egko'.'egko_building'.
Кластеризация двух таблиц
-- получаем объекты из двух таблиц WITH all_geometries AS ( SELECT 'osnova_marketing_2024' as schema_name, 'building' as table_name, id as id, wkb_geometry as wkb_geometry FROM osnova_marketing_2024.building WHERE wkb_geometry IS NOT NULL UNION ALL SELECT 'egko' as schema_name, 'egko_building' as table_name, id as id, wkb_geometry as wkb_geometry FROM egko.egko_building WHERE wkb_geometry IS NOT NULL ), -- проводим кластеризацию clustered AS ( SELECT ST_ClusterDBSCAN(wkb_geometry, eps := 100, minpoints := 5) OVER() AS cluster_id, -- ST_ClusterKMeans(wkb_geometry, k := 5, max_radius := 100) OVER() AS cluster_id, -- ST_ClusterIntersectingWin(wkb_geometry) OVER() AS cluster_id, -- ST_ClusterWithinWin(wkb_geometry, distance:= 50) OVER() AS cluster_id, schema_name, table_name, id, ST_Transform(wkb_geometry, 4326) as wkb_geometry FROM all_geometries ), -- производим фильтрацию filtered_clustered AS ( SELECT cluster_id, schema_name, table_name, id, wkb_geometry FROM clustered WHERE cluster_id IS NOT NULL ), -- формируем геометрии кластеров clusters AS ( SELECT cluster_id, NULL as schema_name, NULL as table_name, NULL as id, ST_ConvexHull(ST_Collect(wkb_geometry)) as wkb_geometry FROM filtered_clustered GROUP BY cluster_id ) -- вывод результатов SELECT cluster_id, schema_name, table_name, id::text, wkb_geometry FROM clusters UNION ALL SELECT cluster_id, schema_name, table_name, id::text, wkb_geometry FROM filtered_clustered
CTE (Common Table Expressions - обобщённые табличные выражения) all_geometries служит для объединения всех объектов двух таблиц. Забираем имена схем и таблиц, идентификаторов объектов и геометрию. Причём не имеющие геометрии объекты нас не интересуют - указываем это в операторе WHERE.
CTE clustered осуществляет кластеризацию всех объектов, полученных в результате предыдущего подзапроса. Здесь выбираем понравившуюся нам функцию, передаём ей аргументы. И с помощью функции ST_Transform переводим объекты в систему координат WGS 84 (World Geodetic System 1984).
В CTE filtered_clustered отбрасывает объекты, которые не попали в кластер.
В CTE clusters получаем геометрии самих кластеров. Группируем объекты по cluster_id, объединяем их геометрии в одну при помощи функции ST_Collect, и для формирования контуров кластеров используем функцию ST_ConvexHull.
Примечание: ST_ConvexHull даёт более чёткое отображение кластеров по сравнению с ST_Envelope, представляющая кластера в виде прямоугольников, но выполняется дольше.
И наконец объединяем геометрии объектов и кластеров. В результате получаем следующую ниже таблицу. В ней для объектов указан номер кластера, из какой таблицы он был взят, его идентификатор и геометрия. У кластеров поля schema_name и table_name имеют значения NULL.



Ура! Кластеризация двух таблиц прошла успешно! Самое время замахнуться на всё.
Кластеризация всех объектов БД
Суммарно в этих двух таблицах 1 500 000 объектов. Они кластеризуются приблизительно за 60 сек.(1 минута). Всего же в базе дан��ых 304 559 489 записей. Если предположить, что время кластеризации линейно зависит от кол-ва объектов, то кластеризация всей БД по геометрии займёт ~ 203 минуты(3-4 часа). Напишем на Python функцию, создающую SQL-запрос, включающий все таблицы. Получим запрос из 77 571 строк, 77 553 из них занимает CTE all_geometries.
Попытка кластеризовать всё в лоб
WITH all_geometries AS ( SELECT '_1059' as schema_name, 'dorogi_federalnogo_znachenija' as table_name, ogc_fid as id, wkb_geometry as wkb_geometry FROM _1059."dorogi_federalnogo_znachenija" WHERE wkb_geometry IS NOT NULL UNION ALL SELECT '_1059' as schema_name, 'dorogi_obschegorodskogo_znachenija' as table_name, ogc_fid as id, wkb_geometry as wkb_geometry FROM _1059."dorogi_obschegorodskogo_znachenija" WHERE wkb_geometry IS NOT NULL . . . UNION ALL SELECT '_1059' as schema_name, 'dorogi_rajonnogo_znachenija' as table_name, ogc_fid as id, wkb_geometry as wkb_geometry FROM _1059."dorogi_rajonnogo_znachenija"), clustered AS (...
Такой километровый запрос PostGIS обработать не в состоянии. Но выход есть. Нужно сначала создать таблицу со всеми объектами и затем использовать её вместо CTE all_geometries. С кластеризацией разобрались. Переходим к поиску объектов.
Поиск одинаковых и похожих объектов
One vs All
Начнём с малого. Пусть нам надо по геометрии из одной таблицы найти такие же или похожие на неё объекты в другой таблице.
В целях уменьшения времени поиска создадим индекс GIST для геометрических полей таблиц.
CREATE INDEX IF NOT EXISTS idx_geom_public_osnova_marketing_2024__building ON public.osnova_marketing_2024__building USING GIST (wkb_geometry); CREATE INDEX IF NOT EXISTS idx_geom_public_egko__egko_building ON public.egko__egko_building USING GIST (wkb_geometry);
Для проверки на равенство геометрий у PostGIS есть булева функция ST_equals(geometry A, geometry B). Но что делать, если нам надо найти не полностью одинаковые, но похожие объекты? Рассмотрим варианты:
Дискретное расстояние Хаусдорфа. Если попытаться объяснить простыми словами, то это наибольшее из всех расстояний от точки в одном множестве до ближайшей точки в другом множестве. В PostGIS вычисляется при помощи вызова функции float ST_HausdorffDistance(geometry g1, geometry g2). Более детально можно ознакомиться с ней здесь.
Дискретное расстояние Фреше ST_FrechetDistance. В Википедии можно найти такое неформальное определение
Представьте себе человека, идущего по конечному криволинейному пути, выгуливающего свою собаку на поводке, причем собака пересекает отдельный конечный изогнутый путь. Каждый из них может изменять свою скорость, чтобы не ослаблять поводок, но ни один из них не может двигаться назад. Расстояние Фреше между двумя кривыми — это длина кратчайшего поводка, достаточная для того, чтобы оба могли пройти свои отдельные пути от начала до конца.
Более детально можно ознакомиться с ней здесь.
Своя метрика на основе расстояния между центроидами, относительных значений площадей и периметров. Понадобятся следующие функции: ST_Centroid - для получения центроида, ST_Area и ST_Perimeter для получения площади и периметра объекта соответственно. А также значения вводимые пользователем для пороговых расстояний между центроидами, относительных значений площадей и периметров. Если все три показателя не превышают своих пороговых значений, то геометрии двух объектов будем считать похожими.
Ниже представлен SQL-запрос для поиска похожих объектов в целевой таблице public.egko__egko_building на объект с идентификатором 52278786 из исходной таблицы public.osnova_marketing_2024__building.
Поиск в другой таблице похожих объектов на данный
-- подзапрос получения исходного валидного геометрического объекта WITH source_object_valid AS ( SELECT * FROM public.osnova_marketing_2024__building WHERE id = 52278786 AND ST_IsValid(wkb_geometry) -- укажите ID вашего объекта ), -- подзапрос получения целевых валидных геометрических объектов target_object_valid AS ( SELECT * FROM public.egko__egko_building WHERE ST_IsValid(wkb_geometry) ), -- подзапрос извлечения исходных данных объекта source_object AS ( SELECT id, wkb_geometry, ST_Centroid(wkb_geometry) AS centroid, ST_Perimeter(wkb_geometry) AS perimeter, ST_Area(wkb_geometry) AS area FROM source_object_valid ), -- подзапрос вычисления метрик для перескающихся по геометрии объектов metrics_for_crossed_geoms AS ( SELECT so.id AS source_id, tb.id AS target_id, so.perimeter as source_perimeter, ST_Perimeter(tb.wkb_geometry) as target_perimeter, so.area as source_area, ST_Area(tb.wkb_geometry) as target_area, so.wkb_geometry as source_geometry, tb.wkb_geometry AS target_geometry, -- Проверка на полное равенство ST_Equals(so.wkb_geometry, tb.wkb_geometry) AS equals, -- Вычисляем метрики для сравнения -- вычисляем расстояние между центроидами ST_Distance(so.centroid, ST_Centroid(tb.wkb_geometry)) AS centroid_distance, -- вычисляем разницу между периметрами ABS(so.perimeter - ST_Perimeter(tb.wkb_geometry)) AS perimeter_diff, -- вычисляем разницу между площадями ABS(so.area - ST_Area(tb.wkb_geometry)) AS area_diff, ST_HausdorffDistance(so.wkb_geometry, tb.wkb_geometry) as HausdorffDistance, ST_FrechetDistance(so.wkb_geometry, tb.wkb_geometry) as FrechetDistance FROM source_object so JOIN target_object_valid tb -- пространственный фильтр для оптимизации ON so.wkb_geometry && tb.wkb_geometry ), -- дополнительно вычисление метрик (арифметические операции над ними, относительные занчения) metrics_range AS ( SELECT source_id, target_id, equals, centroid_distance, perimeter_diff, area_diff, source_perimeter, target_perimeter, source_area, target_area, perimeter_diff / source_perimeter as pd_sp, perimeter_diff / target_perimeter as pd_tp, area_diff / source_area as ad_sa, area_diff / target_area as ad_ta, HausdorffDistance, FrechetDistance FROM metrics_for_crossed_geoms ), pre_result AS ( SELECT source_id, target_id, equals, CASE WHEN equals THEN false WHEN centroid_distance <= 1000 -- порог расстояния между центроидами(метры) AND pd_sp<= 0.10 -- разница периметров <= 10% AND pd_tp <= 0.10 -- разница периметров <= 10% AND ad_sa <= 0.10 -- разница площадей <= 10% AND ad_ta <= 0.10 -- разница площадей <= 10% THEN true ELSE false END AS similarity, centroid_distance, perimeter_diff, source_perimeter, target_perimeter, area_diff, source_area, target_area, pd_sp, pd_tp, ad_sa, ad_ta, HausdorffDistance, FrechetDistance FROM metrics_range ), final_display AS ( SELECT *, ROUND((centroid_distance+pd_sp+pd_tp+ad_sa+ad_ta)::numeric, 2) as score FROM pre_result WHERE equals or similarity ) SELECT * FROM final_display ORDER BY score ASC
В CTE source_object_valid получаем исходный объект в том случае, если его геометрия валидна т.е. если она не имеет самопересечений, у неё не нулевая и не отрицательная площадь, и нет прочих топологических ошибок.
В CTE target_object_valid аналогично предыдущему запросу, но только для целевой таблицы.
В CTE source_object находим центроид, вычисляем площадь и периметр исходного объекта.
В CTE metrics_for_crossed_geoms вычисляем предварительные значения метрик. Используем в строке 47 пространственный фильтр && для оптимизации запроса. Каждая геометрия заключается в прямоугольник - так называемый bounding_box. И отбираются только те объекты, которые пересекаются по bounding_box. Это оказывается намного быстрее, чем искать точное пересечение при помощи функции ST_Intersects. Считаем, что объекты помимо соответствия по пороговым значениям метрик должны ещё и хоть как-то пересекаться по bounding_box.
В CTE metrics_range проводятся дополнительные вычисления над метриками.
В CTE pre_result определяем на основе соответствия метрик пороговым значениям значения equals и similarity.
В CTE final_display отбираем только одинаковые или похожие объекты. Вводим показатель похожести score как сумму всех значений метрик. Чем меньше score - тем более похожи объекты (у одинаковых объектов score равен нулю).
И наконец выводим результаты, отсортировав их по возрастанию score.
С одним исходным объектом получилось. Сразу перейдём к попарному сравнению объектов в двух таблицах.
All vs All
Необходимо внести изменения только в CTE source_object_valid, в котором будем получать все объекты из исходной таблицы.
WITH source_object_valid AS ( SELECT * FROM public.osnova_marketing_2024__building WHERE ST_IsValid(wkb_geometry) ), ...


Посмотрим на геометрии двух одинаковых объектов. Заметим, что объекты находятся в разных таблицах, но имеют совпадающие значения идентификаторов. Также обратим внимание, что их геометрии в формате WKB не совсем совпадают.
SELECT wkb_geometry FROM public.osnova_marketing_2024__building WHERE id=112 UNION ALL SELECT wkb_geometry FROM public.egko__egko_building WHERE id=112


Теперь взглянем на найденные геометрии похожих объектов и увидим отличия.
SELECT wkb_geometry FROM public.osnova_marketing_2024__building WHERE id=294 UNION ALL SELECT wkb_geometry FROM public.egko__egko_building WHERE id=294

Снова объекты находятся в двух разных таблицах, но имеют одинаковое значение идентификатора.
Поиск среди всех объектов БД.
Необходимо сделать то же самое только для 1200 таблиц, содержащих суммарно три сотни миллионов объектов. Но если осуществлять последовательный поиск по всем таблицам, то это займёт много времени. Нужна другая идея.
Был разработан следующий подход:
1) Проводим кластеризацию каждой таблицы. В результате каждый объект будет иметь составной ключ - (schema_name, table_name, cluster_id, id). Геометрии самих объектов нам не понадобятся. Эту информацию по кластеризованным объектам будем заносить в таблицу clustering_objects. Также получаем геометрии кластеров, их будет характеризовать кортеж из (schema_name, table_name, cluster_id, wkb_geometry), и добавляем их в таблицу clustering_clusters.
В этом подходе не стоит использовать алгоритм кластеризации на основе DBSCAN т.к. он выбрасывает некоторые объекты (outliers). У таких объектов cluster_id имеет значение NULL. Его имеет смысл применять для анализа данных, но не для поиска в нашем случае.
Ниже представлены запросы по созданию двух таблиц и примеры данных в них.
CREATE TABLE public.clustering_objects ( schema_name VARCHAR(255), table_name VARCHAR(255), id VARCHAR(255), cluster_id INTEGER);


CREATE TABLE public.clustering_clusters ( schema_name VARCHAR(255), table_name VARCHAR(255), cluster_id INTEGER, wkb_geometry GEOMETRY(MultiPolygon, 4326)); CREATE INDEX idx_clustering_clusters ON public.clustering_clusters USING GIST(wkb_geometry);
Кластеризация таблицы
BEGIN; WITH source_data AS ( SELECT 'osnova_marketing_2024' as schema_name, 'building' as table_name, ogc_fid as id, wkb_geometry as wkb_geometry FROM osnova_marketing_2024."building" WHERE wkb_geometry IS NOT NULL ), clustered AS ( SELECT ST_ClusterDBSCAN(wkb_geometry, eps := 100, minpoints := 50) OVER() as cluster_id, schema_name, table_name, id, wkb_geometry FROM source_data ), clusters_geom AS ( SELECT cluster_id, 'osnova_marketing_2024' as schema_name, 'building' as table_name, ST_Transform(ST_ConvexHull(ST_Collect(wkb_geometry)), 4326) as wkb_geometry FROM clustered WHERE cluster_id IS NOT NULL GROUP BY cluster_id ), -- Первая вставка в таблицу с кластерами insert_clusters AS ( INSERT INTO public.clustering_clusters ( schema_name, table_name, cluster_id, wkb_geometry ) SELECT schema_name, table_name, cluster_id, wkb_geometry FROM clusters_geom RETURNING 1 ) -- Вторая вставка в таблицу объектов INSERT INTO public.clustering_objects ( schema_name, table_name, id, cluster_id ) SELECT schema_name, table_name, id, cluster_id FROM clustered WHERE cluster_id IS NOT NULL; COMMIT;
2) Пересекаем исходный объект с кластерами из таблицы clustering_clusters. Тем самым находим записи с тройкой (schema_name, table_name, cluster_id). Например, теперь мы знаем, что искомые объекты находятся в 20-и кластерах 10-и таблиц.
Скрытый текст
WITH source_object AS ( SELECT ogc_fid as id, wkb_geometry as wkb_geometry FROM _25818._25835 LIMIT 1 ), suitable_clusters AS ( SELECT ac.* FROM clustering_clusters as ac, source_object as so WHERE ac.wkb_geometry && so.wkb_geometry ) SELECT * FROM suitable_clusters -- отправить на поиск похожих объектов с учётом полного имени таблицы и кластера
3) Отправляем исходную геометрию на поиск одинаковых и похожих объектов среди выявленных объектов. То есть переходим уже к самим потенциальным объектам в их исходных таблицах т.к. знаем тройку (schema_name, table_name, cluster_id).
В общем этот подход основан на постепенном сужении области поиска. Но вначале потребуется создать модель кластеризации проведя кластеризацию всех таблиц БД. Это займёт 3-4 часа, что вполне приемлемо. При осуществлении изменений (вставка, удаление, изменение геометрии) в какой-либо одной таблице нужно будет лишь произвести кластеризацию конкретной этой таблицы и внести соответствующие изменения в таблицы clustering_tables и clustering_clusters. Например, это можно осуществить с помощью триггеров.
Заключение
В этой статье были изложены идеи и описан функционал системы для работы с геометрией объектов. Таким образом, мы используем кластеризацию не только для анализа географических объектов, но также и для сокращения времени поисковых запросов.
В следующей части применим тематическое моделирование для анализа наших данных и полнотекстовый поиск для нахождения объектов по их описанию.
Спасибо за внимание!
