Jak tworzyć linie, aby wizualizować różnice między elementami wielokątów w PostGIS?

15

Mam tabelę PostGIS polygon_bz kilkoma funkcjami wielokąta. Istnieje również tabela, polygon_aktóra zawiera te same wielokąty, polygon_bale z niewielkimi zmianami. Teraz chcę utworzyć linie w celu wizualizacji różnic między elementami wielokąta.

wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj wprowadź opis zdjęcia tutaj

Przypuszczam, że ST_ExteriorRingi ST_Differencebędzie wykonać zadanie, ale WHERE wydaje się być dość trudne.

CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, yourSRID) AS geom
FROM 
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE 
    -- ?
    ) AS g;

Czy ktoś może mi pomóc?

EDYCJA 1

Jak napisałem przez „tilt”, próbowałem, ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom)ale wynik nie jest zgodny z oczekiwaniami.

CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, your_SRID) AS geom
FROM 
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE 
    ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom))
     AS g;

wprowadź opis zdjęcia tutaj

EDYCJA 2

workupload.com/file/J0WBvRBb (przykładowy zestaw danych)


Próbowałem zamienić wielokąty w multilinie przed użyciem ST_Difference, ale wyniki są nadal dziwne.

CREATE VIEW multiline_a AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_a.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_a;

CREATE VIEW multiline_b AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_b.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_b;

CREATE VIEW line_difference AS SELECT
row_number() over() as gid,
g.geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(multiline_a.geom, multiline_b.geom)))).geom::geometry(linestring, 4326) AS geom
    FROM
    multiline_a, multiline_b)
As g;

wprowadź opis zdjęcia tutaj

Morze Księżycowe
źródło
Wygląda bardziej jak pytanie topologiczne. Chcesz zidentyfikować segmenty, które nie są objęte inną warstwą. Nie pracowałem zbyt wiele z topologią PostGIS i nie mogę udzielić bezpośredniej odpowiedzi, ale sugeruję przyjrzeć się temu.
Thomas
Ciekawe, czy masz przykładowy zestaw danych do pobrania?
huckfinn

Odpowiedzi:

10

Oto kilka nowych sztuczek, wykorzystujących:

  • EXCEPTaby usunąć geometrie z obu tabel, które są takie same, abyśmy mogli skupić się tylko na geometriach unikalnych dla każdej tabeli ( A_onlyi B_only).
  • ST_Snap aby uzyskać dokładne węzły dla operatorów nakładek.
  • Użyj ST_SymDifferenceoperatora nakładki, aby znaleźć różnicę symetryczną między dwoma zestawami geometrii, aby pokazać różnice. Aktualizacja : ST_Differencepokazuje ten sam wynik dla tego przykładu. Możesz wypróbować dowolną funkcję, aby zobaczyć, co otrzymują.

To powinno uzyskać to, czego oczekujesz:

-- CREATE OR REPLACE VIEW polygon_SymDifference AS
SELECT row_number() OVER () rn, *
FROM (
  SELECT (ST_Dump(ST_SymDifference(ST_Snap(A, B, tol), ST_Snap(B, A, tol)))).*
  FROM (
    SELECT ST_Union(DISTINCT A_only.geom) A, ST_Union(DISTINCT B_only.geom) B, 1e-5 tol
    FROM (
      SELECT ST_Boundary(geom) geom FROM polygon_a
      EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b
    ) A_only,
    (
      SELECT ST_Boundary(geom) geom FROM polygon_b
      EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_a
    ) B_only
  ) s
) s;

 rn |                                        geom
----+-------------------------------------------------------------------------------------
  1 | LINESTRING(206.234028204842 -92.0360704110685,219.846021625456 -92.5340701703592)
  2 | LINESTRING(18.556700448873 -36.4496098325257,44.44438533894 -40.5104231486146)
  3 | LINESTRING(-131.974995802602 -38.6145334122719,-114.067738329597 -39.0215165366584)
(3 rows)

trzy linie


Aby jeszcze bardziej rozpakować tę odpowiedź, pierwszym krokiem ST_Boundaryjest określenie granicy każdego wielokąta, a nie tylko jego zewnętrznej powierzchni. Na przykład, jeśli byłyby dziury, byłyby one śledzone przez granicę.

