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
SET @x = `employees`.`vincenty`(45.0, 0.0, 45.0, 1.0);
SELECT @x;
This yields the answer .707…
I would expect it to yield 1.0
Your example fetches the distance from the point at 45°N, 0°E (a point just to the west of Bordeaux, France) to 45°N, 1°E. When you’re at 45°N or S, and you move one degree east, you don’t go the full 60 nautical miles east. Instead you go 60 * (1/ √2) nm east.
If your example were the distance from 45°N, 0°E to 65°N, 0°E you would go the full 60nm.
thanks for sharing this article. grate work!!.
i want to know how change this SP to return distance from Meters?
This function returns its results in degrees, not km, miles, furlongs, parsecs or any other distance measure.
To convert to km, multiply your degrees by 111.045. For statute miles, use 69.0. For nautical miles, use 60.0
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
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.
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;
Can you add another sample code of how you would use this function in a SQL query please?
Take a look at this post. http://www.plumislandmedia.net/mysql/stored-function-haversine-distance-computation/. It does what you ask using a different stored function with the same formal parameters.