Punkt (eines Linienstreifens) innerhalb von Polygon mit ogr und Python

10

Ich arbeite derzeit an einem Projekt, in dem ich aus den Geometrie-Features, die ich in Shapefiles finde, ein topologisches Netzwerk erstellen muss. Bisher habe ich es mit dem Open-Source-Projekt von Ben Reilly geschafft, Linestrings in Networkx-Kanten umzuwandeln, enge Features zu erkennen (sagen andere Linestrings) und sie dem nächstgelegenen Punkt hinzuzufügen, damit ich Algorithmen für kürzeste Wege ausführen kann.

Aber das ist gut für ein Shapefile. Jetzt muss ich jedoch Features aus verschiedenen Shapefiles in einem großen networkx-Diagramm verbinden. Wenn sich ein Punkt beispielsweise innerhalb eines Polygons befindet, würde ich ihn verbinden (mit Verbinden meine ich das Hinzufügen einer networkx-Kante - add_edge (g.GetPoint (1), g.GetPoint (2)) mit dem Punkt im nächsten Shapefile befindet sich auch innerhalb eines Polygons, das ein ähnliches Attribut (z. B. ID) aufweist. Beachten Sie, dass die Polygone in den verschiedenen Shps nur dieselben IDs und nicht die Koordinaten verwenden. Die Punkte, die in die Polygone fallen, haben auch nicht dieselben Koordinaten.

Meine Lösung für dieses Problem bestand darin, den Punkt zu identifizieren, der sich in einem Polygon befindet, ihn zu speichern, den Punkt im nächsten Shapefile zu finden, der sich im Polygon mit derselben ID befindet, und dann die networkx-Kante zwischen ihnen hinzuzufügen.

Wie finde ich heraus, ob sich ein Punkt in einem Polygon befindet? Nun, es gibt einen bekannten Algorithmus: den RayCasting- Algorithmus, der das macht. Hier blieb ich jedoch stecken, denn um den Algorithmus zu implementieren, benötige ich die Koordinaten des Polygons und weiß nicht, wie ich jetzt darauf zugreifen soll, selbst nachdem ich eine Dokumentation der Geometrie des OGR durchgesehen habe . Die Frage, die ich stelle, ist also, wie auf die Polygonpunkte oder Koordinaten zugegriffen werden kann ODER ob es eine einfachere Möglichkeit gibt, festzustellen, ob ein Punkt in ein Polygon fällt. Unter Verwendung von Python mit der Bibliothek osgeo.ogr habe ich Folgendes codiert:

 if g.GetGeometryType() == 3: #polygon
                c = g.GetDimension()
                x = g.GetPointCount()
                y = g.GetY()
                z = g.GetZ()

Sehen Sie sich das Bild an, um mein Problem besser zu verstehen. Alt-Text

[BEARBEITEN] Bisher habe ich versucht, alle Polygonobjekte in einer Liste zu speichern, mit der ich dann den ersten und letzten Linestring vergleichen würde . Das Beispiel von Paolo bezieht sich jedoch auf die Verwendung der Punktobjektreferenz und der Polygonobjektreferenz, die mit einer Linienobjektreferenz nicht funktionieren würden, da sich nicht die gesamte Linie innerhalb des Polygons befindet, sondern der erste oder letzte Punkt seines Linienstreifens.

[EDIT3] Das Erstellen eines neuen Geometriepunktobjekts aus den Koordinaten des ersten und letzten Punkts des Linienstreifens und das anschließende Vergleichen mit den in einer Liste gespeicherten Polygongeometrieobjekten scheint einwandfrei zu funktionieren:

for findex in xrange(lyr.GetFeatureCount()):
    f = lyr.GetFeature(findex)
    flddata = getfieldinfo(lyr,f,fields)
    g = f.geometry()
    if g.GetGeometryType() == 2:
        for j in xrange(g.GetPointCount()):
            if j == 0 or j == g.GetPointCount():
                point = ogr.Geometry(ogr.wkbPoint)
                point.AddPoint(g.Getx(j),g.GetY(j))
                if point.Within(Network.polygons[x][0].GetGeometryRef()):
    print g.GetPoint(j)

Danke Paolo und Chris für die Hinweise.

user39901230
quelle

Antworten:

11

Shapely ist cool und elegant, aber warum nicht noch ogr mit seinen räumlichen Operatoren (in der OGRGeometry-Klasse) verwenden?

Beispielcode:

from osgeo import ogr
driver = ogr.GetDriverByName('ESRI Shapefile')
polyshp = driver.Open('/home/pcorti/data/shapefile/multipoly.shp')
polylyr = polyshp.GetLayer(0)
pointshp = driver.Open('/home/pcorti/data/shapefile/point.shp')
pointlyr = pointshp.GetLayer(0)
point = pointlyr.GetNextFeature()
polygon = polylyr.GetNextFeature()
point_geom = point.GetGeometryRef()
polygon_geom = polygon.GetGeometryRef()
print point_geom.Within(polygon_geom)

Beachten Sie, dass Sie GDAL mit GEOS-Unterstützung kompilieren müssen.

Capooti
quelle
grzie Paolo, aber was ich wirklich brauche, ist nicht, die Objektreferenzen eines Punktes mit denen eines Polygons zu vergleichen, sondern den letzten und den ersten Punkt eines Linestrings, auf den ich derzeit mit GetPoint zugreife. Es gibt kein GetPointRef, das mir aufgefallen ist. Wie würde ich das umsetzen?
user39901230
1
linestring ist eine iterierbare Geometriesammlung von Punkten (Scheitelpunkten), auf die Sie problemlos zugreifen können.
Capooti
Ich habe versucht, die Polygonobjekte in einer Liste zu speichern, die ich später verwenden würde, um sie mit dem ersten und letzten Punkt eines Linienstreifens zu vergleichen. Aus irgendeinem Grund werden der Liste der Punkte innerhalb eines Polygons jedoch keine Punkte hinzugefügt: Schauen Sie hier: codepad.org/Cm2BV5mp
user39901230
8

Ich bin nicht vertraut mit NetworkX aber wenn ich deine Frage richtig verstanden, könnten Sie shapely und OGR lib Punkt im Polygon von Shape - Datei für die Suche. Hier ist ein Beispiel, wie es funktioniert, um festzustellen, ob ein Punkt (2000, 1200) mit einem Polygon aus einem Shapefile fehlschlägt. Für das Ergebnis werden Koordinaten dieses Polygons gedruckt.

from shapely.geometry import Point, Polygon
from shapely.wkb import loads
from osgeo import ogr

file1 = ogr.Open("d:\\fileWithData.shp")
layer1 = file1.GetLayerByName("fileWithData")

point1 = Point(2000,1200)

polygon1 = layer1.GetNextFeature()

while polygon1 is not None:
    geomPolygon = loads(polygon1.GetGeometryRef().ExportToWkb())
    if geomPolygon.contains(point1):
        xList,yList = geomPolygon.exterior.xy
        print xList
        print yList
    polygon1 = layer1.GetNextFeature()

Ich hoffe, es hilft.

Mario Miler
quelle