Добрый день! Я — Виктор, разработчик в Gems development. Ежедневно наша команда работает с пространственными данными разной сложности и качества. При выполнении операции пространственного пересечения с помощью Postgis в СУБД Postgresql мы столкнулись со следующей ошибкой:
XX000: GEOSIntersects: TopologyException: side location conflict at 10398.659 3844.9200000000001Запрос, приводящий к ошибке, выглядит так:
select q1.key,st_asGeoJson(geoloc)
from usahalinsk.V_GEO_OOPT q1
where ST_Intersects(geoloc,
ST_GeomFromGeoJSON('{"type":"Polygon","coordinates":
[[[11165.15,2087.5],[11112,2066.6],[11127.6,2022.5],
[11122.6,2020.7],
[11122.25,2021.2],[11107.07,2015.7],
[11121,1947],[11123.48,1922.99],[11128.42,1874.4],
[11131.5,1875],[11140.96,1876.81],[11160.73,1880.59],
[11201.04,1888.3],[11194.2,1908],[11221.93,1916.57],
[11223.3,1917],[11165.15,2087.5]]]}'))
Решение этой проблемы блокирует работу пользователей, так как не позволяет строить отчёты по данным и замедляет работу по оказанию услуг. Многие действия в разрабатываемой нами системе, такие как: подготовка схемы расположения земельного участка, подготовка градостроительного плана земельного участка и другие, используют пространственные операции подобные этой.
Предположим, что дело в некорректной геометрии. Часто подобную ошибку выдает операция пересечения, если участвующие в запросе объекты имеют самопересечения или дубли точек. Пример данных ошибок в геометрии можно увидеть ниже. (Граница полигона пересекает сама себя и в линии имеется две одинаковые координаты)
Мы провели собственное расследование по поиску причин ошибки и хотим рассказать об этом.
В настоящий момент мы используем Postgis 2.4 и Postgresql 9.6. Перейдем сразу к практике. Проверим константную геометрию на валидность и находим, что все работает корректно.
Можно предположить, что дело в таблице (представлении) usahalinsk.V_GEO_OOPT в котором мы ищем пересечения. Для того, чтобы подтвердить гипотезу, проверим эти данные тоже.
Но и здесь ошибок не находим. Кроме того, данные в выборку не попали вообще. Если бы они были, то задача решалась бы исправлением найденных записей через функцию Postgis st_makeValid.
Но ошибок в представлении нет, а запрос не выполняется. Предлагаем посмотреть его план.
Примечание: в реальной модели мы используем три столбца для геометрии (для полигонов, линий и точек), но для краткости будем называть это полем geoloc – в нем хранится геометрия и оно выводится в представление.
Наше представление usahalinsk.V_GEO_OOPT построено как выборка из таблицы с пространственными данными usahalinsk.d_geometry и по полю с геометрией создан пространственный индекс.
Значит, при выполнении запроса идет чтение индекса и где-то в таблице, не попадая в нашу выборку, есть невалидные пространственные данные, которые попали в индекс, т.к. он построен по всей таблице.
Давайте попробуем удалить индекс:
DROP INDEX usahalinsk.d_geometry_cs1_all_sx;И попробуем выполнить проблемный запрос.
Он выполнился без ошибок. Подтверждаем, что дело в индексе. Можно вернуть индекс, но с условием на корректную геометрию:
CREATE INDEX d_geometry_cs1_all_sx
ON usahalinsk.d_geometry
USING gist(geoloc)
where st_isvalid(geoloc)=true;
Проверим выполнение и посмотрим план.
Запрос выполнился без ошибок, индекс в плане также используется. Из минусов такого решения может быть замедление вставки/обновления, т.к. дополнительно будет проверяться условие при перестроении индекса.
Вернем это изменение назад и попробуем всё-таки найти из-за каких объектов в индексе наш запрос не может выполниться.
DROP INDEX usahalinsk.d_geometry_cs1_all_sx;
CREATE INDEX d_geometry_cs1_all_sx
ON usahalinsk.d_geometry
USING gist
(geoloc);
Напомню, что у нас есть координаты места ошибки:
XX000: GEOSIntersects: TopologyException: side location conflict at 10398.659 3844.9200000000001Но если мы поищем в данных или в результате функции IsValidReason, которая выдает причину ошибки, то ничего похожего не найдем.
select key,ST_IsValidReason(geoloc)
from usahalinsk.d_geometry
where st_isvalid(geoloc)!=true
and ST_AsText(geoloc) like '%3844.9200000000001%';
select key,ST_IsValidReason(geoloc)
from usahalinsk.d_geometry
where st_isvalid(geoloc)!=true
and ST_IsValidReason(geoloc) like '%3844.9200000000001%';
Чтобы найти объекты, которые влияют на запрос, можно использовать следующий скрипт. Будем проверять каждый объект таблицы и пересекать его с искомой константой. В ходе выполнения ловим исключения и проверяем их содержимое. Если ошибка будет содержать нужные нам координаты, значит это наша проблемная геометрия.
do
$$
declare
tKey bigint;
rec record;
error_text text;
--Тест ошибки
error_info text:='GEOSIntersects: TopologyException: side location conflict at 10398.659 3844.9200000000001';
begin
--Перебираем все данные в таблице
for rec in(select key from usahalinsk.d_geometry)
loop
begin
select key into tKey
from (select * from usahalinsk.d_geometry q1
--сравнение по первичному ключу
where q1.key=rec.key
and ST_Intersects(geoloc,
--константная геометрия
ST_GeomFromGeoJSON('{"type":"Polygon","coordinates":[[[11165.15,2087.5],
[11112,2066.6],[11127.6,2022.5],[11122.6,2020.7],
[11122.25,2021.2],[11107.07,2015.7],[11121,1947], [11123.48,1922.99],[11128.42,1874.4],
[11131.5,1875],[11140.96,1876.81], [11160.73,1880.59],[11201.04,1888.3],
[11194.2,1908],[11221.93,1916.57],[11223.3,1917],
[11165.15,2087.5]]]}'))) geoQ;
exception when others then
--получаем ошибку если она есть
GET STACKED DIAGNOSTICS error_text = MESSAGE_TEXT;
--Если её тест равен искомому, то выводим ключ
if error_text=error_info then
raise info '%',rec.key;
end if;
end;
end loop;
end$$;
В результате получаем три ключа геометрии, которые легко исправить:
update usahalinsk.d_geometry
set cs1_geometry_polygone=st_collectionextract(st_makevalid(geoloc),3)
where key in(
1000010001988961,
1000010001989399,
1000010004293508);
Отвечу на возникающий вопрос: «почему нельзя исправить всю ошибочную геометрию в таблице, чтобы не искать выборочно причины?».
Дело в том, что пространственные данные поступают к нам в систему из различных источников (в т.ч. из Росреестра) и мы не можем выполнять исправление (как правило, оно сопровождается искажением) всех данных. Получив нужные ключи, мы анализируем какие данные они представляют и можно ли их исправить.
Тривиальная задача поиска причины ошибки может превратиться в целое расследование со скриптом исправления в конце.
Более сложный вариант задачи: как быть, если выполняется пересечение не с константой, а с другой таблицей? Как вариант, выполнять пересечение каждого из участвующих объектов первой таблицы с каждым объектом второй. И отлавливать исключения.
А как часто вы сталкиваетесь с проблемами геометрии и как обеспечиваете качество пространственных данных?
