[GIS] Testing if geometry is inside another using OGR in Python

Question

I need to tell whether geometries from a shapefile are inside other geometries or not.
For example, I want to test if there are trees (represented as points in trees.shp) inside the urban areas (reprensented as polygons in area.shp)

What's the easiest way to proceed?

I loop over the points, and within this loop, I have a nested loop over the polygons.
Then, I test if my point geometry is inside the polygon with the Within(Geometry geo) method from OGR.
It's very long, not optimized at all.

What should I do to fasten the computation?

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

You could try only doing the Within() test against points that are within the polygon features envelopes, i.e.

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