# Calculate distance in km to nearest points (given in lat/long) using ArcGIS DEsktop and/or R?

I have two point datasets in ArcGIS, both of which are given in WGS84 lat/lon co-ordinates and the points are spread across the whole world. I would like to find the nearest point in Dataset A to each point in Dataset B, and get the distance between them in kilometres.

This seems like a perfect use of the Near tool, but that gives me results in the co-ordinate system of the input points: that is, decimal degrees. I know I could re-project the data, but I gather (from this question) that it is difficult (if not impossible) to find a projection that will give accurate distances all over the world.

The answers to that question suggest using the Haversine formula to calculate distances using the latitude-longitude co-ordinates directly. Is there a way to do this and get a result in km using ArcGIS? If not, what is the best way to approach this?

Although this isn't an ArcGIS solution, your problem can be solved in R by exporting your points from Arc and using the `spDists` function from the `sp` package. The function finds the distances between a reference point(s) and a matrix of points, in kilometers if you set `longlat=T`.

Here's a quick and dirty example:

``````library(sp)
## Sim up two sets of 100 points, we'll call them set a and set b:
a <- SpatialPoints(coords = data.frame(x = rnorm(100, -87.5), y = rnorm(100, 30)), proj4string=CRS("+proj=longlat +datum=WGS84"))
b <- SpatialPoints(coords = data.frame(x = rnorm(100, -88.5), y = rnorm(100, 30.5)), proj4string=CRS("+proj=longlat +datum=WGS84"))

## Find the distance from each point in a to each point in b, store
##    the results in a matrix.
results <- spDists(a, b, longlat=T)
``````
• Thanks - this seems like the most realistic solution. Looking at the docs it seems that I can only do this between a reference point and a set of other points, so I would have to do it in a loop to go through all of my points. Do you know of a more efficient way to do this in R? – robintw Sep 14 '12 at 10:07
• No looping needed, you can give the function two sets of points and it will return a matrix with distances between each combination of points. Edited answer to include example code. – allen Sep 17 '12 at 22:42

It's not an ArcGIS solution, but using a Round Earth data model in a spatial database would do the trick. Calculating earth distance in database supporting this would be pretty easy. I can suggest you two readings:

http://postgis.net/workshops/postgis-intro/geography.html

http://blog.safe.com/2012/08/round-earth-data-in-oracle-postgis-and-sql-server/

You need a distance calculation that works with Lat/Long. Vincenty is the one I would use (0.5mm accuracy). I have played with it before, and it is not too hard to use.

The code is a bit long, but it works. Given two points in WGS, it will return a distance in meters.

