Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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)
share|improve this question

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.

share|improve this answer
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...

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 =
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.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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