EXCEPTWarunek służy do usuwania geometrii z A, które stanowią część B i rzędów z B, które stanowią część A. zmniejsza liczbę wierszy, które są tylko część A i część B tylko. Na przykład, aby uzyskać A_only:

SELECT ST_Boundary(geom) geom FROM polygon_a
EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b

Oto 6 wierszy A_only i 3 wiersze B_only: Tylko B_only

Następnie ST_Union(DISTINCT A_only.geom)służy do łączenia linii w jedną geometrię, zazwyczaj MultiLineString.

ST_Snap służy do przyciągania węzłów z jednej geometrii do drugiej. Na przykład ST_Snap(A, B, tol)weźmie geometrię A i doda więcej węzłów z geometrii B lub przeniesie je do geometrii B, jeśli znajdują się w tolodległości. Prawdopodobnie istnieje kilka sposobów korzystania z tych funkcji, ale chodzi o uzyskanie współrzędnych z każdej geometrii, które są do siebie dokładne. Dwie geometrie po przyciągnięciu wyglądają tak:

Warknął B warknął

I pokazać różnice, można wybrać albo ST_SymDifferencealbo ST_Difference. Oba pokazują ten sam wynik dla tego przykładu.

Mike T.
źródło
Niezła odpowiedź. Zastanawiałem się, co wykorzystałeś do wizualizacji wyników swoich zapytań pośrednich. Nie od razu wyglądał jak qgis, a może jest to coś, co renderuje się trochę szybciej?
RoperMaps
1
Używam JTS Testbuilder do przeglądania i przetwarzania geometrii. Jest to powiązany silnik geometrii z GEOS i Shapely, ale ma graficzny interfejs użytkownika oparty na Javie.
Mike T
Czy istnieje sposób na zignorowanie / pominięcie problemów „Nie-węzłowe skrzyżowanie między LINESTRING”? Wszystkie wejściowe wielokąty wydają się być w porządku (sprawdzane za pomocą narzędzia sprawdzania geometrii QGIS).
eclipsed_by_the_moon
1
„ST_Boundary (ST_SnapToGrid (geom, 0,001))” zamiast „ST_Boundary (geom)” rozwiązuje problem.
eclipsed_by_the_moon
6

Myślę, że jest to trochę skomplikowane, z powodu różnych zestawów węzłów obu twoich wielokątów (zielony wielokąt A, czerwony różne segmenty polionu B). Porównanie segmentów obu wielokątów daje wskazówkę, które segmenty wielokąta B zostaną zmodyfikowane.

Węzły wielokąta A

poli a

Węzły „różnych” segmentów wielokąta B

seg diff

Niestety pokazuje to tylko różnicę w strukturze segmentu, ale mam nadzieję, że jest to punkt wyjścia i działa tak:

Po pobraniu i rozpakowaniu zaimportowałem zestaw danych za pomocą PostgrSQL 9.46, PostGIS 2.1 pod Debian Linux Jessie z poleceniami.

$ createdb gis-se
$ psql gis-se < /usr/share/postgis-2.1/postgis.sql
$ psql gis-se < /usr/share/postgis-2.1/spatial_ref_sys.sql
$ shp2pgsql -S polygon_a | psql gis-se
$ shp2pgsql -S polygon_b | psql gis-se

Zakładając, że segmentów wielokąta A nie ma w B i vice vera, staram się budować różnicę między segmentami obu zestawów wielokątów, pomijając przynależność segmentu do wielokątów w każdej grupie (A lub B). Z powodów dydaktycznych formułuję SQL w kilku widokach.

Odpowiadając temu postowi GIS-SE , rozkładam oba wielokąty na tabele segmentów segments_aisegments_b

-- Segments of the polygon A
CREATE VIEW segments_a AS SELECT sp, ep
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
    -- extract the individual linestrings
     (SELECT (ST_Dump(ST_Boundary(geom))).geom
      FROM polygon_a
     ) AS linestrings
    -- be sure that nothing is scrambled
    ORDER BY sp, ep
) AS segments;

Wielokąt A tabeli segmentów:

SELECT 
  st_astext(sp) AS sp, 
  st_astext(ep) AS ep 
FROM segments_a 
LIMIT 3;
                    sp                     |                 ep
