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

Answer 2

You could add a field to your polygons and fill with some random but traceable value. Then do a Spatial join with the points, and all the points residing within a specific polygon will get the corresponding value.
Not sure if this will be any faster, but at least it's an alternative approach.

Answer 3

What you need is a Point-In-Polygon test, which is a standard algorithm and has been covered extensively in other sources.

See here or here or even this post here on GIS-SE