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 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
add comment

5 Answers

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:

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)
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
add comment

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
add comment

To calculate distances you need a equidistant map projection.

The problem is: There is no projection for the whole world which is everywhere/in all direction equidistant. Here's a good explanation of this: http://www.progonos.com/furuti/MapProj/Normal/CartProp/DistPres/distPres.html

You can split the data in different regions. Reproject the data of each area into a equidistant map projection which is good for that region. Make the distance calculation in each area.

share|improve this answer
    
Yes, that seems to be the only way to do it in ArcGIs. However, it seems that there are mathematical ways to do it that use the geometry of a sphere (or geoid, depending how accurate they are) that are better ways to do it, and don't require splitting the data up. It seems strange that some of these more standard geometry calculations can't be done in ArcGIS. –  robintw Sep 14 '12 at 10:03
    
(-1) This solution is not a good one because equidistant projections tend to be really bad for general-purpose distance calculations. They are "equidistant" only relative to a small number of fixed base points (usually just one or two). In order to give them this property, the projections usually introduce a lot of distortion between other pairs of points. –  whuber Sep 18 '12 at 16:27
add comment

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://workshops.opengeo.org/postgis-intro/geography.html

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

share|improve this answer
add comment

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 (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
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.