Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I am looking for a python script which would do similar as the ArcGIS 3D Analyst NEAR_3D function.

Basically I would like to programatically in python find nearest points of one point feature class to another point feature class and calculate their delta X, Y, Z.

Any help much appreciated!

Thanks meryloo

share|improve this question
add comment

1 Answer 1

You are doing what is known as a nearest neighbor search. This is quite easy to do iteratively: for each point, loop through all candidate neighbors and calculate their distance. What's challenging is how to make it efficient when you want to get the N nearest neighbors of P points and N and P are very large.

In my experience, the best balance of complexity and performance is to eschew arcpy algother and just use the spatial.KDTree module in scipy. A k-d tree is a space partioning concept that subdivides k-dimensional space into regions that drastically reduces the number of distance formula calculations required in the nearest neighbor search. If the nearest corner of Region A is farther away from the starting point than the farthest point in Region B, then all points in Region B must be closer to the starting point than all points in Region A. Therefore, we shouldn't bother calculating a distance formula from the starting point to points in Region A.

Here's a code sample. I'm not at work so I can't test this, but it's adapted from a project where I had to find the 100 nearest neighbors to approximately 400,000 2-dimensional points. The total processing time with cKDTree was approximately 4 minutes.

import numpy as np
#Better to use cKDTree when not doing complicated queries
#from scipy.spatial import KDTree
from scipy.spatial import cKDTree as KDTree

#Input data
MyData = [[0,0],[1,6],[2,9]]

#Convert data to numpy array, requirement for scipy classes
ArrayOfPoints = np.array(MyData)

#Create KDTree from numpy array
KDTreeOfPoints = KDTree(ArrayOfPoints)

#Iterate through points
for PointNum in range(len(MyData)):

    #Query KDTree, return distance to nearest n points and their point numbers
    #Need to have second parameter = 2 because "closest" neighbor in KDTree is itself
    Distances,Indices = KDTreeOfPoints.query(ArrayOfPoints[PointNum],2)

    #Get info about nearest neighbor.
    CurrentPoint = str(MyData[PointNum])
    NearestNeighbor = str(MyData[Indices[1])
    NeighborDist = float(Distances[1])

    print "Nearest to %s is %s, dist = %f" % (CurrentPoint,NearestNeighbor,NeighborDist)

Note that scipy is not part of the standard python library. It can be downloaded here. Use the "superpack" installer for windows as you probably don't want to build the source.

share|improve this answer
    
many thanks - this is what I am looking for. I will need to add a 3rd dimension to it, but it is a really good start for me.meryloo –  meryloo Nov 26 '12 at 11:44
add comment

Your Answer

 
discard

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.