Efficient computation of nearest distance for each point to closest polygon

Question

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)

Answer 1

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.

Answer 2

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')
prj_txt = prj_file.read()
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)

shp_file = shapefile.Reader(shapefile_path_without_dot_shape)
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[1])
    lons = np.append(lons, np.array(coordTrans.TransformPoints(shapes[i].points)).T[0])

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)

min_dist = distances.min() # your answer

Now granted, your shape file needs a high resolution of coordinates to be fairly accurate. And this would certainly be better built as a function, but hopefully its of use to somebody.