The Vincenty great-circle distance formula

This Vincenty formula is a more numerically stable version of the spherical cosine law formula (commonly and wrongly known as the Haversine formula) for computing great circle distances. The question of numerical stability comes up specifically when the distances between points are small. In those cases the cosine is very close to 1, so the inverse cosine function is not as precise as the inverse tangent function used in this formula.

See http://en.wikipedia.org/wiki/Great-circle_distance for more information.

DELIMITER $$
DROP FUNCTION IF EXISTS vincenty$$
CREATE FUNCTION vincenty(
        lat1 FLOAT, lon1 FLOAT,
        lat2 FLOAT, lon2 FLOAT
     ) RETURNS FLOAT
    NO SQL
    DETERMINISTIC
    COMMENT 'Returns the distance in degrees on the
             Earth between two known points
             of latitude and longitude
             using the Vincenty formula
             from http://en.wikipedia.org/wiki/Great-circle_distance'
BEGIN
    RETURN  DEGREES(
    ATAN2(
      SQRT(
        POW(COS(RADIANS(lat2))*SIN(RADIANS(lon2-lon1)),2) +
        POW(COS(RADIANS(lat1))*SIN(RADIANS(lat2)) -
             (SIN(RADIANS(lat1))*COS(RADIANS(lat2)) *
              COS(RADIANS(lon2-lon1))) ,2)),
      SIN(RADIANS(lat1))*SIN(RADIANS(lat2)) +
      COS(RADIANS(lat1))*COS(RADIANS(lat2))*COS(RADIANS(lon2-lon1))));
END$$
DELIMITER ;

Vincenty formula for Oracle PL/SQL

Here is the same formula written in the Oracle stored function language.

create or replace FUNCTION vincenty  (
  lat1 IN REAL, lon1 IN REAL, 
  lat2 IN REAL, lon2 IN REAL)
RETURN REAL
DETERMINISTIC 
IS

  DEGREES REAL := 57.2957795;
  RADIANS REAL := 0.0174532925;

BEGIN
 RETURN  DEGREES * (
    ATAN2(
      SQRT(
        POWER(COS(RADIANS * (lat2))*SIN(RADIANS * (lon2-lon1)),2) +
        POWER(COS(RADIANS * (lat1))*SIN(RADIANS * (lat2)) -
             (SIN(RADIANS * (lat1))*COS(RADIANS * (lat2)) *
              COS(RADIANS * (lon2-lon1))) ,2)),
      SIN(RADIANS * (lat1))*SIN(RADIANS * (lat2)) +
      COS(RADIANS * (lat1))*COS(RADIANS * (lat2))*COS(RADIANS * (lon2-lon1))));
END;

Vincenty formula in T-SQL for Microsoft SQL server

Here is the same formula written in the T-SQL stored function language.

IF OBJECT_ID (N'dbo.Vincenty', N'FN') IS NOT NULL
    DROP FUNCTION Vincenty;
GO
CREATE FUNCTION dbo.Vincenty(
  @lat1 float, @lon1 float,
  @lat2 float, @lon2 float )
RETURNS float 
AS 
-- Returns distance in degrees
BEGIN
   RETURN  DEGREES(
    ATN2(
      SQRT(
        POWER(COS(RADIANS(@lat2))*SIN(RADIANS(@lon2-@lon1)),2) +
        POWER(COS(RADIANS(@lat1))*SIN(RADIANS(@lat2)) -
             (SIN(RADIANS(@lat1))*COS(RADIANS(@lat2)) *
              COS(RADIANS(@lon2-@lon1))) ,2)),
      SIN(RADIANS(@lat1))*SIN(RADIANS(@lat2)) +
      COS(RADIANS(@lat1))*COS(RADIANS(@lat2))*COS(RADIANS(@lon2-@lon1))));
