Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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
share|improve this question
add comment

3 Answers

up vote 2 down vote accepted

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
share|improve this answer
    
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
add comment

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.

share|improve this answer
    
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
add comment

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

share|improve this answer
    
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
1  
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
1  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.