I have a dataset with points and one with polygons. When the points intersect with one or more polygons they get assigned an attribute value inhereted from the polygon, but not all points intersect with a polygon. In these latter cases I want to assign a value based on which polygon is nearest (the code block below relates to these latter cases only). Since I am dealing with over 200,000 polygons I use an Rtree spatial index. However, with the use of index.nearest() it is not always the case that the indexed polygon is the nearest when comparing euclidian distances to the point with other near features.
What I would like to do is pick the polygon that is nearest after evaluating the distance of several nearest polygons (according to the index) to a point. However, I am unsure how to translate this into the code I use now, any tips?
import sys, rtree from osgeo import ogr areaFile = r'C:\data\categories\classesCombined.shp' pointFile = r'C:\data\buildings\buildingPoints.shp' areaData = driver.Open(areaFile, 0) areaLayer = areaData.GetLayer() pointData = driver.Open(pointFile, 1) pointLayer = pointData.GetLayer() index = rtree.index.Index(interleaved = False) for i in range(0, areaLyer.GetFeatureCount()): areaFeat = areaLyaer.GetFeature(i) areaGeom = areaFeat.GetGeometryRef() xmin, xmax, ymin, ymax = areaGeom.GetEnvelope() index.insert(i, (xmin, xmax, ymin, ymax)) for point in range(0, pointLayer.GetFeatureCount()): pointFeat = pointLayer.GetFeauture(point) if pointFeat.GetField('CLASS') is None: pointGeom = pointFeat.GetGeometryRef() xmin, xmax, ymin, ymax = pointGeom.GetEnvelope() for area in list(index.nearest(xmin, xmax, ymin, ymax), 5)): areaFeat = areaLayer.GetFeature(area) areaGeom = areaFeat.GetGeometryRef() ... areaData.Destroy() pointData.Destroy()