How to find out which polygon contains the most points using Python OGR?

Question

I have two shapefiles: point (railway stations) and polygon (regions). My task is to write a script using OGR module in Python that can determine in which region the biggest amount of railway stations is.

I don't have any problem loading shp and reading its content, it's just the mechanics of solving the task that I can't seem to figure out.

I thought if I iterate over the regions (kind of abusing SetAttributeFilter), get their geometry, then in another loop get geometry of the points and finally do intersect (or maybe within, I am not sure, seemed to un/work the same), I get the result but I must be doing something wrong. Funny thing is, the first result I get is correct but following on it doesn't work as expected.

You can check the crucial part of my code (until here it works fine), I am a beginner so it may be ugly.

filters = []

for code in codes:
    filter1 = str("KOD_KRAJ='" + code + "'")
    filter2 = str('"'+filter1+'"')
    filters.append(filter2)

stations = 0

for filter in filters:
    regions.SetAttributeFilter(filter)

    for i in range(0,numberRe):
        regionFeat = regions.GetNextFeature()
        regionGeometry = regionFeat.GetGeometryRef()

        for r in range(0,numberRa):
            railwayFeature = railway.GetFeature(r)
            railwayGeometry = railwayFeature.GetGeometryRef()

            if railwayGeometry.Within(regionGeometry):
                stations = stations + 1
        print filter, str(stations)

Note: KOD_KRAJ is the header in shapefile's attribute table and variable "codes" is a list containing its instances (administrative code of each region). SetAttributeFilter should be used like this ("NAME='something'") which is why the weird first loop.

I just want to get it working to the stage when I get the numbers for each region, after that it will not be hard for me to do the rest.

Answer 1

You cannot mix GetNextFeature() and a for loop (look at Python GDAL/OGR: retrieved feature is NoneType)

The solution is

for i in range(layer.GetFeatureCount()): #no need of range(0,layer.GetFeatureCount())
    regionFeat = layer.GetFeature(i)
    regionGeometry = regionFeat.geometry()

And you cannot use directly nested loops with OGR (look at Python-OGR: nested loop only loops once)

from osgeo import ogr
poly = ogr.Open("poly_stack.shp")
polys = poly.GetLayer()
point = ogr.Open("pt_stack.shp")
points = point.GetLayer()

If I directly use a nested loop

tot = {}
for i, pol in enumerate(polys):
     #iterate through polygons
     pol_geom = pol.geometry()
     tot[i] = 0
     for j, pt in enumerate(points):
         #iterate through points
         point_geom = pt.geometry()
         if pol_geom.Contains(point_geom):
             tot[i]= tot[i] + 1
 print tot         
 {0: 3, 1: 0, 2: 0}

Only the first loop of points is taken into consideration and you cannot use points.ResetReading() at the end of the loop

If you want only the biggest amount of points in polygons, you can use the Clone() function of OGR (ogr: what is the advantage of using geom.Clone() over geom?)

Make a list of the geometries

pol_geoms =[pol.geometry().Clone() for pol in polys]
point_geoms = [pt.geometry().Clone()  for pt in points] 

Now you can use a nested loop:

tot = {}
#iterate through polygons
for i, pol in enumerate(pol_geoms):
    tot[i] = 0
    for j, pt in enumerate(point_geoms):
         #iterate through points
         if pol.Contains(pt):
             tot[i]= tot[i] + 1
print tot         
{0: 3, 1: 2, 2: 6}

I use here the enumerate function because you can find the original features

But it is much easier with Fiona (More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (1) or GeoPandas (More Efficient Spatial join in Python without QGIS, ArcGIS, PostGIS, etc (2)