Hide
Know what you're getting – Unlike many sites, all our code is clearly licensed. Join Siafoo Now or Learn More

PL/pgSQL Earth (GPS) distance formula Atom Feed 0

In Brief A function for calculating the distance (in Km) between two GPS coordinates, along with some support functions like Heaver Sin/Arcsin and Earth radius calculation.... more
# 's
 1-- Heaver Sine
2CREATE OR REPLACE FUNCTION heaver_sine(FLOAT) RETURNS FLOAT AS $$
3 SELECT power(sin($1/2.0), 2);
4$$ LANGUAGE SQL;
5
6-- Heaver Arcsine
7CREATE OR REPLACE FUNCTION arc_heaver_sin(FLOAT) RETURNS FLOAT AS $$
8 SELECT 2.0*asin(sqrt($1));
9$$ LANGUAGE SQL;
10
11-- Provides the Earth radius
12CREATE OR REPLACE FUNCTION earth_radius(lat FLOAT) RETURNS float AS $$
13 DECLARE
14 rad_equatorial FLOAT := 6378.137;
15 rad_polar FLOAT := 6356.752;
16 BEGIN
17 RETURN sqrt(power(rad_equatorial, 2) * power(cos(lat), 2) + power(rad_polar, 2)* power(sin(lat), 2)/rad_equatorial*power(cos(lat), 2) + rad_polar*power(sin(lat), 2));
18 END;
19$$ LANGUAGE plpgsql;
20
21-- Give two sets of coordinates in degree form
22CREATE OR REPLACE FUNCTION distance_in_km(lat1 float, lon1 float, lat2 float, long2 float)
23 RETURNS float AS $$
24DECLARE
25 lat1_r FLOAT := radians(lat1);
26 lon1_r FLOAT := radians(lon1);
27 lat2_r FLOAT := radians(lat2);
28 lon2_r FLOAT := radians(lon2);
29BEGIN
30
31 RETURN earth_radius((lat1_r + lat2_r) / 2.0) * arc_heaver_sin(heaver_sin(abs(lat1_r-lat2_r)) + cos(lat1_r)*cos(lat2_r)*heaver_sin(abs(lon1_r - lon2_r)));
32 END;
33$$ LANGUAGE plpgsql;

A function for calculating the distance (in Km) between two GPS coordinates, along with some support functions like Heaver Sin/Arcsin and Earth radius calculation.

You should check out the PostgreSQL Earth Distance contrib module http://www.postgresql.org/docs/8.3/interactive/earthdistance.html.

This is based on A couple of GIS functions in Python.

Comments

over 8 years ago (06 May 2008 at 01:23 PM) by David Isaacson
Awesome.
over 6 years ago (08 Dec 2010 at 08:41 AM) by bricklen
Here is a slightly different version making use of Perl's Math::Trig module, and gives roughly the same results as the contrib Earthdistance & Cube version.

Note: You need to have plperlu installed as a language in your postgresql db.

-- this example using the contrib module "Earthdistance" yields the same results as the Math::Trig value below
select ROUND( ('(33.0, 97.1)'::point <@> '(24.0, 88.6)')::NUMERIC * 1.609344 ) as km;

create or replace function distance_in_km (numeric,numeric,numeric,numeric) returns numeric as
$body$
use strict;
use Math::Trig qw(great_circle_distance deg2rad);
# from http://perldoc.perl.org/Math/Trig.html
my $lat1 = shift(@_);
my $lon1 = shift(@_);
my $lat2 = shift(@_);
my $lon2 = shift(@_);

#elog WARNING, "$lat1, $lon1, $lat2, $lon2\n";

# Notice the 90 - latitude: phi zero is at the North Pole.
sub NESW { deg2rad($_[0]), deg2rad(90 - $_[1]) }
my @L = NESW( $lat1, $lon1 );
my @T = NESW( $lat2, $lon2 );
my $km = great_circle_distance(@L, @T, 6378);

return $km;
$body$ language plperlu strict security definer;

grant execute on function distance_in_km(numeric,numeric,numeric,numeric) to public,nbaffnet;

select ROUND(km) as km from distance_in_km(33.0,97.1,24.0,88.6) as km;
km
-----
945