Teilen Sie ein komplexes Shapefile in ein Raster

Ich habe ein anständig detailliertes Shapefile mit Polygon- / Multipolygon-Features (die Datei ist ungefähr 500 MB groß). Es ist eigentlich ein Shapefile der ganzen Welt, dessen Merkmale die Küsten darstellen. Ich muss diese Daten durch ein Raster teilen. Um es klar zu sagen, ich möchte die Daten nicht "sortieren", sondern die Polygone in Kacheln schneiden. Mir ist klar, dass diese Frage schon einmal gestellt wurde, aber die Lösungen, die ich gefunden habe, haben bei mir nicht funktioniert.

Ich habe es versucht:

• Wenn Sie QGIS verwenden und meinen Shapefile-Inhalt mit einem Vektorgitter schneiden, sind die Ergebnisse schrecklich. Der größte Teil der großen Landmasse verschwindet auf magische Weise, obwohl es manchmal so aussieht, als ob kleinere Landstücke es schaffen. Ich sollte beachten, dass diese Methode mit viel einfacheren Daten (dh weniger Punkten) sehr gut funktioniert.

• Verwenden der OGR-Schnittwerkzeuge. Ich habe es sowohl mit ogr2ogr als auch mit meinem eigenen C ++ - Tool versucht. Beide haben das gleiche Problem wie QGIS. Sie weisen dieses Problem auch bei einfachen Dateien nicht auf, scheitern jedoch bei den komplexeren. Als Referenz verwende ich ein Shapefile aus Australien und Neuseeland mit einer Größe von weniger als 20 MB, und sowohl QGIS als auch OGR können es nicht "Gridifizieren".

Jemand schlug vor, PostGIS an einer Stelle zu verwenden, da es eine Schnittfunktion hat - aber STGIIntersect von PostGIS verwendet dasselbe GEOS-Backend wie OGR. Tatsächlich rufen beide, soweit ich das beurteilen kann, dieselbe Funktion auf, sodass ich nicht glaube, dass PostGIS zu unterschiedlichen Ergebnissen führen wird.

Ich suchte nach Vorschlägen, was ich sonst noch versuchen könnte. Ich benötige eine robuste Anwendung oder ein Toolkit, mit dem sehr detaillierte Shapefiles in Kacheln unterteilt werden können.

BEARBEITEN: Weitere Informationen hinzufügen

Als Antwort auf Simbamangu:

• Das Shapefile besteht im Wesentlichen aus Küstendaten von OpenStreetMap. Es ist eine zusammengeführte Version der 'verarbeiteten_p'-Datei (also nicht in Kacheln aufgeteilt), die ich per E-Mail an ihre Entwicklerliste erhalten habe. Beachten Sie, dass die Aufteilung der Kacheln (in 100 km x 100 km große Blöcke mit Überlappung) nicht unbedingt das ist, was ich möchte - ich möchte keine Überlappung und ich möchte die Freiheit, die Rastergröße zu wählen, oder ich würde einfach die verwenden Standard verarbeitet_p.

