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?

Thanks

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