PostGis. Как найти ошибку в пространственном запросе?

Добрый день! Я — Виктор, разработчик в 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); 

Отвечу на возникающий вопрос: «почему нельзя исправить всю ошибочную геометрию в таблице, чтобы не искать выборочно причины?».

Дело в том, что пространственные данные поступают к нам в систему из различных источников (в т.ч. из Росреестра) и мы не можем выполнять исправление (как правило, оно сопровождается искажением) всех данных. Получив нужные ключи, мы анализируем какие данные они представляют и можно ли их исправить.

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

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

А как часто вы сталкиваетесь с проблемами геометрии и как обеспечиваете качество пространственных данных?