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)

2 Answers 2


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, 2012 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. Dec 23, 2012 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, 2012 at 23:45
  • I was cautious, since it depends a lot on how big the datasets are (not specified). There is no silver bullet. Dec 24, 2012 at 8:18

Maybe this will help...


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()
# 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.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.