You can use this as a Python script in ArcGIS, or wrap it around another script that simply iterates over the two Point Shapefiles and builds a distance matrix for you. Or, it is probably easier to feed the results of GENERATE_NEAR_TABLE with finding the 2-3 nearest features (to avoid complications of earth's curvature).

``````import math

ellipsoids = {
#name        major(m)   minor(m)            flattening factor
'WGS-84':   (6378137,   6356752.3142451793, 298.25722356300003),
'GRS-80':   (6378137,   6356752.3141403561, 298.25722210100002),
'GRS-67':   (6378160,   6356774.5160907144, 298.24716742700002),

}

def distanceVincenty(lat1, long1, lat2, long2, ellipsoid='WGS-84'):
"""Computes the Vicenty distance (in meters) between two points
on the earth. Coordinates need to be in decimal degrees.
"""
# Check if we got numbers
# Removed to save space
# Check if we know about the ellipsoid
# Removed to save space
major, minor, ffactor = ellipsoids[ellipsoid]
# Define our flattening f
f = 1 / ffactor
# Find delta X
deltaX = y2 - y1
# Calculate U1 and U2
U1 = math.atan((1 - f) * math.tan(x1))
U2 = math.atan((1 - f) * math.tan(x2))
# Calculate the sin and cos of U1 and U2
sinU1 = math.sin(U1)
cosU1 = math.cos(U1)
sinU2 = math.sin(U2)
cosU2 = math.cos(U2)
# Set initial value of L
L = deltaX
# Set Lambda equal to L
lmbda = L
# Iteration limit - when to stop if no convergence
iterLimit = 100
while abs(lmbda) > 10e-12 and iterLimit >= 0:
# Calculate sine and cosine of lmbda
sin_lmbda = math.sin(lmbda)
cos_lmbda = math.cos(lmbda)
# Calculate the sine of sigma
sin_sigma = math.sqrt(
(cosU2 * sin_lmbda) ** 2 +
(cosU1 * sinU2 -
sinU1 * cosU2 * cos_lmbda) ** 2
)
if sin_sigma == 0.0:
# Concident points - distance is 0
return 0.0
# Calculate the cosine of sigma
cos_sigma = (
sinU1 * sinU2 +
cosU1 * cosU2 * cos_lmbda
)
# Calculate sigma
sigma = math.atan2(sin_sigma, cos_sigma)
# Calculate the sine of alpha
sin_alpha = (cosU1 * cosU2 * math.sin(lmbda)) / (sin_sigma)
# Calculate the square cosine of alpha
cos_alpha_sq = 1 - sin_alpha ** 2
# Calculate the cosine of 2 sigma
cos_2sigma = cos_sigma - ((2 * sinU1 * sinU2) / cos_alpha_sq)
# Identify C
C = (f / 16.0) * cos_alpha_sq * (4.0 + f * (4.0 - 3 * cos_alpha_sq))
# Recalculate lmbda now
lmbda = L + ((1.0 - C) * f * sin_alpha * (sigma + C * sin_sigma * (cos_2sigma + C * cos_sigma * (-1.0 + 2 * cos_2sigma ** 2))))
# If lambda is greater than pi, there is no solution
if (abs(lmbda) > math.pi):
raise ValueError("No solution can be found.")
iterLimit -= 1
if iterLimit == 0 and lmbda > 10e-12:
raise ValueError("Solution could not converge.")
# Since we converged, now we can calculate distance
# Calculate u squared
u_sq = cos_alpha_sq * ((major ** 2 - minor ** 2) / (minor ** 2))
# Calculate A
A = 1 + (u_sq / 16384.0) * (4096.0 + u_sq * (-768.0 + u_sq * (320.0 - 175.0 * u_sq)))
# Calculate B
B = (u_sq / 1024.0) * (256.0 + u_sq * (-128.0 + u_sq * (74.0 - 47.0 * u_sq)))
# Calculate delta sigma
deltaSigma = B * sin_sigma * (cos_2sigma + 0.25 * B * (cos_sigma * (-1.0 + 2.0 * cos_2sigma ** 2) - 1.0/6.0 * B * cos_2sigma * (-3.0 + 4.0 * sin_sigma ** 2) * (-3.0 + 4.0 * cos_2sigma ** 2)))
# Calculate s, the distance
s = minor * A * (sigma - deltaSigma)
# Return the distance
return s
``````

I made similar experiences with small datasets using the Point Distance tool. Doing so, you cannot automatically find the nearest points in your Dataset A, but at least get a table output with useful km or m results. In a next step you could select the shortest distance to each point of Dataset B out of the table.

But this approach would depend on the amount of points in your datasets. It might not work properly with large datasets.

• Thanks for the suggestion. However, I can't see how that will help me. According to the docs (help.arcgis.com/en/arcgisdesktop/10.0/help/index.html#//…) "the distance is in the linear unit of the input features coordinate system.", which as my input features are in lat/lon will surely give me results in decimal degrees? (I haven't got a machine with ArcGIS on it here to test) – robintw Sep 7 '12 at 22:29
• In this case, I would probably use a "quick and dirty" solution by adding X and Y fields in your datatable and click on Calculate Geometry choosing X and Y in meter. If not possible to pick this option, change the coordinate system of your MXD. I was working on a project before, where my client wanted long/lat, X/Y and Gauss-Krueger R/H values all in each Shape file. To avoid complicated calculation, simply changing projections and calculate geometry was the most easy way it worked out. – basto Sep 10 '12 at 6:45

If you need high-precision and robust geodesic measurements, use GeographicLib, which is natively written in several programming languages, including C++, Java, MATLAB, Python, etc.

See C. F. F. Karney (2013) "Algorithms for geodesics" for a literary reference. Note that these algorithms are more robust and accurate than Vincenty's algorithm, for instance near antipodes.

To calculate distance in metres between two points, get the `s12` distance attribute from the inverse geodesic solution. E.g., with the geographiclib package for Python

``````from geographiclib.geodesic import Geodesic
g = Geodesic.WGS84.Inverse(-41.32, 174.81, 40.96, -5.50)
print(g)  # shows:
{'a12': 179.6197069334283,
'azi1': 161.06766998615873,
'azi2': 18.825195123248484,
'lat1': -41.32,
'lat2': 40.96,
'lon1': 174.81,
'lon2': -5.5,
's12': 19959679.26735382}
``````

Or make a convenience function, which also converts from metres to kilometres:

``````dist_km = lambda a, b: Geodesic.WGS84.Inverse(a, a, b, b)['s12'] / 1000.0
a = (-41.32, 174.81)
b = (40.96, -5.50)
print(dist_km(a, b))  # 19959.6792674 km
``````

Now to find the closest point between lists `A` and `B`, each with 100 points:

``````from random import uniform
from itertools import product
A = [(uniform(-90, 90), uniform(-180, 180)) for x in range(100)]
B = [(uniform(-90, 90), uniform(-180, 180)) for x in range(100)]
a_min = b_min = min_dist = None
for a, b in product(A, B):
d = dist_km(a, b)
if min_dist is None or d < min_dist:
min_dist = d
a_min = a
b_min = b

print('%.3f km between %s and %s' % (min_dist, a_min, b_min))
``````

22.481 km between (84.57916462672875, 158.67545706102192) and (84.70326937581333, 156.9784597422855)