Как лучше всего исправить проблему отсутствия пересечения в PostGIS?

38

Я использую PL/Rфункцию и PostGISдля генерации вороных многоугольников вокруг множества точек. Функция, которую я использую, определяется здесь . Когда я использую эту функцию для определенного набора данных, я получаю следующее сообщение об ошибке:

Error : ERROR:  R interpreter expression evaluation error
DETAIL:  Error in pg.spi.exec(sprintf("SELECT %3$s AS id,   
st_intersection('SRID='||st_srid(%2$s)||';%4$s'::text,'%5$s') 
AS polygon FROM %1$s WHERE st_intersects(%2$s::text,'SRID='||st_srid(%2$s)||';%4$s');",  
:error in SQL statement : Error performing intersection: TopologyException: found non-noded 
intersection between LINESTRING (571304 310990, 568465 264611) and LINESTRING (568465 
264611, 594406 286813) at 568465.05533706467 264610.82749605528
CONTEXT:  In R support function pg.spi.exec In PL/R function r_voronoi

Из изучения этой части сообщения об ошибке:

Error performing intersection: TopologyException: found non-noded intersection between
LINESTRING (571304 310990, 568465 264611) and LINESTRING (568465 264611, 594406 286813) 
at 568465.05533706467 264610.82749605528

Вот как выглядит проблема, перечисленная выше:

введите описание изображения здесь

Сначала я подумал, что это сообщение может быть вызвано наличием идентичных точек, и попытался решить эту проблему с помощью st_translate()функции, используемой следующим образом:

ST_Translate(geom, random()*20, random()*20) as geom 

Это решает проблему, но я обеспокоен тем, что теперь я перевожу все точки на расстояние до ~ 20 м в направлении x / y. Я также не могу сказать, какая сумма перевода нужна. Например, в этом наборе данных методом проб и ошибок все в 20m * random numberпорядке, но как я могу определить, должно ли это быть больше?

Основываясь на изображении выше, я думаю, что проблема в том, что точка пересекается с линией, а алгоритм пытается пересечь точку с многоугольником. Я не уверен, что мне следует делать, чтобы точка находилась внутри многоугольника, а не пересекалась с линией. Ошибка происходит в этой строке:

"SELECT 
  %3$s AS id, 
  st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'') AS polygon 
FROM 
  %1$s 
WHERE 
  st_intersects(%2$s::text,''SRID=''||st_srid(%2$s)||'';%4$s'');"

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

djq
источник
Если ваши входные данные недопустимы для начала, запустите для них ST_MakeValid (). Если они действительны, добавление энтропии, как вы делаете, является следующей доступной уловкой и, возможно, последней уловкой на данный момент.
Пол Рэмси
Да, я использую WHERE ST_IsValid(p.geom)для фильтрации точек изначально.
DJQ

Ответы:

30

По моему опыту, эта проблема почти всегда вызвана:

  1. Высокая точность в ваших координатах (43.231499999999996), в сочетании с
  2. Линии, которые почти совпадают, но не идентичны

Подход «подталкивания» ST_Bufferрешений позволяет вам избежать проблем с №2, но все, что вы можете сделать для устранения этих основных причин, например привязка вашей геометрии к сетке 1е-6, сделает вашу жизнь проще. Буферизованная геометрия обычно подходит для промежуточных вычислений, таких как область перекрытия, но вы должны быть осторожны с их сохранением, потому что они могут усугубить ваши проблемы, но не совсем, в долгосрочной перспективе.

Возможность обработки исключений в PostgreSQL позволяет вам писать функции-оболочки для обработки этих особых случаев, буферизуя только при необходимости. Вот пример для ST_Intersection; Я использую аналогичную функцию для ST_Difference. Вам нужно решить, приемлемы ли буферизация и потенциальный возврат пустого многоугольника в вашей ситуации.

CREATE OR REPLACE FUNCTION safe_isect(geom_a geometry, geom_b geometry)
RETURNS geometry AS
$$
BEGIN
    RETURN ST_Intersection(geom_a, geom_b);
    EXCEPTION
        WHEN OTHERS THEN
            BEGIN
                RETURN ST_Intersection(ST_Buffer(geom_a, 0.0000001), ST_Buffer(geom_b, 0.0000001));
                EXCEPTION
                    WHEN OTHERS THEN
                        RETURN ST_GeomFromText('POLYGON EMPTY');
    END;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;

Еще одним преимуществом этого подхода является то, что вы можете точно определить геометрию, которая на самом деле вызывает ваши проблемы; просто добавьте несколько RAISE NOTICEоператоров в EXCEPTIONблок для вывода WKT или что-то еще, что поможет вам отследить проблему.

dbaston
источник
Это умная идея. Я часто обнаруживал, что проблемы пересечения возникают из-за линий линий, возникающих во время сочетаний союзов, различий, буферов и т. Д., Которые могут быть исправлены путем буферизации всего или сброса всего и только выбора полигонов / многолигонов. Это интересный подход.
Джон Пауэлл
Вы упоминаете привязку геометрии к сетке 1e-6, но я сижу здесь и думаю, будет ли лучше привязка к степени 2. PostGIS (и GEOS) используют числа с плавающей запятой, поэтому привязка к степени 10 может на самом деле не сильно урезать координаты, поскольку число может не иметь двоичного представления конечной длины. Но если вы скажете 2 ^ -16, я думаю, что будет гарантированно обрезать любую дробную часть до 2 байтов. Или я не так думаю?
jpmc26
12

Через много проб и ошибок я в конечном итоге понял, что это non-noded intersectionстало результатом проблемы самопересечения. Я обнаружил, что ST_buffer(geom, 0)для решения проблемы можно использовать ветку, в которой предлагается использовать (хотя в целом это делает ее намного медленнее). Затем я попытался использовать ST_MakeValid()и при применении непосредственно к геометрии перед любой другой функцией. Это, кажется, решает проблему надежно.

ipoint <- pg.spi.exec(
        sprintf(
            "SELECT 
                    %3$s AS id, 
                    st_intersection(ST_MakeValid(''SRID=''||st_srid(%2$s)||'';%4$s''::text), ST_MakeValid(''%5$s'', 0)) AS polygon 
            FROM %1$s 
            WHERE 
                ST_Intersects(ST_MakeValid(%2$s::text),ST_MakeValid(''SRID=''||st_srid(%2$s)||'';%4$s''));",
            arg1,
            arg2,
            arg3,
            curpoly,
            buffer_set$ewkb[1]
        )
    )

Я пометил это как ответ, так как кажется, что это единственный подход, который решает мою проблему.

djq
источник
11

Я столкнулся с этой же проблемой (Postgres 9.1.4, PostGIS 2.1.1), и единственное, что мне помогло, - это обернуть геометрию очень маленьким буфером.

SELECT ST_Intersection(
    (SELECT geom FROM table1), ST_Union(ST_Buffer(geom, 0.0000001))
) FROM table2

ST_MakeValidне работает для меня, а также не комбинация ST_Nodeи ST_Dump. Похоже, что буфер не приводит к ухудшению производительности, но если я уменьшу его, я все равно получу ошибку пересечения без узлов.

Уродливо, но это работает.

Обновить:

Стратегия ST_Buffer, кажется, работает хорошо, но я столкнулся с проблемой, когда она приводила к ошибкам при приведении геометрии к географии. Например, если точка изначально имеет значение -90,0 и буферизуется значением 0,0000001, теперь она имеет значение -90,0000001, что является недопустимой географией.

Это означало , что даже если ST_IsValid(geom)было t, ST_Area(geom::geography)вернулись NaNдля многих функций.

Чтобы избежать проблемы пересечения без узлов, сохраняя правильную географию, я в конечном итоге использовал ST_SnapToGridвот так

SELECT ST_Union(ST_MakeValid(ST_SnapToGrid(geom, 0.0001))) AS geom, common_id
    FROM table
    GROUP BY common_id;
jczaplew
источник
6

В postgis ST_Node должен разбивать ряд линий на пересечениях, что должно решить проблему пересечения без узлов. Обтекание этого в ST_Dump генерирует составной массив из пунктирных линий.

Немного связано, есть отличная презентация PostGIS: Советы для опытных пользователей, в которой четко обозначены проблемы и решения такого рода.

WolfOdrade
источник
Это отличная презентация (спасибо @PaulRamsey). Как я должен использовать ST_Nodeи ST_Dump? Я полагаю, что мне нужно использовать их рядом с этой частью функции, но я не уверен: st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'')in
djq
Хм, я не заметил, что две линии имеют одинаковую координату, что должно быть хорошо. Если вы нанесете эти координаты, точка пересечения находится на расстоянии около 18 см от пересечения. Не совсем решение, просто наблюдение.
WolfOdrade
Не совсем понятно, как я использую st_nodeздесь - могу ли я использовать его раньше st_intersection?
DJQ
1
Презентация больше не доступна. Я застрял с той же проблемой, когда пытаюсь ST_Clip (раст, полигон)
Джеки
1
@Jackie: я исправил ссылку на презентацию в ответе: PostGIS: Советы для опытных пользователей .
Пит
1

Исходя из своего опыта, я решил свою non-noded intersectionошибку, используя функцию St_SnapToGrid, которая решила проблему с высокой точностью в координатах вершины многоугольников.

SELECT dissolve.machine, dissolve.geom FROM (
        SELECT machine, (ST_Dump(ST_Union(ST_MakeValid(ST_SnapToGrid(geom,0.000001))))).geom 
        FROM cutover_automatique
        GROUP BY machine) as dissolve
WHERE ST_isvalid(dissolve.geom)='t' AND ST_GeometryType(dissolve.geom) = 'ST_Polygon';
CptGasse
источник