# Efficient computation of nearest distance for each point to closest polygon

I'm looking for a more efficient way to compute the nearest distance between features from 2 layers, in this example for each point the nearest distance to the closest polygon. Is there a more efficient way compared to below where I compute all polygon distances for each point? (Python 2.7.2, GDAL/OGR 1.9.1)

``````nPts = pointsLayer.GetFeatureCount()
nPolys = polysLayer.GetFeatureCount()

for pt in range(0,nPts):
minDist = 1000000
pointFeature = pointsLayer.GetFeature(pt)
for poly in range(0,nPolys):
polyFeature = polysLayer.GetFeature(poly)
Dist = pointFeature.GetGeometryRef().Distance( polyFeature.GetGeometryRef() )
if (Dist < minDist):
minDist = Dist
print (pt,minDist)
``````

The first trivial optimisation that comes to mind is storing `pointFeature.GetGeometryRef()` in a temporary variable instead of repeating the lookup for each inner loop iteration.

But that is likely completely insignificant. Since you want just the distance to the nearest polygon, you don't need two loops. Make a union out of all the polygons / dissolve the layer and then check the distance against that.

• I'm looking for a more generic efficient solution, for example what if both layers are point feature type? Thanks. – Dave Dec 23 '12 at 13:58
• Well, why didn't you ask that way then? In that case you can /maybe/ achieve a speed gain by buffering the points incrementally and checking if you've made an intersection with the other layer. If so, do a range check only on the contained points. Playing with a minimal nQuadSegs value can improve the check speed, but then you need to do it twice to make sure you get the right point, since the diagonal of a square is longer than the side. – lynxlynxlynx Dec 23 '12 at 14:26
• Sorry I was not clear in my original post. However, you answer including "Maybe", "Playing with", "do it twice" gives me the impression that this may or may not help..please supply example code that you know is definitely substantially more efficient. THANK YOU! – Dave Dec 23 '12 at 23:45
• I was cautious, since it depends a lot on how big the datasets are (not specified). There is no silver bullet. – lynxlynxlynx Dec 24 '12 at 8:18

Maybe this will help...

https://stackoverflow.com/questions/13489835/gps-positioning-with-python-geopy-shapely

I would be nice if there were a GDAL function for this. What I did for my self (not needing a high level of efficiency) was convert the shape file x, y coordinates to lat, lon coordinates using osr. Then I looped through the polygons and calculated the vincenty distances using geopy. for example...

``````import numpy as np
import osr
import shapefile
from geopy.distance import vincenty

# get the source projection parameters
prj_file = open(path_to_dot_prj_file, 'r')
sourceSR = osr.SpatialReference()
sourceSR.ImportFromESRI([prj_txt])
# get the destination projection parameters
targetSR = osr.SpatialReference()
targetSR.ImportFromEPSG(4326) # or whatever spatial reference you want

coordTrans = osr.CoordinateTransformation(sourceSR,targetSR)

shapes = shp_file.shapes()

lats = []
lons = []
for i in range(0, len(shapes) - 1):
lats = np.append(lats, np.array(coordTrans.TransformPoints(shapes[i].points)).T)
lons = np.append(lons, np.array(coordTrans.TransformPoints(shapes[i].points)).T)

point_coord = (30.0, 90.0) # could an array of locations if you modified code

distances = []
for i in range(0,len(lats) - 1):
distances = np.append(distances, vincenty(point_coord, (lats[i],lons[i])).km)