In marketing we often use distance to the buyer as a parameter for data mining. There are plenty of ways and sources how to calculate it, including Google Maps, Yandex Maps, spatial databases etc. Here I present some solution to calculate it having only longitude and latitude of the place.
My data are:
These are the coordinates of New York and Washington. I want to calculate the distance between these cities as easy as possible and as accurate as possible.
At once let's install the earthdistance PostreSQL contrib module.
Code:
| CREATE EXTENSION cube; |
| CREATE EXTENSION earthdistance; |
Secondly, let's check the custom function of Stou S. from here.
Code:
| -- Heaver Sine |
| CREATE OR REPLACE FUNCTION heaver_sin(FLOAT) |
| RETURNS FLOAT |
| AS |
| $$ |
| SELECT power(sin($1/2.0), 2); |
| $$ |
| LANGUAGE SQL; |
| |
| -- Heaver Arcsine |
| CREATE OR REPLACE FUNCTION arc_heaver_sin(FLOAT) |
| RETURNS FLOAT |
| AS |
| $$ |
| SELECT 2.0*asin(sqrt($1)); |
| $$ |
| LANGUAGE SQL; |
| |
| -- Provides the Earth radius |
| CREATE OR REPLACE FUNCTION earth_radius(lat FLOAT) |
| RETURNS float |
| AS |
| $$ |
| DECLARE rad_equatorial FLOAT := 6378.137; |
| rad_polar FLOAT := 6356.752; |
| BEGIN |
| 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)); |
| END; |
| $$ |
| LANGUAGE plpgsql; |
| |
| -- Give two sets of coordinates in degree form |
| CREATE OR REPLACE FUNCTION distance_in_km(lat1 float, lon1 float, lat2 float, lon2 float) |
| RETURNS float |
| AS |
| $$ |
| DECLARE lat1_r FLOAT := radians(lat1); |
| lon1_r FLOAT := radians(lon1); |
| lat2_r FLOAT := radians(lat2); |
| lon2_r FLOAT := radians(lon2); |
| BEGIN 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))); |
| END; |
| $$ |
| LANGUAGE plpgsql; |
Thirdly, we can create our own function based of the very well-known expression.
Code:
| CREATE OR REPLACE FUNCTION CalculateDistance(lat1 float8, long1 float8, lat2 float8, long2 float8) |
| RETURNS FLOAT |
| AS |
| $$ |
| SELECT atan2((sqrt(power(cos(lat2 * pi() / 180) * sin(long2 * pi() / 180 - long1 * pi() / 180), 2) + power(cos(lat1 * pi() / 180) * sin(lat2 * pi() / 180) - sin(lat1 * pi() / 180) * cos(lat2 * pi() / 180) * cos(long2 * pi() / 180 - long1 * pi() / 180), 2))), (sin(lat1 * pi() / 180) * sin(lat2 * pi() / 180) + cos(lat1 * pi() / 180) * cos(lat2* pi() / 180) * cos(long2 * pi() / 180 - long1 * pi() / 180))) * 6372.795477598; |
| $$ |
| LANGUAGE SQL; |
Sorry if not comprehensive but I tried not to use PL/pgSQL. First you must notice is the different values of Earth radius which dramatically influences on the final result. But let's calculate the distance at last.
Code:
| SELECT ROUND( ('(-74.00583, 40.71417)'::point'(-77.03611, 38.895)')::NUMERIC * 1.609344 ) AS "km (earthdistance extension)" |
| , distance_in_km(40.71417, -74.00583, 38.895, -77.03611) AS "km (distance_in_km custom function)" |
| , CalculateDistance(40.71417, -74.00583, 38.895, -77.03611) AS "km (CalculateDistance custom function)"; |
The results are: 328 (earthdistance module), 252.65 (distance_in_km function) and 328.57(CalculateDistance function). The second function looks like unreliable but you always can edit it.
Now I am moving to SQL SERVER and I want to use the embedded STDistance function for this task. The MSDN says it "returns the shortest distance between a point in a geometry instance and a point in another geometry instance".
Code:
| DECLARE @g geography,@h geography |
| SET @g=geography::Point(-74.00583,40.71417,4326) |
| SET @h=geography::Point(-77.03611,38.895,4326) |
| SELECT @h.STDistance(@g)/1000 AS sqlserver_stdistance |
| GO |
The result is 342, which is slightly higher than the other functions' yield.
Google Maps uses roads to calculate the distance and shows 368.3.
That's at no means the full list of possible function but I hope it helps.