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 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?

share|improve this question
up vote 4 down vote accepted

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:

## 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)
share|improve this answer
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. – Sirenology 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:

share|improve this answer

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]
    # Convert degrees to radians
    x1 = math.radians(lat1)
    y1 = math.radians(long1)
    x2 = math.radians(lat2)
    y2 = math.radians(long2)
    # 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
share|improve this answer

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.

share|improve this answer
Thanks for the suggestion. However, I can't see how that will help me. According to the docs (…) "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[0], a[1], b[0], b[1])['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)

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.