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

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.

[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.

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.

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.