Identyfikacja „długich i wąskich” wielokątów w PostGIS

10

Mam zestaw wielokątów reprezentujących duże obszary, powiedzmy dzielnice miast. Chcę zidentyfikować duże nakładające się obszary między nimi.

Ale jest problem: czasami wielokąty zachodzą na siebie na obwodzie (ponieważ zostały narysowane z małą precyzją). Spowoduje to generowanie długich i wąskich nakładek, na których mnie nie obchodzi.

Ale innym razem będą istnieć duże nakładki solidnych wielokątów, co oznacza duże obszary, w których wielokąt sąsiedztwa nakłada się na inny. Chcę wybrać tylko te.

Zobacz zdjęcie poniżej nakładających się na siebie. Wyobraź sobie, że chciałem wybrać tylko niebieski wielokąt w lewym dolnym rogu.

zakładka

Mogłem patrzeć na obszary, ale czasami wąskie są tak długie, że ostatecznie mają obszary tak duże jak niebieski wielokąt. Próbowałem zrobić stosunek powierzchni do obwodu, ale to również dało mieszane wyniki.

Próbowałem nawet użyć ST_MinimumClearance, ale czasami duże obszary będą miały przymocowaną wąską część lub dwa bardzo bliskie wierzchołki.

Jakieś pomysły na inne podejścia?


Ostatecznie najlepiej działało mi użycie bufora ujemnego, jak sugerują poniższe @Cyril i @FGreg.

Użyłem czegoś takiego jak:

ST_Area(ST_Buffer(geom, -10)) as neg_buffer_area

W moim przypadku jednostkami były metry, więc 10 m bufor ujemny.

W przypadku wąskich wielokątów obszar ten zwrócił zero (geometria byłaby pusta). Następnie użyłem tej kolumny do odfiltrowania wąskich wielokątów.

bplmp
źródło
4
Z pewnością można do tego wykorzystać stosunek powierzchni do obwodu.
Vince
Trudno powiedzieć, gdzie są wyraźne wielokąty z obrazu, ale zrobienie czegoś takiego, jak ten gis.stackexchange.com/a/265233/64838, może działać? Oblicz minimalną obróconą obwiednię, a następnie odrzuć te o małej szerokości lub wysokości.
FGreg
Możesz także spróbować użyć bufora ujemnego, jak opisano tutaj: Jak mogę zidentyfikować naprawdę cienkie wielokąty w moim pliku kształtu?
FGreg

Odpowiedzi:

5

Spróbowałbym utworzyć bufor ujemny, jeśli zjada cienkie wielokąty, to dobrze, jeśli nie je wielokąta, to jest mój ... :-)

uruchom ten skrypt, wcześniej ustawiając 2/3 szerokości wielokątów liniowych ...

create table name_table as
SELECT ST_Buffer(
(ST_Dump(
(ST_Union(
ST_Buffer(
(geom),-0.0001))))).geom,
0.0001)) as geom from source_table

System operacyjny: -) ...

Cyryl Michalchenko
źródło
ostatecznie twoja sugestia jest dla mnie najlepsza. Skończyłem używać czegoś takiego ST_Area(ST_Buffer(geom, -10)), w moim przypadku -10 to -10 metrów. Jeśli cokolwiek zwróciło 0 z tego wyrażenia, mógłbym je odfiltrować.
bplmp
9

Zamiast obszaru / obwodu lepiej jest użyć obszaru podzielonego przez kwadrat obwodu (lub jego odwrotność).

Nazywa się to również „wskaźnikiem kształtu”. Kwadrat obwodu podzielonego przez obszar ma minimalną wartość 4 * Pi () (w przypadku dysku, który jest najbardziej kompaktową geometrią 2D), więc można go znormalizować za pomocą 4 * Pi () interpretacja (znormalizowane wartości bliskie 1 oznaczają wtedy, że masz bardzo zwarte obiekty, a kwadraty mają wartości około 1,27).

EDYCJA: Próg na tym obszarze byłby przydatny do usunięcia bardzo małych artefaktów, które mogłyby być zwarte. Wtedy wskaźnik kształtu pokazałby lepszy kontrast. EDYCJA: oprócz tej odpowiedzi użycie ST_Snap może pomóc w rozwiązaniu problemu przed jego wystąpieniem.

