Ponto (de uma cadeia de linhas) no Polygon usando ogr e Python

Atualmente, estou trabalhando em um projeto no qual preciso construir uma rede topológica a partir dos recursos de geometria encontrados nos shapefiles. Até agora, usando o projeto de código aberto de Ben Reilly, eu consegui transformar cadeias de linhas em bordas da redex, além de detectar recursos próximos (dizem outras cadeias de linhas) e adicioná-los ao ponto mais próximo para que eu possa executar algoritmos de caminho mais curto.

Mas isso é bom para um shapefile. No entanto, agora preciso conectar recursos de diferentes shapefiles a um grande gráfico de redex. Assim, por exemplo, se um ponto estiver dentro de um polígono, eu o conectaria (por conectá-lo, quero dizer adicionar uma borda de redex - add_edge (g.GetPoint (1), g.GetPoint (2)) com o ponto no próximo shapefile também está dentro de um polígono que compartilha um atributo semelhante (digamos, ID). Observe que os polígonos nos diferentes shps compartilham apenas os mesmos IDs e não as coordenadas.Os pontos que se enquadram nos polígonos também não compartilham as mesmas coordenadas.

Minha solução para esse problema foi identificar o ponto que reside em um polígono, armazená-lo, encontrar o ponto no próximo shapefile que reside no polígono com o mesmo ID e adicionar a borda da redex entre eles.

Como descobrir se um ponto reside dentro de um polígono? Bem, existe um algoritmo bem conhecido: o algoritmo RayCasting que faz isso. Foi aqui que fiquei realmente preso, porque, para implementar o algoritmo, preciso das coordenadas do polígono e não sei como acessá-las agora, mesmo depois de percorrer uma documentação da Geometria do OGR. Então, a pergunta que estou perguntando é como acessar os pontos poligonais ou coordenadas OU existe uma maneira mais fácil de detectar se um ponto se enquadra em um polígono? Usando python com a biblioteca osgeo.ogr, codifiquei o seguinte:

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

veja a imagem para uma melhor compreensão do meu problema.

[EDIT] Até agora, tentei armazenar todos os objetos de polígono em uma lista com a qual compararia o primeiro e o último ponto da cadeia de linhas . Mas o exemplo de Paolo está relacionado ao uso da referência de objeto de ponto e de referência de objeto de polígono, que não funcionaria com uma referência de objeto de linha, pois nem toda a linha está dentro do polígono, mas o primeiro ou o último ponto de sua cadeia de linhas.

[EDIT3] Criar um novo objeto de ponto de geometria a partir das coordenadas do primeiro e último ponto da cadeia de linhas e, em seguida, usá-lo para comparar com os objetos de geometria de polígono salvos em uma lista, parece funcionar bem:

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)

Obrigado Paolo e Chris pelas dicas.

Shapely é legal e elegante, mas por que não usar o ogr imóvel, com seus operadores espaciais (na classe OGRGeometry)?

Código de amostra:

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)

Observe que você precisa compilar o GDAL com o suporte do GEOS.

Eu não estou familiarizado com o networkx, mas se entendi corretamente sua pergunta, você pode usar a forma shapely e a OGR lib para encontrar pontos no polígono do shapefile. Aqui está um exemplo de como ele funciona para descobrir se um ponto (2000, 1200) falha com qualquer polígono de um shapefile. Para o resultado, ele imprime coordenadas desse polígono.

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

Espero que ajude.