Vincenty or (Project + Euclidean) for accurate distances

Question

I posted a similar question on stackoverflow, however was directed here for more help.

If my main goal was to calculate distances such that my error was minimised to the nearest metre assuming that the points I am working are bounded on a country-level (for example UK or France ... those kind of medium-sized countries).

Assuming I start with geographic WGS84 co-ordinates would it be better to:

1) Use the vincenty-formula on the epgs:4326 coordinates

2) Project the co-ordinates (e.g. from epsg:4326 to epsg:27700 for the UK) and calculate the hypotenuse

The latter option would be much quicker; however I am not sure which one would be more accurate.

I ran a nearest-neighbour test using (both unvectorised, haversine rather than the more accurate vincenty ... ideally I would compare both functions vectorised using numpy):

nplen = 1000000
# WGS84 lat/long
point = [51.349,-0.19]

# This contains WGS84 lat/long
points = np.ndarray.tolist(np.column_stack(
        [np.round(np.random.randn(nplen)+51,5),
         np.round(np.random.randn(nplen),5)]))

def proj_list(points,
              inproj = Proj(init='epsg:4326'),
              outproj = Proj(init='epsg:27700')):
    """ Projected geo coordinates"""
    return [list(transform(inproj,outproj,x,y)) for y,x in points]

proj_points = proj_list(points)
proj_point = proj_list([point])[0]       

# Run tests
from math import *
def haversine(origin,
              destination):
    """
    Find distance between a pair of lat/lng coordinates
    """
    lat1, lon1, lat2, lon2 = map(radians, [origin[0],origin[1],destination[0],destination[1]])
    dlon = lon2 - lon1
    dlat = lat2 - lat1
    a = sin(dlat / 2) ** 2 + cos(lat1) * cos(lat2) * sin(dlon / 2) ** 2
    c = 2 * asin(sqrt(a))
    r = 6371000  # Metres
    return (c * r)

def closest_math_unproj(points,point):
    """ Haversine on unprojected """
    return (min((haversine(point,pt),pt[0],pt[1]) for pt in points))

def closest_math_proj(points,point):
    """ Simple angle since projected"""
    return (min((hypot(x2-point[1],y2-point[0]),y2,x2) for y2,x2 in points)) 

closest_math_unproj - 22.1 seconds with a distance of 138.491

closest_math_proj - 3.66 seconds with a distnace of 138.904

Where the two points were:

(51.34892,-0.18801)

(51.349,-0.19)

Answer 1

The scale factor – a measure of linear distortion – of national projected coordinates typically (e.g., for UTM) range between 0.9996 and 1.0010. In places this would amount to more than 1m, over large distances.

However, these distortions can be calculated and hence applied. See influence-of-the-scale-factor-on-the-projection, best-formula-for-calculating-short-distances-in-utm, calculating-distance-scale-factor-by-latitude-for-mercator, for example.

Answer 2

Per your points, your baseline is about 100m. I've made tests about this before, and on such short distances it barely makes a difference (in fact, sphere-Earth calculations such as haversine give more error than either).

At this scale, the local topography affects your distance much more than the difference in planimetric measurement methods - if you have Z values accurate enough, calculating topographic distance will be more useful to you than fretting over methods.

Having said that, Vincenty is more accurate. Not only that, your original coordinates are in lat/lon, when you convert them to a projected system, you are imbuing some amount of error in them, it's not a 1:1 translation. The error is pretty small for first order translation, but it's there nonetheless, and it builds up with the distortions that projections naturally give. Unless speed is absolutely crucial to you (say, calculating thousends of baselines for a webservice), I see no reason to convert.

Answer 3

GeographicLib has the highest accuracy for measuring distances on ellipsoids of revolution, with a published accuracy of less than 15 nm.

For example, using GeodSolve to solve the inverse geodesic problem in the question, the distance been the two points (s12) is 138.92830618 m.

Comparision to other distance methods

• Vincenty's method is less accurate, and fails to converge for antipodal points.
• Haversine formula only applies to a sphere, and not (e.g.) to the WGS84 ellipsoid.
• All flat-map projections (e.g. EPSG:27700) have distortions of distance.