Podziel złożony plik kształtu na siatkę

Mam dość szczegółowy plik kształtu z funkcjami wielokąta / wielokąta (plik ma około 500 MB). W rzeczywistości jest to plik kształtu całego świata, z funkcjami reprezentującymi linie brzegowe. Muszę podzielić te dane za pomocą siatki. Żeby było jasne, nie chcę „sortować” danych, ale przecinam wielokąty na kafelki. Zdaję sobie sprawę, że to pytanie zostało już zadane, ale znalezione przeze mnie rozwiązania nie zadziałały.

Próbowałem:

• Używanie QGIS i przecinanie zawartości mojego pliku kształtu z siatką wektorową - wyniki są okropne. Większość dużego lądu magicznie znika, choć wydaje się, że czasami robią to mniejsze kawałki ziemi. Powinienem zauważyć, że ta metoda działa naprawdę dobrze z dużo prostszymi danymi (tj. Mniej punktów)

• Korzystanie z narzędzi Przecięcia OGR. Wypróbowałem to zarówno przez ogr2ogr, jak i nawet uruchamiając własne narzędzie C ++. Oba mają ten sam problem co QGIS. Nie wykazują również tego problemu w przypadku prostych plików, ale zawodzą te bardziej złożone. Dla porównania używam pliku kształtu Australii i Nowej Zelandii, o wielkości poniżej 20 MB, a zarówno QGIS, jak i OGR nie potrafią go „ujednolicić”.

Ktoś zasugerował użycie PostGIS w pewnym momencie, ponieważ ma funkcję przecięcia - ale ST_Intersect PostGIS używa tego samego zaplecza GEOS, co OGR. W rzeczywistości oba wywołują tę samą funkcję, o ile mogę powiedzieć, więc nie sądzę, że PostGIS przyniesie różne wyniki.

Szukałem sugestii, co jeszcze mogę spróbować. Potrzebuję solidnej aplikacji lub zestawu narzędzi, który może dzielić bardzo szczegółowe pliki kształtów na kafelki.

EDYCJA: Dodanie dodatkowych informacji

W odpowiedzi na Simbamangu:

• Plik shapefile to w zasadzie dane wybrzeża z OpenStreetMap. Jest to scalona wersja pliku „przetwarzanego_p” (więc nie jest podzielony na kafelki), którą otrzymałem, wysyłając e-mail do ich listy deweloperów. Pamiętaj, że ich dzielenie płytek (na kawałki o wymiarach 100 km x 100 km z nakładaniem się) niekoniecznie jest tym, czego chcę - nie chcę nakładania się i chcę swobodę wyboru rozmiaru siatki, lub po prostu skorzystam z domyślnie przetworzone_p.

