Punkt (z linii) w obrębie wielokąta za pomocą ogr i Python

Obecnie pracuję nad projektem, w którym muszę zbudować sieć topologiczną z elementów geometrii znajdujących się w plikach kształtów. Do tej pory, korzystając z projektu open source Bena Reilly, udało mi się przekształcić znaczniki linii w krawędzie networkx, a także wykryć bliskie cechy (mówią inne linie) i dodać je do najbliższego punktu, aby móc uruchomić algorytmy najkrótszej ścieżki.

Ale to dobrze dla jednego pliku kształtu. Jednak teraz muszę połączyć funkcje z różnych plików kształtów w duży wykres Networkx. Na przykład, jeśli punkt znajduje się w obrębie wielokąta, podłączę go (przez połączenie to znaczy dodaj krawędź networkx - add_edge (g.GetPoint (1), g.GetPoint (2)) z punktem w następnym pliku kształtu, który również znajduje się w wielokącie, który ma podobny atrybut (powiedzmy, ID). Zauważ, że wielokąty w różnych shps mają tylko te same identyfikatory, a nie współrzędne. Punkty, które mieszczą się w wielokątach, również nie mają tych samych współrzędnych.

Moim rozwiązaniem tego problemu było zidentyfikowanie punktu znajdującego się w wielokącie, przechowanie go, znalezienie punktu w następnym pliku kształtu znajdującym się w wielokącie o tym samym identyfikatorze, a następnie dodanie krawędzi networkx między nimi.

Jak sprawdzić, czy punkt znajduje się w wielokącie? Istnieje dobrze znany algorytm: algorytm RayCasting , który to robi. W tym momencie jednak utknąłem, ponieważ do implementacji algorytmu potrzebuję współrzędnych wielokąta i nie wiem, jak uzyskać do nich dostęp teraz, nawet po przejrzeniu dokumentacji Geometrii OGR. Pytanie, które zadaję, brzmi: jak uzyskać dostęp do punktów wielokąta lub współrzędnych LUB czy istnieje łatwiejszy sposób wykrycia, czy punkt mieści się w wielokącie? Używając Pythona z biblioteką osgeo.ogr kodowałem:

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

zobacz obraz, aby lepiej zrozumieć mój problem.

[EDYCJA] Do tej pory próbowałem przechowywać wszystkie obiekty wielokątów na liście, z którymi porównywałbym najpierw i na końcu linijkę . Ale przykład Paolo dotyczy używania odniesienia do obiektu punktu i odniesienia do obiektu wielokąta, który nie działałby z odniesieniem do obiektu linii, ponieważ nie cała linia znajduje się w obrębie wielokąta, a raczej pierwszy lub ostatni punkt jego linii.

[EDIT3] Utworzenie nowego obiektu punktu geometrii ze współrzędnych pierwszego i ostatniego punktu linii, a następnie użycie go do porównania z obiektami geometrii wielokąta zapisanymi na liście wydaje się działać dobrze:

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)

Dziękuję Paolo i Chrisowi za podpowiedzi.

Shapely jest fajny i elegancki, ale dlaczego nie używać jeszcze ogr, z operatorami przestrzennymi (w klasie OGRGeometry)?

przykładowy kod:

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)

Pamiętaj, że musisz skompilować GDAL z obsługą GEOS.

Nie znam siecix, ale jeśli dobrze zrozumiałem twoje pytanie, możesz użyć foremnej i OGR lib do znalezienia punktu w wielokącie z pliku kształtu . Oto jeden przykład, w jaki sposób działa sprawdzanie, czy jeden punkt (2000,1200) zawiedzie w przypadku dowolnego wielokąta z jednego pliku kształtu. W rezultacie drukuje współrzędne tego wielokąta.

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()

Mam nadzieję, że to pomoże.