• Standardmäßig weisen die Küstendaten Geometriefehler auf, die von QGIS gemeldet werden. Ich behebe diese Fehler mit einem kleinen Tool, das ich mithilfe eines Codes zusammengestellt habe, der speziell für dieses Problem entwickelt wurde (Reparieren von Geometriefehlern in Küstendaten: https://github.com/tudelft-gist/prepair ). Das Überfahren der Dateien mit diesem Tool behebt praktisch alle Fehler, die QGIS auffängt. Ich versuche erst, die Kreuzung nach dem Bereinigen der Dateien durchzuführen.

• Genau das, was ich mit QGIS gemacht habe: Öffnen Sie die Daten, um sicherzustellen, dass sie in QGIS gut aussehen. Versuchen Sie, es in Kacheln zu unterteilen, indem Sie mit Vector Grid eine Ebene mit Kacheln mit einem bestimmten Abstand erstellen und dann die beiden Ebenen schneiden - no go. Versuchen Sie es mit einem kleineren Datensatz. Wählen Sie Features in Ozeanien (Aus, NZ) aus, um einen kleineren Datensatz zu verwenden. Diese Formdatei ist <20 MB groß. Versuchen Sie es erneut zu teilen, funktioniert nicht.

• Was ich mit OGR gemacht habe: ogr2ogr direkt mit den Optionen '-spat' und '-clipsrc' mit spat_extent. Außerdem wurde ein kleines C ++ - Tool geschrieben, das mit WKT funktioniert. Daher konvertiere ich das Shapefile mit ogr2ogr in WKT und füge die Textdatei meiner Anwendung zu. Es durchläuft die Datei und ruft die hier dokumentierte Intersection () -Methode auf: http://www.gdal.org/ogr/classOGRGeometry.html . Ich denke, es macht genau das Gleiche wie die direkte Verwendung von ogr2ogr.

Als Antwort auf Brent:

1. Es tut. Alles ist in WGS84 Lat / Lon
2. Ich hätte gedacht, dass das Gegenteil der Fall ist - dass es für einen bestimmten Satz von Gitterkacheln viel länger dauern würde, ein riesiges Multipolygon zu schneiden, als eine Reihe fragmentierter Merkmale, die räumlich auf jede Kachel lokalisiert werden könnten, aber das ist es Ein interessanter Vorschlag - ich werde es versuchen und zurückmelden.
3. Während des Prozesses werden keine Attributfelder gespeichert, ich interessiere mich nur für Geometrie.
4. Ich bin mir nicht sicher, aber ich denke, Sie sagen, ich sollte die Polygone auswählen, die eine bestimmte Gitterkachel überlappen, und dann den Schnittpunkt ausführen. Dies ist mit QGIS manuell zu umständlich. Mein Tool macht dies bis zu einem gewissen Grad bereits mit einem Begrenzungsrahmen-Check. Es gibt eine gewisse Beschleunigung, aber das Endergebnis ist immer noch schlecht und nicht merklich anders.
5. Dies ist keine Option. Im Moment versuche ich, die Daten so aufzuteilen, dass sie 1 ° lat x 1 ° lon sind, und ich suche nach einer allgemeinen / robusten Methode, die in allen Fällen funktioniert. Ich habe versucht, die Rastergröße (dh 10 x 10) zu erhöhen, um bessere Ergebnisse zu erzielen, und ich sehe keine Korrelation zwischen der Rastergröße und der Qualität der Ausgabe.

Edit # 2:

Ich habe versucht, mehr damit herumzuspielen, und im Allgemeinen scheint es nur so, als ob die Ergebnisse sowohl mit GEOS als auch mit QGIS (das fTools verwendet, ich weiß nicht, ob das wiederum GEOS verwendet) unzuverlässig sind. Ich habe mich geirrt, als ich angegeben habe, dass die Größe des Rasters nichts mit den Ergebnissen zu tun hat - je größer das Raster ist, desto besser sind die Ergebnisse (das ist gut zu wissen, aber immer noch keine Lösung). Hier ist ein Screenshot eines wirklich verteilten Gitters, das größtenteils funktioniert hat, aber teilweise in einer Kachel versagt hat:

Die Geometrie ist sauber - QGIS zeigt 0 Fehler mit dem Werkzeug "Gültigkeit prüfen" an. Ich möchte dieses Problem nicht Schritt für Schritt angehen. Es ist nicht praktikabel zu überprüfen, ob bestimmte Features den Schnittpunkt eines so großen Datensatzes nicht bestanden haben, wenn er nicht visuell sichtbar ist (und dies bei kleineren Kacheln nicht der Fall ist).

Ich habe gerade meine eigenen Tools dafür erstellt.

Ich habe die Clipper-Bibliothek ( http://www.angusj.com/delphi/clipper.php ) zusammen mit OGR verwendet, um meine Daten zu teilen. Zu beachten ist, dass das naive Durchführen von Schnittpunkten mit dieser Bibliothek sehr lange dauert. Daher habe ich stattdessen einen Quadtree-Ansatz verwendet ... dh in vier Gitterzellen teilen, jede in vier weitere teilen usw., bis Sie die gewünschte Auflösung erhalten. Die Bibliothek funktioniert jedoch hervorragend. Ich habe einen Screenshot angehängt, der die Ergebnisse auf der östlichen Hemisphäre zeigt:

Das obige Ergebnis dauerte auf einem 1,33-GHz-Prozessor etwa 4,5 Stunden.

Hier sind die Tools für den Fall, dass jemand in Zukunft auf ein ähnliches Problem stößt. Bitte beachten Sie, dass sie Proof-of-Concepts zusammen gehackt haben und Sie sie wahrscheinlich nicht direkt verwenden sollten (könnte jedoch als guter Ausgangspunkt für etwas dienen):

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

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

Es klingt definitiv so, als hätten Sie Geometrieprobleme. Es ist unwahrscheinlich, dass unabhängig von der verwendeten Software saubere Ergebnisse aus einer verschmutzten Eingabedatei erzielt werden können, es sei denn, Sie gehen zuerst auf Ihre Geometrieprobleme ein. Sobald Sie Ihre Geometrieprobleme gelöst haben, können Sie Folgendes versuchen, wenn Sie immer noch Probleme haben:

1) Stellen Sie sicher, dass Ihr Raster-Dataset dieselbe Projektion wie Ihr Welt-Polygon-Dataset hat. Wenn nicht, erstellen Sie es in der richtigen Projektion neu.

2) Konvertieren Sie alle Features in ein Teil - viel einfacher zu verarbeiten

3) Entfernen Sie alle überflüssigen Felder, und behalten Sie nur das ID-Feld bei, damit Sie Ihre Attribute nach dem Ausführen der Kreuzung wieder zusammenfügen können - wiederum viel einfacher zu verarbeiten

4) Anstatt das gesamte Raster-Dataset mit dem gesamten Welt-Polygon-Dataset zu schneiden, versuchen Sie, Ihre Raster-Polygone zu durchlaufen, die sich überschneidenden Polygone in Ihrem Welt-Dataset auszuwählen und einen Clip basierend auf Ihrem Raster-Polygon auszuführen. Auf diese Weise können Sie alle Probleme eingrenzen und am Ende die Ergebnisse zusammenführen, um Ihr ursprüngliches Ziel zu erreichen.

5) Versuchen Sie, größere Gitterpolygone zu verwenden.

Ein anderer Ansatz könnte darin bestanden haben, eine Vektor-Raster-Konvertierung zum Erstellen eines Punktdatensatzes zu versuchen und dann den Punktdatensatz als Grundlage für das Schreiben von Code zum Erstellen Ihrer Kacheln zu verwenden.