Verwenden von Shapely: Übersetzen zwischen Polygonen und MultiPolygonen

12

[EDIT: Die Lösung hierfür bestand einfach darin, OGR zum Lesen von Shapefiles zu verwenden. Siehe das Beispiel von geographika.]

In einem ESRI-Shapefile gibt es keinen Unterschied zwischen Polygonen und MultiPolygonen. Darüber hinaus gibt es keine explizite Unterscheidung zwischen Innenlöchern und Außenringen (abgesehen von der "Händigkeit" eines bestimmten Polygons).

Nachdem ich ein Shapefile gelesen habe, habe ich eine Liste von Koordinatensequenzen, die Ringe beschreiben, aber ohne eine intensivere Verarbeitung kann ich nicht unterscheiden, welche dieser Ringe Außenringe, Innenlöcher oder zusätzliche Polygone sind.

Es scheint , dass für wohlgeformte ‚s Polygon und Multipolygon Bauer, es muss eine klare Unterscheidung zwischen Außen- und Innenringen sein, so wie soll ich bewege sich von einer unklaren Liste der Ringe zu einem geordneten Satz von getrennten Polygone mit klar bezeichneten Innen- und Außenringe ?

Zusammenfassend: Wenn ich eine Liste von Polygonringen habe, aber nicht weiß, welche Ringe Löcher im Inneren oder separate Polygone sind, wie sollte ich sie am besten in separate Polygone mit bestimmten inneren Löchern sortieren?

Ich suche nach einer einfachen algorithmischen Lösung, die ich in Python implementieren kann, mit der Hunderte von Polygonen in einer Minute oder weniger verarbeitet werden können, und ich mache dies, um eine große Anzahl von Schnittpunkten auszuführen.

BenjaminGolder
quelle
In dieser Frage fehlen die entscheidenden Informationen darüber, was Sie zum Lesen des Shapefiles verwenden.
Inc42
@ inc42 Ich habe Python verwendet, um die Datei direkt zu lesen.
BenjaminGolder
Ah, dann ist das Shapely-Bit irreführend. Ihr eigentliches Problem bestand darin, herauszufinden, wie Sie die "Art" des Rings im Shapefile-Format bestimmen können. :)
inc42

Antworten:

10

Um die Antwort von relet zu erhalten, wie einzelne Polygone erhalten werden, können Sie dann einen Schnittpunkt für alle Polygone ausführen, um die Löcher zu erstellen. Wenn Ihr Datensatz überlappende Polygone enthält, haben Sie kein Glück.

Erklären Sie noch einmal, was mit vorhandenen Shapefile-Lesern nicht stimmt.

Wäre es nicht einfacher, Feature-IDs und M-Werte aus dem Shapefile zu exportieren und sie dann nach Verwendung eines vorhandenen Shapefile-Readers wieder mit den Polygonen zu verbinden?

Für Multipatches können Sie dieselbe Technik verwenden, um einer "Patch-ID" Polygon-IDs zuzuweisen und dieses Attribut dann wieder zu den Features hinzuzufügen.

Bearbeiten: Während Sie sagen, dass Sie OGR nicht verwenden möchten, nur für den Fall, dass Sie Ihre Meinung ändern ..

import ogr
# Get the driver
driver = ogr.GetDriverByName('ESRI Shapefile')
# Open a shapefile
shapefileName = "D:/temp/myshapefile.shp"
dataset = driver.Open(shapefileName, 0)

layer = dataset.GetLayer()
for index in xrange(layer.GetFeatureCount()):
    feature = layer.GetFeature(index)
    geometry = feature.GetGeometryRef()
    #geometry for polygon as WKT, inner rings, outer rings etc. 
    print geometry

Die Geometrie sollte wie folgt ausgegeben werden:

POLYGON ((79285 57742,78741 54273...),(76087 55694,78511 55088,..))

Die erste Klammer enthält die Koordinaten des Außenrings, die nachfolgenden Klammern die Koordinaten der Innenringe. Wenn Sie Z-Werte haben, sollten Punkte das Format 79285 57742 10 haben (wobei die letzte Koordinate eine Höhe ist).

Andernfalls können Sie die Funktionen Shapely Contains und Within verwenden, um jedes Polygon miteinander zu bewerten und zuvor einen räumlichen Index anzuwenden - http://pypi.python.org/pypi/Rtree/ , um die Verarbeitung zu beschleunigen.

