Jak najlepiej naprawić problem skrzyżowania bez węzłów w PostGIS?

Używam PL/Rfunkcji i PostGISdo generowania wielokątów voronoi wokół zestawu punktów. Funkcja, której używam, jest tutaj zdefiniowana . Gdy korzystam z tej funkcji w określonym zestawie danych, pojawia się następujący komunikat o błędzie:

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

Po zbadaniu tej części komunikatu o błędzie:

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

Tak wygląda powyższy problem:

Początkowo myślałem, że ten komunikat może być spowodowany istnieniem identycznych punktów, i próbowałem rozwiązać to za pomocą st_translate()funkcji użytej w następujący sposób:

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

To rozwiązuje problem, ale obawiam się, że tłumaczę teraz wszystkie punkty do ~ 20 m w kierunku x / y. Nie umiem też powiedzieć, jaka jest odpowiednia kwota tłumaczenia. Na przykład w tym zestawie danych metodą prób i błędów a 20m * random numberjest w porządku, ale skąd mam wiedzieć, czy to musi być większe?

Na podstawie powyższego obrazu myślę, że problem polega na tym, że punkt przecina się z linią, podczas gdy algorytm próbuje przeciąć punkt wielokątem. Nie jestem pewien, co powinienem zrobić, aby upewnić się, że punkt znajduje się w wielokącie, zamiast przecinać się z linią. Błąd występuje w tym wierszu:

"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'');"

Przeczytałem poprzednie pytanie: Co to jest „skrzyżowanie bez węzła”? starając się lepiej zrozumieć ten problem i chętnie skorzystamy z wszelkich porad, jak najlepiej go rozwiązać.

Z mojego doświadczenia wynika, że ​​ten problem prawie zawsze jest spowodowany:

1. Wysoka precyzja współrzędnych (43.231499999999996) w połączeniu z
2. Linie, które są prawie zbieżne, ale nie identyczne

Podejście „podchwytliwe” ST_Bufferrozwiązań pozwala ci uciec się do punktu 2, ale wszystko, co możesz zrobić, aby rozwiązać te podstawowe przyczyny, takie jak przyciąganie geometrii do siatki 1e-6, ułatwi ci życie. Buforowane geometrie są zwykle odpowiednie do obliczeń pośrednich, takich jak obszar nakładania się, ale należy zachować ostrożność, aby je zachować, ponieważ mogą pogorszyć bliskie, ale niezupełnie problemy na dłuższą metę.

Funkcja obsługi wyjątków PostgreSQL pozwala pisać funkcje otoki w celu obsługi tych specjalnych przypadków, buforując tylko w razie potrzeby. Oto przykład ST_Intersection; Używam podobnej funkcji dla ST_Difference. Musisz zdecydować, czy buforowanie i potencjalny zwrot pustego wielokąta są dopuszczalne w twojej sytuacji.

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;

Inną zaletą tego podejścia jest to, że możesz wskazać geometrie, które faktycznie powodują problemy; wystarczy dodać kilka RAISE NOTICEzdań w EXCEPTIONbloku, aby wyświetlić WKT lub coś innego, co pomoże wyśledzić problem.

Przez wiele prób i błędów w końcu zdałem sobie sprawę, że non-noded intersectionwynikało to z problemu skrzyżowania. Znalazłem wątek, który sugeruje użycie ST_buffer(geom, 0)może być użyty do rozwiązania problemu (chociaż ogólnie sprawia, że ​​jest znacznie wolniejszy). Następnie próbowałem użyć ST_MakeValid()i po zastosowaniu bezpośrednio do geometrii przed jakąkolwiek inną funkcją. Wydaje się, że to zdecydowanie rozwiązuje problem.

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]
        )
    )

Oznacziłem to jako odpowiedź, ponieważ wydaje się, że jest to jedyne podejście, które rozwiązuje mój problem.

Zetknąłem się z tym samym problemem (Postgres 9.1.4, PostGIS 2.1.1) i jedyną rzeczą, która działała dla mnie, było owinięcie geometrii bardzo małym buforem.

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

ST_MakeValidnie działało dla mnie, podobnie jak kombinacja ST_Nodei ST_Dump. Bufor nie powodował żadnego pogorszenia wydajności, ale jeśli go zmniejszyłem, nadal pojawiał się błąd skrzyżowania bez węzłów.

Brzydkie, ale działa.

Aktualizacja:

Wydaje się, że strategia ST_Buffer działa dobrze, ale natknąłem się na problem, który powodował błędy podczas rzutowania geometrii na geografię. Na przykład, jeśli punkt ma pierwotnie wartość -90,0 i jest buforowany o 0,0000001, to teraz ma wartość -90,0000001, co jest nieprawidłową lokalizacją geograficzną.

Oznaczało to, że choć ST_IsValid(geom)było t, ST_Area(geom::geography)wrócił NaNdo wielu funkcji.

Aby uniknąć problemu skrzyżowania bez węzłów, przy zachowaniu prawidłowej lokalizacji geograficznej, skończyło się na tym, ST_SnapToGridże tak

SELECT ST_Union(ST_MakeValid(ST_SnapToGrid(geom, 0.0001))) AS geom, common_id
    FROM table
    GROUP BY common_id;

W Postgis ST_Node powinien przerwać serię linii na skrzyżowaniach, co powinno rozwiązać problem skrzyżowania bez węzłów. Zawijanie tego w ST_Dump generuje złożony układ linii przerywanych.

Nieco spokrewniona, jest niesamowita prezentacja PostGIS: Wskazówki dla zaawansowanych użytkowników, która wyraźnie przedstawia tego rodzaju problemy i rozwiązania.

Z mojego doświadczenia wynika , że mój non-noded intersectionbłąd rozwiązałem za pomocą funkcji St_SnapToGrid, która rozwiązała problem wysokiej precyzji współrzędnych wierzchołka wielokątów.

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';