radouxju
źródło
Dzięki! Ale nie jestem pewien, jak ST_Snap może pomóc w tym przypadku ... Jeśli mam rację, sugerujesz coś takiego (o.overlap_perimeter^2 / o.overlap_area) / (4 * Pi()) as overlap_ratio? To ma dla mnie gorsze wyniki niż tylko obszar / obwód.
bplmp
Teraz używam o.overlap_perimeter / (4 * sqrt(o.overlap_area)) as overlap_ratiozgodnie z tym artykułem, ale wciąż gorsze wyniki (chociaż trudno jest oszacować, co mam na myśli przez gorsze) isprs-ann-photogramm-remote-sens-spatial-inf-sci.net/I-7/135/… , strona 183.
bplmp
2
Dziękuję za to, nigdy nie słyszałem o „indeksie kształtów”. Zawsze myślałem, że użycie minimalnego prostokąta ograniczającego jest najlepszym sposobem na odpowiedź na tego rodzaju pytanie. Znalazłem to, repository.asu.edu/attachments/111230/content/… , co jest interesujące.
John Powell,
@JohnPowell intersting paper, dzięki. Widzę, że to, co nazywam indeksem kształtu, nazywa się w indeksie kołowym. Mój problem z minimalnymi prostokątami ograniczającymi polega na tym, że nie działa z bardzo wklęsłymi obiektami (np. W kształcie litery U)
radouxju
@bplmp ST_Snap pomoże ci przyciągnąć wierzchołki „prawie” sąsiadujących wielokątów, aby nie nakładały się już na siebie. Na twoich figurach nie ma skali, ale twój artefakt wygląda jak linie, więc myślę, że możesz użyć wartości tolerancji, która wystarcza, aby uniknąć artefaktów, ale nie wpływa na duże wielokąty.
radouxju
5

Jedną z opcji byłoby użycie stosunku powierzchni wielokąta do najdłuższej linii, którą można narysować za pomocą jego kończyn. Identyfikacja długich wąskich wielokątów.

select * from polygons where ST_Length(ST_LongestLine(geom, geom)) < ST_Area(geom) * 4

Działa to całkiem dobrze w przypadku wielokątów płaskich. Możesz dopasować współczynnik (przez pomnożenie obszaru) do własnych potrzeb i projekcji.

HeikkiVesanto
źródło
1

Wygląda na to, że może to pasować do Twojego przypadku użycia: Wyeliminuj wybrane wielokąty

Łączy wybrane wielokąty warstwy wejściowej z pewnymi sąsiadującymi wielokątami, usuwając ich wspólną granicę. Przylegający wielokąt może być albo tym, który ma największy lub najmniejszy obszar, albo tym, który dzieli największą wspólną granicę z wielokątem, który ma zostać wyeliminowany.

Eliminacja jest zwykle stosowana do pozbycia się wielokątów płaskich, tj. Małych wielokątów, które są wynikiem procesów przecinania się wielokątów, w których granice danych wejściowych są podobne, ale nie identyczne.

Wygląda na to, że chcesz wypróbować opcję „Największa wspólna granica”.

FGreg
źródło
Zdaję sobie sprawę, że teraz prosiliście o rozwiązania postgis, a nie rozwiązania qgis. Przepraszam, nie sądzę, że postgis ma równoważną funkcję, ale pozostawię to potomności.
FGreg
0

Wygląda mi to na idealny przypadek użycia dla rozszerzenia topologii PostGIS . Parametr tolerancji topologii określa, jak daleko pozwalasz wierzchołkom przyciągać się do innych istniejących wielokątów, radzić sobie z niską precyzją danych źródłowych i czyścić je.

Krótko mówiąc, strategia jest następująca:

1. Włącz rozszerzenie topologii

CREATE EXTENSION postgis_topology;

2. Utwórz nową pustą topologię

SELECT topology.CreateTopology('neighborhoods_topo', 4326, 1e-7);

Trzecim parametrem jest tolerancja w jednostkach CRS; wybierz mądrze. Idealnie, potrzebujesz CRS, gdzie jednostką jest metry. Jeśli jednostka CRS nie jest licznikiem, tak jak w przypadku WGS 84 aka 4326, użyj ST_Transformponownie rzutować swoje wielokąty.

3. Dodaj kolumnę TopoGeometry do tabeli wielokątów

SELECT topology.AddTopoGeometryColumn('neighborhoods_topo', 'public', 'neighborhoods', 'topogeom', 'POLYGON');

To zwraca nowy layer_id. Zapisz to, będzie potrzebne później. Będzie to warstwa, 1jeśli zaczniesz od zera, i będzie rosła przy każdym nowym połączeniu.

4. Dodaj wszystkie wielokąty do topologii

UPDATE public.neighborhoods
SET topogeom = topology.toTopoGeom(geom, 'neighborhoods_topo', 1, 1e-7);

Duży zestaw danych może potrwać kilka godzin, bądź cierpliwy. 1jest zwrócony wcześniej identyfikator_warstwy.

5. Znajdź twarze pojawiające się w kilku dzielnicach

Znajdź wszystkie twarze z topologii, które są obecne w 2 lub więcej topogeometrii. Zostawię zapytanie jako ćwiczenie. Najłatwiej jest prawdopodobnie z GetTopoGeomElementsfunkcją, a następnie pogrupuj według identyfikatora twarzy i spójrz na te z liczbą 2 lub więcej. Alternatywnie, możesz utworzyć nową tabelę z oczyszczoną geometrią z kolumny topogeom, po prostu rzutuj ją na standardową geometrię topogeom::geometryi powtórz to, co już masz, ale teraz z czystym zestawem danych bez nakładania się taśmy.

François
źródło