geographika
quelle
Danke, das ist nicht das, wonach ich suche, aber ich habe einige weitere Klarstellungen gesehen, die ich zu meiner Frage machen könnte. Um Ihre Fragen zu beantworten: 1, weil ich sie sortiere, um überhaupt Tonnen von Kreuzungen zu machen. 2, nichts, ich brauche nur eine, die leicht ist und die ich leicht mit Python verwenden kann, und ich habe ogr nicht gelernt. 3, nein, in diesem Fall wäre es nicht einfacher.
BenjaminGolder
1
Beeindruckend! Danke geographika! Ich bin jetzt zu 90% davon überzeugt, dass ogr hier die Lösung ist. Ich habe einige Installationsprobleme mit ogr festgestellt, aber es sieht so aus, als ob es sich lohnt, sie durchzuarbeiten. Also organisiert ogr die Ringe in seiner wkt-Ausgabe? und ist es gut mit 3D-Shapefiles (ich benutze viel 3D)?
BenjaminGolder
Es sieht so aus, als ob die Antworten auf meine Fragen ja und ja sind. Und danke, dass Sie mich bei rtree vorgestellt haben. OGR scheint mein Problem vollständig zu lösen - solange ich die Installation dieses Mal richtig konfigurieren kann.
BenjaminGolder
Neuinstallation von OGR - Denken Sie daran, dass Sie die Binärdateien für Windows und die entsprechenden Python-Bindungen benötigen.
Geographika
9

Verwenden Sie zuerst ogr, um das Shapefile zu öffnen:

from osgeo import ogr
source = ogr.Open("mpolys.shp")
layers =  source.GetLayerByName("mpoly")
len(layers)
1

Shapefile-Geometrien in formschöne Geometrien umwandeln

from shapely.wkb import loads
element=layers[0] #(because lenght of layer =1, else you need "for element in layers: ...")
geom = loads(element.GetGeometryRef().ExportToWkb())
geom.geom_type
'MultiPolygon'
print geom
MULTIPOLYGON ((..... # the geometry in shapely wkt format

Für die Polygone im Multipolygon:

poly=[]
for pol in geom:
    poly.append(pol)
poly[0]
<shapely.geometry.polygon.Polygon object at 0x00B82CB0>
poly[0].geom_type
'Polygon'
print(poly[1])
POLYGON ((.... # the geometry in shapely wkt format

Und jetzt können Sie alle Funktionen von Shapely ( Shapely ) nutzen.

Constantinius
quelle
1

Ich bin nicht allzu vertraut damit, wie Polygone tatsächlich in Formdateien gespeichert werden, aber - sollte ein Polygonring nicht genau dann eine geschlossene Schleife sein, wenn die Startkoordinate wiederholt wird? Wenn Sie also jede nachfolgende Koordinate mit der Startkoordinate vergleichen, finden Sie den ersten Punkt, an dem ein Polygon geschlossen wird. Wenn dies die letzte Koordinate des Polygons ist, handelt es sich um ein einfaches Polygon. Wenn nicht, handelt es sich um ein Multipolygon, für das die anderen Schleifen verarbeitet werden müssen.

Das ist vielleicht die "intensivere Verarbeitung", die Sie vermeiden möchten, aber es ist wirklich nur eine Iteration durch die Koordinaten, die kostenlos ist, wenn Sie sie trotzdem einlesen müssen.

relet
quelle
Ich entschuldige mich - meine Frage war etwas irreführend und ich habe sie bearbeitet. Ich habe die Polygonringe und weiß, wann ich eine Liste mit mehr als einem Polygonring habe, aber abgesehen von der Reihenfolge der Punkte im oder gegen den Uhrzeigersinn weiß ich nicht, ob es sich um Innen- oder Außenringe handelt. Wenn es sich um Innenringe handelt, kann ich nicht sicher sein, zu welchem ​​Außenring sie gehören, abgesehen von der Messung ihrer Position.
BenjaminGolder
Aha. Ich bin auch froh zu sehen, dass du deinen Weg zu ogr gefunden hast. ;)
Relet
-2

Wie durch @ inc42 angegeben , finden Sie auf Seite 8 der Technischen Beschreibung des ESRI-Shapefiles: Ein ESRI-Weißbuch - Juli 1998 (Seite 12 von 34 im PDF) eine Diskussion über den Inhalt von Polygondatensätzen , nach denen Sie möglicherweise suchen :

Ein Polygon kann mehrere äußere Ringe enthalten. Die Reihenfolge der Eckpunkte oder die Ausrichtung eines Rings gibt an, auf welcher Seite des Rings sich das Innere des Polygons befindet.

PolyGeo
quelle
S.10 hat es "Ein Polygon kann mehrere äußere Ringe enthalten. Die Reihenfolge der Eckpunkte oder die Ausrichtung für einen Ring gibt an, welche Seite des Rings das Innere des Polygons ist."
Inc42
@ inc42 entsprechend aktualisiert, obwohl ich denke, dass die Seitenzahl 8 (oder 12 von 34) sein sollte.
PolyGeo