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

Question

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?

Answer 1

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)

Answer 2

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/

Answer 3

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

Answer 4

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.

Answer 5

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)

Answer 6

Just in case somebody wants to exchange the existing sp solution for one making use of sf:

library(sf)
#> Linking to GEOS 3.9.1, GDAL 3.4.3, PROJ 7.2.1; sf_use_s2() is TRUE

set.seed(42)

# create two sets of 100 points each
a <- data.frame(x = rnorm(100, -87.5), y = rnorm(100, 30)) |> st_as_sf(coords = c("x", "y"), crs = 4326)
b <- data.frame(x = rnorm(100, -88.5), y = rnorm(100, 30.5)) |> st_as_sf(coords = c("x", "y"), crs = 4326)

# get the full distance matrix
results <- st_distance(a, b)