Punto (de una cadena lineal) dentro de Polygon usando ogr y Python

10

Actualmente estoy trabajando en un proyecto en el que necesito construir una red topológica a partir de las características de geometría que encuentro en los archivos de forma. Hasta ahora, utilizando el proyecto de código abierto de Ben Reilly, he logrado transformar cadenas lineales en bordes networkx, así como detectar características cercanas (dicen otras cadenas lineales) y agregarlas al punto más cercano para que pueda ejecutar algoritmos de ruta más corta.

Pero eso está bien para un archivo de formas. Sin embargo, ahora necesito conectar características de diferentes archivos de forma en un gran gráfico networkx. Entonces, por ejemplo, si un punto está dentro de un polígono, lo conectaría (al conectarlo me refiero a agregar un borde networkx - add_edge (g.GetPoint (1), g.GetPoint (2)) con el punto en el siguiente archivo de forma que también está dentro de un polígono que comparte un atributo similar (por ejemplo, ID). Tenga en cuenta que los polígonos en los diferentes shps solo comparten las mismas ID y no las coordenadas. Los puntos que se encuentran dentro de los polígonos tampoco comparten las mismas coordenadas.

Mi solución a este problema fue identificar el punto que reside en un polígono, almacenarlo, encontrar el punto en el siguiente archivo de forma que reside en el polígono con la misma identificación y luego agregar el borde de networkx entre ellos.

¿Cómo encontrar si un punto reside dentro de un polígono? Bueno, hay un algoritmo bien conocido: el algoritmo RayCasting que hace eso. Sin embargo, aquí es donde realmente me quedé atascado, porque para implementar el algoritmo necesito las coordenadas del polígono, y no sé cómo acceder a ellas en este momento, incluso después de hojear una documentación de la Geometría de OGR. Entonces, la pregunta que hago es cómo acceder a los puntos del polígono, o las coordenadas O ¿hay una manera más fácil de detectar si un punto cae dentro de un polígono? Usando python con la biblioteca osgeo.ogr codifiqué lo siguiente:

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

Vea la imagen para una mejor comprensión de mi problema. texto alternativo

[EDITAR] Hasta ahora he intentado almacenar todos los objetos poligonales en una lista con la que luego compararía el primer y el último punto de la cadena lineal . Pero el ejemplo de Paolo está relacionado con el uso de la referencia de objeto de punto y la referencia de objeto de polígono, que no funcionaría con una referencia de objeto de línea ya que no toda la línea está dentro del polígono, sino más bien el primer o el último punto de su cadena lineal.

[EDITAR3] Crear un nuevo objeto de punto de Geometría a partir de las coordenadas del primer y último punto de la cadena de líneas y luego usarlo para compararlo con los objetos de geometría de polígono guardados en una lista parece funcionar bien:

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)

Gracias Paolo y Chris por las pistas.

python shapefile geometry ogr networkx usuario39901230
fuente

Respuestas:

11

Shapely es genial y elegante, pero ¿por qué no usar todavía ogr, con sus operadores espaciales (en la clase OGRGeometry)?

Código de muestra:

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)

Tenga en cuenta que necesita compilar GDAL con soporte GEOS.

capooti
fuente
Grzie Paolo, pero lo que realmente necesito no es comparar las referencias de objeto de un punto con la de un polígono, sino el último y el primer punto de una cadena lineal, a la que actualmente accedo con GetPoint. No hay GetPointRef que noté. ¿Cómo haría para implementar eso?
user39901230
1
La cadena lineal es una colección de geometría iterable de punto (vértice), puede acceder fácilmente al primero y al último de ellos.
capooti
Intenté almacenar los objetos poligonales en una lista que luego usaría para verificar el primer y último punto de una cadena lineal. Sin embargo, no se agregan puntos a la lista de puntos dentro de un polígono por alguna razón: eche un vistazo aquí: codepad.org/Cm2BV5mp
user39901230
8

No estoy familiarizado con networkx, pero si entendí correctamente su pregunta, podría usar lib bien shapely y OGR para encontrar puntos en el polígono desde shapefile. Aquí hay un ejemplo de cómo funciona para encontrar si un punto (2000,1200) falla dentro de cualquier polígono de un archivo de forma. Para el resultado, imprime las coordenadas de ese 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 ayude.

Mario Miler
fuente