# How to test if geometry is inside another?

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
``````
-

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
``````
-
Thx for answer, but instead of using a for loop u should have used GetNextFeature(), otherwise ur SpatialFilter is useless – eouti May 13 '13 at 9:11

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.

-
I think the OP wants a solution in python and ogr/gdal. The tags don't include ArcPy. – Curlew Apr 26 '13 at 13:59

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

-
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... – eouti Apr 26 '13 at 11:15
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 . – Curlew Apr 26 '13 at 11:37
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... – eouti Apr 26 '13 at 12:26