# Vincenty or (Project + Euclidean) for accurate distances

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)

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.

• I see - so it would appear that projecting makes it computationally easier (compared to haversine or vincenty); however sacrifices some accuracy Commented Feb 8, 2016 at 10:29

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.