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.

• Curlew: I think the OP wants a solution in python and ogr/gdal. The tags don't include ArcPy.

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

• eouti: I think the most efficient algo for my case is situated in the Within() method. I just need to use it efficiently, because I might be doing useless tests...
• Curlew: Ahh, i must have missed that you already used the Within() method. What do you mean with efficiency ? It is too slow? How many points do you have? Maybe you could preprocess your point shape before (identical points, exclude all points not within the extent of the polygon). Or maybe you can make use of threading .
• eouti: Yes it's very slow. I have around 15k polygons and 20k points. I'm looking for clever tricks to speed it up ;) Because right now, for a point, i apply Within() method n times. With n number of polygons...