END;
GO

  5 comments for “The Vincenty great-circle distance formula

  1. Fred
    April 10, 2019 at 6:09 am

    Hi Ollie,
    Thank you sharing your valuable knowledge.

    I use to use the Haversine formula but its output was never accurate. I started to look at the Vicenty formula until I came across your post. The plSql Vicenty function you’re sharing returns the distance in latitude degrees, right? if yes, then once I convert that to KM it doesn’t give the correct km distance for all the points I am looking for.

    I believe that the Vicenty formula is a bit more complex than this and involves an itereation at the end in order to come up with the extact distance. Could you please comment on this?

    NB: I am using this for airports coordinates

    Regards,
    Fred

    • Ollie
      April 10, 2019 at 2:59 pm

      Can you characterize the inaccuracy? Is it Tierra del Fuego inaccurate, or across-the-street inaccurate?

      None of these great-circle formulae are suitable for high-accuracy applications like IFR navigation or surveying. They assume the Earth is spherical, which is OK for store finders but not for high accuracy: The Earth’s radius at the equator is a little bit larger than it is at the poles.

      Why is Vincenty’s formula better than the cosine law (haversine) formula? Because the cosine law depends on arccos(x) where x is near 1. The smaller the distance involved, the nearer 1 it gets. So, tiny differences in x get magnified by arccos(). Also, the way the rest of the formula works, it sometimes can generate x values slightly larger than 1, which is out of the domain of arccos().

      On the other hand, Vincenty’s formula depends on arctan (y,x), which doesn’t have that problem.

      • Fred
        April 11, 2019 at 12:21 pm

        Hi Ollie,

        Thank you for your quick feedback.
        Below is the script am using to calculate the distance between 2 airports using Haversine formula.
        I use the distance to then calculate a “Circuity factor” which is the ratio bewteen Dist[A to B]+Dist[B to C] / Dist[A to C].
        Now when I compare the distance that I get vs. the one from greatcirclemapper.net/ or travelmath.com/ the gap is a bit small for short distances but becomes big for long ones.
        I have the Vincenty formula in lua script but I’ll need to convert it and test it in plsql and share the findings with you.

        SELECT ROUND(GCMAP_DIST(‘LHR’,’AMS’),4) D FROM DUAL; –(369.0651 vs 369.77)
        SELECT ROUND(GCMAP_DIST(‘LHR’,’CDG’),4) D FROM DUAL; –(347.1075 vs 347.20)
        SELECT ROUND(GCMAP_DIST(‘LHR’,’HND’),4) D FROM DUAL; –(9601.9238 vs 9614.78)

        CREATE OR REPLACE FUNCTION GCMAP_DIST(ORIG VARCHAR2,DEST VARCHAR2) RETURN FLOAT AUTHID CURRENT_USER AS

        LAT1 NUMBER; LON1 NUMBER; LAT2 NUMBER; LON2 NUMBER;
        EARTHRADIUS NUMBER:=6378.137; VPI NUMBER:= 3.14159265359;

        BEGIN
        SELECT CASE WHEN COORDINATE1=’N’ THEN (VPI*(SUBSTR(LATITUDE,1,2)+(SUBSTR(LATITUDE,4,2)/60))/180)
        ELSE -1*(VPI*(SUBSTR(LATITUDE,1,2)+(SUBSTR(LATITUDE,4,2)/60))/180) END INTO LAT1 FROM STATION WHERE STATION =ORIG;
        SELECT CASE WHEN COORDINATE2=’E’ THEN (VPI*(SUBSTR(LONGITUDE,1,3)+(SUBSTR(LONGITUDE,5,2)/60))/180)
        ELSE -1*(VPI*(SUBSTR(LONGITUDE,1,3)+(SUBSTR(LONGITUDE,5,2)/60))/180) END INTO LON1 FROM STATION WHERE STATION =ORIG;
        SELECT CASE WHEN COORDINATE1=’N’ THEN (VPI*(SUBSTR(LATITUDE,1,2)+(SUBSTR(LATITUDE,4,2)/60))/180)
        ELSE -1*(VPI*(SUBSTR(LATITUDE,1,2)+(SUBSTR(LATITUDE,4,2)/60))/180) END INTO LAT2 FROM STATION WHERE STATION =DEST;
        SELECT CASE WHEN COORDINATE2=’E’ THEN (VPI*(SUBSTR(LONGITUDE,1,3)+(SUBSTR(LONGITUDE,5,2)/60)) /180)
        ELSE -1*(VPI*(SUBSTR(LONGITUDE,1,3)+(SUBSTR(LONGITUDE,5,2)/60))/180) END INTO LON2 FROM STATION WHERE STATION =DEST;

        RETURN 2*EARTHRADIUS*ASIN(SQRT((POWER(SIN((LAT1 – LAT2) / 2),2) + COS(LAT1)*COS(LAT2)*(POWER(SIN((LON1 – LON2) / 2) , 2)))));
        END GCMAP_DIST;

  2. Andrew
    June 13, 2014 at 6:54 pm

    Can you add another sample code of how you would use this function in a SQL query please?

Leave a Reply

Your email address will not be published. Required fields are marked *