• Domyślnie dane linii brzegowej zawierają błędy geometrii zgłoszone przez QGIS. Naprawiam te błędy za pomocą małego narzędzia, które przygotowałem przy użyciu kodu, który został zaprojektowany specjalnie w celu rozwiązania tego problemu (naprawa błędów geometrii w danych linii brzegowej: https://github.com/tudelft-gist/prepair ). Uruchomienie plików za pomocą tego narzędzia naprawia praktycznie wszystkie błędy, które pobiera QGIS. Próbuję wykonać skrzyżowanie dopiero po wyczyszczeniu plików.

• Dokładnie to, co zrobiłem przy użyciu QGIS: Otwórz dane, aby upewnić się, że w QGIS wygląda dobrze. Spróbuj podzielić go na kafelki, tworząc warstwę kafelków za pomocą siatki wektorowej o określonych odstępach, a następnie przecinając dwie warstwy - nie ma mowy. Spróbuj użyć mniejszego zestawu danych - wybierz funkcje w Oceanii (Aus, NZ), aby wypróbować mniejszy zestaw danych - ten plik kształtu ma rozmiar <20 MB. Ponownie spróbuj go podzielić, nie działa.

• Co zrobiłem z OGR: ogr2ogr bezpośrednio przy użyciu opcji „-spat” i „-clipsrc” z spat_extent. Napisałem również małe narzędzie C ++, które działa na WKT, więc przekonwertowałem plik shapefile na WKT za pomocą ogr2ogr, a następnie podaję plik tekstowy do mojej aplikacji. Przebiega przez plik i wywołuje udokumentowaną tutaj metodę Intersection (): http://www.gdal.org/ogr/classOGRGeometry.html . Myślę, że ostatecznie robi to dokładnie to samo, co bezpośrednie używanie ogr2ogr.

W odpowiedzi na Brent:

1. To robi. Wszystko jest w WGS84 Lat / Lon
2. Wydawało mi się, że jest odwrotnie - że dla danego zestawu kafelków siatki potrzeba więcej czasu, by przeciąć jeden gigantyczny wielokąt, niż kilka fragmentów fragmentów, które mogłyby być bardziej przestrzennie zlokalizowane na każdym kafelku, ale to jest ciekawa propozycja - spróbuję i zdam raport.
3. Podczas procesu nie są przechowywane pola atrybutów, interesuje mnie tylko geometria.
4. Nie jestem pewien, ale myślę, że mówisz, że powinienem wybrać wielokąty, które nachodzą na dany kafelek siatki, a następnie wykonać skrzyżowanie. Jest to zbyt uciążliwe ręcznie w QGIS. Moje narzędzie już to robi w pewnym stopniu, zaznaczając obwiednię. Jest trochę przyspieszenia, ale wynik końcowy jest nadal słaby i nie jest zauważalnie inny.
5. To nie jest opcja. W tej chwili próbuję podzielić dane tak, aby były one 1 stopień lat x 1 stopień l i szukam ogólnej / solidnej metodologii, która działałaby we wszystkich przypadkach. Próbowałem zwiększyć rozmiar siatki (tj. 10x10), aby zobaczyć, czy uzyskałem lepsze wyniki i nie widzę żadnej korelacji między rozmiarem siatki a jakością wydruku.

Edytuj # 2:

Próbowałem grać z tym więcej i ogólnie wydaje się, że wyniki są niewiarygodne zarówno przy użyciu GEOS, jak i QGIS (który używa fTools, nie wiem, czy to z kolei ponownie używa GEOS). Myliłem się, twierdząc, że rozmiar siatki nie ma nic wspólnego z wynikami - im większa siatka, tym lepsze wyniki (dobrze wiedzieć, ale nadal nie jest to rozwiązanie). Oto zrzut ekranu z naprawdę rozłożoną siatką, która w większości działała, ale częściowo zawiodła w jednym kafelku:

Geometria jest czysta - QGIS pokazuje 0 błędów za pomocą narzędzia „Sprawdź ważność”. Nie zamierzam podchodzić do tego problemu krok po kroku; sprawdzanie, czy pewne funkcje nie udały się przeciąć tak dużego zestawu danych, gdy nie jest to widoczne (i nie będzie z mniejszymi kafelkami), nie jest praktyczne.

Właśnie skończyłem tworzyć własne narzędzia do tego.

Korzystałem z biblioteki Clippera ( http://www.angusj.com/delphi/clipper.php ) wraz z OGR, aby podzielić konfigurację danych. Warto zauważyć, że naiwne wykonywanie skrzyżowań z tą biblioteką zajmuje bardzo dużo czasu, więc zamiast tego zastosowałem podejście z poczwórnym drzewem ... tj. Podzielę na cztery komórki siatki, podziel każdą z nich na jeszcze cztery itd., Aż uzyskasz pożądaną rozdzielczość. Biblioteka działa świetnie, załączyłem zrzut ekranu pokazujący wyniki na wschodniej półkuli:

Powyższy wynik zajął około 4,5 godziny na procesorze 1,33 GHz.

Oto narzędzia na wypadek, gdyby ktoś napotkał podobny problem w przyszłości. Należy pamiętać, że zostały one zhakowane razem proof-of-concept i prawdopodobnie nie należy ich używać bezpośrednio (może to być dobry punkt wyjścia do czegoś):

https://github.com/preet/scratch/tree/master/gis/polytoolkit

https://github.com/preet/scratch/tree/master/gis/shapefiles/shptk

To na pewno brzmi, jakbyś miał problemy z geometrią. Jest mało prawdopodobne, że będzie w stanie uzyskać czyste wyniki z brudnego pliku wejściowego niezależnie od używanego oprogramowania, chyba że najpierw rozwiążesz problemy z geometrią. Po rozwiązaniu problemów z geometrią możesz spróbować wykonać następujące czynności, jeśli nadal występują problemy:

1) Upewnij się, że zestaw danych siatki ma taką samą projekcję jak zestaw danych wielokąta świata. Jeśli nie, odtwórz go we właściwej projekcji.

2) Konwertuj wszystkie funkcje na jedną część - znacznie łatwiejsze do przetworzenia

3) Usuń wszystkie obce pola, zachowując tylko pole identyfikatora, co pozwoli ci połączyć atrybuty z powrotem po wykonaniu skrzyżowania - ponownie znacznie łatwiejsze do przetworzenia

4) Zamiast przecinać cały zestaw danych siatki z całym zestawem danych wielokąta świata, spróbuj zapętlić wielokąty siatki, wybierając przecinające się wielokąty w swoim zestawie danych świata i wykonując klip na podstawie wielokąta siatki. Umożliwi to izolowanie wszelkich problemów, a na koniec możesz połączyć wyniki razem, aby osiągnąć swój pierwotny cel.

5) Spróbuj użyć większych wielokątów siatki.

Innym podejściem może być próba konwersji wektora na raster w celu utworzenia zestawu danych punktów, a następnie wykorzystanie zestawu danych punktów jako podstawy do napisania kodu do utworzenia kafelków.