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.