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

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.

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;

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;

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)

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:

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:

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

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

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

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;

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;

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

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.