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

38

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:

wprowadź opis zdjęcia tutaj

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ć.

djq
źródło
Jeśli twoje dane wejściowe nie są poprawne na początek, uruchom na nich ST_MakeValid (). Jeśli są poprawne, dodanie entropii, tak jak robisz, jest kolejną dostępną lewą i być może ostatnią lewą na razie.
Paul Ramsey,
Tak, WHERE ST_IsValid(p.geom)początkowo używam do filtrowania punktów.
djq

Odpowiedzi:

30

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.

dbaston
źródło
To sprytny pomysł. Często stwierdziłem, że problemy z przecięciem wynikają z linii pojawiających się podczas kombinacji związków, różnic, buforów itp., Które można naprawić, buforując wszystko lub zrzucając wszystko i wybierając tylko Wielokąty / Mutlipolygons. To ciekawe podejście.
John Powell,
Wspomniałeś przyciąganie geometrii do siatki 1e-6, ale siedzę tutaj, zastanawiając się, czy przyciąganie do potęgi 2 byłoby lepsze. PostGIS (i GEOS) przy użyciu liczb zmiennoprzecinkowych, więc przyciąganie do potęgi 10 może w rzeczywistości nie bardzo obciąć współrzędne, ponieważ liczba może nie mieć reprezentacji binarnej o skończonej długości. Ale jeśli zgodzisz się powiedzieć 2 ^ -16, uważam, że gwarantowałoby to skrócenie dowolnej części ułamkowej do zaledwie 2 bajtów. A może źle myślę?
jpmc26
12

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.

djq
źródło
11

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;
jczaplew
źródło
6

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.

WolfOdrade
źródło
To świetna prezentacja (dzięki @PaulRamsey). Jak powinienem używać ST_Nodei ST_Dump? Wyobrażam sobie, że musiałbym użyć ich w pobliżu tej części funkcji, ale nie jestem pewien: st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'')in
djq
Hmmm nie zauważyłem, że dwie linie mają identyczną współrzędną, co powinno być w porządku. Jeśli narysujesz te współrzędne, punkt przecięcia znajduje się około 18 cm od skrzyżowania. Nie do końca rozwiązanie, tylko obserwacja.
WolfOdrade
Nie do końca jasne, jak st_nodetu korzystam - czy mogę go używać wcześniej st_intersection?
djq
1
Prezentacja nie jest już dostępna. Utknąłem z tym samym problemem, gdy próbuję ST_Clip (rast, wielokąt)
Jackie
1
@Jackie: Naprawiłem link do prezentacji w odpowiedzi: PostGIS: Wskazówki dla zaawansowanych użytkowników .
Pete,
1

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';
CptGasse
źródło