NEAR_3D function without the need for 3D Analyst, Arcpy

Question

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

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.