-------------------------------------------+--------------------------------------------
POINT(-292.268907321861 95.0342877387557)  | POINT(-287.118411917425 99.4165242769195)
POINT(-287.118411917425 99.4165242769195)  | POINT(-264.62129248575 93.2470010145007)
POINT(-277.459563916327 -44.5629543976138) | POINT(-292.268907321861 95.03428773875

Tę samą procedurę zastosowano do wielokąta B.

-- Segments of the polygon B
CREATE VIEW segments_b AS SELECT sp, ep
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
    -- extract the individual linestrings
     (SELECT (ST_Dump(ST_Boundary(geom))).geom
      FROM polygon_b
     ) AS linestrings
    -- be sure that nothing is scrambled
    ORDER BY sp, ep
) AS segments;

Wielokąt B tabeli segmentów

SELECT
  st_astext(sp) AS sp, 
  st_astext(ep) AS ep 
FROM segments_b 
LIMIT 3;
                    sp                     |                    ep
-------------------------------------------+-------------------------------------------
POINT(-292.268907321861 95.0342877387557)  | POINT(-287.118411917425 99.4165242769195)
POINT(-287.118411917425 99.4165242769195)  | POINT(-264.62129248575 93.2470010145007)
POINT(-277.459563916327 -44.5629543976138) | POINT(-292.268907321861 95.0342877387557)
...                        

Mogę zbudować widok tabeli różnic o nazwie segments_diff_{a,b}. Różnica wynika z braku posortowanych punktów początkowych lub końcowych w zestawie segmentów A i B.

CREATE VIEW segments_diff_a AS
SELECT st_makeline(b.sp, b.ep) as geom
FROM segments_b as b
LEFT JOIN segments_a as a ON (a.sp=b.sp and a.ep = b.ep)
-- filter segments without corresponding stuff in polygon A
WHERE a.sp IS NULL;

segs diff b

I rzeczy uzupełniające:

CREATE VIEW segments_diff_b AS
SELECT st_makeline(a.sp, a.ep) as geom
FROM segments_a as a
LEFT JOIN segments_b as b ON (a.sp=b.sp and a.ep = b.ep)
-- filter segments without corresponding stuff in polygon B
WHERE b.sp IS NULL;

segs diff a

Wniosek: Aby uzyskać właściwy wynik dla małych małych segmentów oznaczonych czerwoną strzałką, oba wielokąty muszą mieć taką samą strukturę węzła i wymagany jest krok przecięcia na poziomie węzła (wstawienie wierzchołków wielokąta A w B). Skrzyżowanie może być wykonane przez:

CREATE VIEW segments_bi AS 
SELECT distinct sp, ep
FROM (
 SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
 FROM (
   SELECT st_difference(b.seg, a.seg) as geom FROM 
      segments_diff_a as a, segments_diff_b as b 
      WHERE st_intersects(a.seg, b.seg)
    ) as cut
  ) as segments
  WHERE sp IS NOT NULL AND ep IS NOT NULL
ORDER BY sp, ep;

Ale z dziwnymi wynikami ...

wycięta wersja

huckfinn
źródło
Dziękuję za Twój wysiłek. Cóż, wyniki są dziwne. Zastanawiam się tylko, czy ST_HausdorffDistance () może pomóc odpowiedzieć na pytanie: gis.stackexchange.com/questions/180593/…
Lunar Sea
Hm, st_haudorffdistance daje miarę podobieństwa, a nie pożądanych segmentów (czerwone strzałki skierowane w stronę).
huckfinn
To tylko pomysł, ST_HausdorffDistance można wykorzystać do porównania geometrii obu tabel. Gdyby wielokąty nie były przestrzennie równe, stworzę linie. Po prostu nie wiem jak to zrobić.
Morze Księżycowe
Wydaje się, że jest to kwestia precyzji i topologii ( gis.stackexchange.com/a/182838/26213 i webhelp.esri.com/arcgisdesktop/9.2/… )
huckfinn
1

Patrząc na przykład, zmiana sugeruje, że funkcje z nowej tabeli, które zostały zmienione, zawsze będą się nakładać na funkcje ze starej tabeli. Dlatego skończyłbyś

ST_Overlaps (geoma, geomb) AND! ST_Touches (geoma, geomb)

Negacja dotyku polega na tym, że obiekty nakładają się również, jeśli tylko ich granice dzielą te same położenia wierzchołków.

przechylenie
źródło