 License New BSD license
 Lines 29
##### Keywords gis (1) gps distance (1) PostgreSQL (7)
##### Related Loading the GeoIP City database into PostgreSQL Removing duplicates from a PostgreSQL database PostgreSQL script for loading GeoIP data Simple PostgreSQL index test code Utility for loading GeoIP data into a PostgreSQL database
##### Permissions Owner: Stou S. Viewable by Everyone Editable by All Siafoo Users Siafoo – the intersection of pastebin, help desk, version control, and social networking

# PL/pgSQL Earth (GPS) distance formula 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
 Language SQL
# '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;1011-- Provides the Earth radius 12CREATE OR REPLACE FUNCTION earth_radius(lat FLOAT) RETURNS float AS \$\$13    DECLARE14        rad_equatorial FLOAT := 6378.137;15        rad_polar FLOAT := 6356.752;16    BEGIN17    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;2021-- Give two sets of coordinates in degree form22CREATE OR REPLACE FUNCTION distance_in_km(lat1 float, lon1 float, lat2 float, long2 float) 23    RETURNS float AS \$\$24DECLARE25    lat1_r FLOAT := radians(lat1);26    lon1_r FLOAT := radians(lon1);27    lat2_r FLOAT := radians(lat2);28    lon2_r FLOAT := radians(lon2);29BEGIN3031    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.

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;
# 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.