¿Probar si una geometría está dentro de otra usando OGR en Python?

Question

Necesito saber si las geometrías de un shapefile están dentro de otras geometrías o no. Por ejemplo, quiero comprobar si hay árboles (representados como puntos en trees.shp ) dentro de las zonas urbanas (representadas como polígonos en area.shp )

¿Cuál es la forma más fácil de proceder?

Hago un bucle sobre los puntos, y dentro de este bucle, tengo un bucle anidado sobre los polígonos. A continuación, compruebo si mi geometría de puntos está dentro del polígono con el Within(Geometry geo) de OGR. Es muy largo, no está optimizado en absoluto.

¿Qué debo hacer para acelerar el cómputo?

def PointsInPolygons(self, pointsLayer, polysLayer):
    nbInside = 0
    total = pointsLayer.GetFeatureCount()
    for i in range(total):
        pointFeat = pointsLayer.GetFeature(i)
        pointGeo = pointFeat.GetGeometryRef()

        for j in range(polysLayer.GetFeatureCount()):
            polyFeat = polysLayer.GetFeature(j)
            polyGeo = polyFeat.GetGeometryRef()

            if pointGeo.Within(polyGeo):
                nbInside+=1
                break

Answer 1

Podrías intentar hacer sólo el Within() prueba contra los puntos que están dentro de las características del polígono sobres es decir

def PointsInPolygons(pointsLayer, polysLayer):
    nbInside = 0
    polyFeat=polysLayer.GetNextFeature()
    while polyFeat:
        polyGeo = polyFeat.GetGeometryRef()
        pointsLayer.SetSpatialFilter(polyGeo) #<----Only test points within feature envelope

        pointFeat = pointsLayer.GetNextFeature()
        while pointFeat:
            pointGeo = pointFeat.GetGeometryRef()

            if pointGeo.Intersects(polyGeo):
                nbInside+=1

            pointFeat = pointsLayer.GetNextFeature()

        polyFeat=polysLayer.GetNextFeature()
    return nbInside

Answer 2

Podrías añadir un campo a tus polígonos y rellenarlo con algún valor aleatorio pero rastreable. Luego hacer un Unión espacial con los puntos, y todos los puntos que residan dentro de un polígono específico obtendrán el valor correspondiente.
No estoy seguro de que esto sea más rápido, pero al menos es un enfoque alternativo.

Answer 3

Lo que necesitas es una prueba de Punto en Polígono, que es un algoritmo estándar y ha sido cubierto ampliamente en otras fuentes.

Ver aquí o aquí o incluso esto publicar aquí en GIS-SE