Skip to content

Instantly share code, notes, and snippets.

@ywjno
Created February 4, 2015 04:18
Show Gist options
  • Save ywjno/10916f389f8733127a01 to your computer and use it in GitHub Desktop.
Save ywjno/10916f389f8733127a01 to your computer and use it in GitHub Desktop.
Calculate distance between two points on a globe(http://www.codecodex.com/wiki/Calculate_distance_between_two_points_on_a_globe)
/// @brief The usual PI/180 constant
static const double DEG_TO_RAD = 0.017453292519943295769236907684886;
/// @brief Earth's quatratic mean radius for WGS-84
static const double EARTH_RADIUS_IN_METERS = 6372797.560856;
/** @brief Computes the arc, in radian, between two WGS-84 positions.
*
* The result is equal to <code>Distance(from,to)/EARTH_RADIUS_IN_METERS</code>
* <code>= 2*asin(sqrt(h(d/EARTH_RADIUS_IN_METERS )))</code>
*
* where:<ul>
* <li>d is the distance in meters between 'from' and 'to' positions.</li>
* <li>h is the haversine function: <code>h(x)=sin²(x/2)</code></li>
* </ul>
*
* The haversine formula gives:
* <code>h(d/R) = h(from.lat-to.lat)+h(from.lon-to.lon)+cos(from.lat)*cos(to.lat)</code>
*
* @sa http://en.wikipedia.org/wiki/Law_of_haversines
*/
double ArcInRadians(const Position& from, const Position& to) {
double latitudeArc = (from.lat - to.lat) * DEG_TO_RAD;
double longitudeArc = (from.lon - to.lon) * DEG_TO_RAD;
double latitudeH = sin(latitudeArc * 0.5);
latitudeH *= latitudeH;
double lontitudeH = sin(longitudeArc * 0.5);
lontitudeH *= lontitudeH;
double tmp = cos(from.lat*DEG_TO_RAD) * cos(to.lat*DEG_TO_RAD);
return 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
}
/** @brief Computes the distance, in meters, between two WGS-84 positions.
*
* The result is equal to <code>EARTH_RADIUS_IN_METERS*ArcInRadians(from,to)</code>
*
* @sa ArcInRadians
*/
double DistanceInMeters(const Position& from, const Position& to) {
return EARTH_RADIUS_IN_METERS*ArcInRadians(from, to);
}
import com.google.android.maps.GeoPoint;
public class DistanceCalculator {
private double Radius;
// R = earth's radius (mean radius = 6,371km)
// Constructor
DistanceCalculator(double R) {
Radius = R;
}
public double CalculationByDistance(GeoPoint StartP, GeoPoint EndP) {
double lat1 = StartP.getLatitudeE6()/1E6;
double lat2 = EndP.getLatitudeE6()/1E6;
double lon1 = StartP.getLongitudeE6()/1E6;
double lon2 = EndP.getLongitudeE6()/1E6;
double dLat = Math.toRadians(lat2-lat1);
double dLon = Math.toRadians(lon2-lon1);
double a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(Math.toRadians(lat1)) * Math.cos(Math.toRadians(lat2)) *
Math.sin(dLon/2) * Math.sin(dLon/2);
double c = 2 * Math.asin(Math.sqrt(a));
return Radius * c;
}
}
# haversine.rb
#
# haversine formula to compute the great circle distance between two points given their latitude and longitudes
#
# Copyright (C) 2008, 360VL, Inc
# Copyright (C) 2008, Landon Cox
#
# http://www.esawdust.com (Landon Cox)
# contact:
# http://www.esawdust.com/blog/businesscard/businesscard.html
#
# LICENSE: GNU Affero GPL v3
# The ruby implementation of the Haversine formula is free software: you can redistribute it and/or modify
# it under the terms of the GNU Affero General Public License version 3 as published by the Free Software Foundation.
#
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Affero General Public
# License version 3 for more details. http://www.gnu.org/licenses/
#
# Landon Cox - 9/25/08
#
# Notes:
#
# translated into Ruby based on information contained in:
# http://mathforum.org/library/drmath/view/51879.html Doctors Rick and Peterson - 4/20/99
# http://www.movable-type.co.uk/scripts/latlong.html
# http://en.wikipedia.org/wiki/Haversine_formula
#
# This formula can compute accurate distances between two points given latitude and longitude, even for
# short distances.
# PI = 3.1415926535
RAD_PER_DEG = 0.017453293 # PI/180
# the great circle distance d will be in whatever units R is in
Rmiles = 3956 # radius of the great circle in miles
Rkm = 6371 # radius in kilometers...some algorithms use 6367
Rfeet = Rmiles * 5282 # radius in feet
Rmeters = Rkm * 1000 # radius in meters
@distances = Hash.new # this is global because if computing lots of track point distances, it didn't make
# sense to new a Hash each time over potentially 100's of thousands of points
=begin rdoc
given two lat/lon points, compute the distance between the two points using the haversine formula
the result will be a Hash of distances which are key'd by 'mi','km','ft', and 'm'
=end
def haversine_distance( lat1, lon1, lat2, lon2 )
dlon = lon2 - lon1
dlat = lat2 - lat1
dlon_rad = dlon * RAD_PER_DEG
dlat_rad = dlat * RAD_PER_DEG
lat1_rad = lat1 * RAD_PER_DEG
lon1_rad = lon1 * RAD_PER_DEG
lat2_rad = lat2 * RAD_PER_DEG
lon2_rad = lon2 * RAD_PER_DEG
# puts "dlon: #{dlon}, dlon_rad: #{dlon_rad}, dlat: #{dlat}, dlat_rad: #{dlat_rad}"
a = Math.sin(dlat_rad/2)**2 + Math.cos(lat1_rad) * Math.cos(lat2_rad) * Math.sin(dlon_rad/2)**2
c = 2 * Math.asin( Math.sqrt(a))
dMi = Rmiles * c # delta between the two points in miles
dKm = Rkm * c # delta in kilometers
dFeet = Rfeet * c # delta in feet
dMeters = Rmeters * c # delta in meters
@distances["mi"] = dMi
@distances["km"] = dKm
@distances["ft"] = dFeet
@distances["m"] = dMeters
end
def test_haversine
lon1 = -104.88544
lat1 = 39.06546
lon2 = -104.80
lat2 = lat1
haversine_distance( lat1, lon1, lat2, lon2 )
puts "the distance from #{lat1}, #{lon1} to #{lat2}, #{lon2} is:"
puts "#{@distances['mi']} mi"
puts "#{@distances['km']} km"
puts "#{@distances['ft']} ft"
puts "#{@distances['m']} m"
if @distances['km'].to_s =~ /7\.376*/
puts "Test: Success"
else
puts "Test: Failed"
end
end
test_haversine
var R = 6371; // km
var dLat = (lat2-lat1)*Math.PI/180;
var dLon = (lon2-lon1)*Math.PI/180;
var a = Math.sin(dLat/2) * Math.sin(dLat/2) +
Math.cos(lat1*Math.PI/180) * Math.cos(lat2*Math.PI/180) *
Math.sin(dLon/2) * Math.sin(dLon/2);
var c = 2 * Math.asin(Math.sqrt(a));
var d = R * c;
// adapted from C++ code on this page
+ (CGFloat)directMetersFromCoordinate:(CLLocationCoordinate2D)from toCoordinate:(CLLocationCoordinate2D)to {
static const double DEG_TO_RAD = 0.017453292519943295769236907684886;
static const double EARTH_RADIUS_IN_METERS = 6372797.560856;
double latitudeArc = (from.latitude - to.latitude) * DEG_TO_RAD;
double longitudeArc = (from.longitude - to.longitude) * DEG_TO_RAD;
double latitudeH = sin(latitudeArc * 0.5);
latitudeH *= latitudeH;
double lontitudeH = sin(longitudeArc * 0.5);
lontitudeH *= lontitudeH;
double tmp = cos(from.latitude*DEG_TO_RAD) * cos(to.latitude*DEG_TO_RAD);
return EARTH_RADIUS_IN_METERS * 2.0 * asin(sqrt(latitudeH + tmp*lontitudeH));
}
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);
my $lat1 = shift(@_);
my $lon1 = shift(@_);
my $lat2 = shift(@_);
my $lon2 = shift(@_);
# 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;
SELECT *, (
6371 * acos(cos(radians(search_lat)) * cos(radians(lat) ) *
cos(radians(lng) - radians(search_lng)) + sin(radians(search_lat)) * sin(radians(lat)))
) AS distance
FROM table
WHERE lat != search_lat AND lng != search_lng AND distance < 25
ORDER BY distance
FETCH 10 ONLY
#coding:UTF-8
"""
Python implementation of Haversine formula
Copyright (C) <2009> Bartek Górny <bartek@gorny.edu.pl>
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
"""
import math
def recalculate_coordinate(val, _as=None):
"""
Accepts a coordinate as a tuple (degree, minutes, seconds)
You can give only one of them (e.g. only minutes as a floating point number) and it will be duly
recalculated into degrees, minutes and seconds.
Return value can be specified as 'deg', 'min' or 'sec'; default return value is a proper coordinate tuple.
"""
deg, min, sec = val
# pass outstanding values from right to left
min = (min or 0) + int(sec) / 60
sec = sec % 60
deg = (deg or 0) + int(min) / 60
min = min % 60
# pass decimal part from left to right
dfrac, dint = math.modf(deg)
min = min + dfrac * 60
deg = dint
mfrac, mint = math.modf(min)
sec = sec + mfrac * 60
min = mint
if _as:
sec = sec + min * 60 + deg * 3600
if _as == 'sec': return sec
if _as == 'min': return sec / 60
if _as == 'deg': return sec / 3600
return deg, min, sec
def points2distance(start, end):
"""
Calculate distance (in kilometers) between two points given as (long, latt) pairs
based on Haversine formula (http://en.wikipedia.org/wiki/Haversine_formula).
Implementation inspired by JavaScript implementation from http://www.movable-type.co.uk/scripts/latlong.html
Accepts coordinates as tuples (deg, min, sec), but coordinates can be given in any form - e.g.
can specify only minutes:
(0, 3133.9333, 0)
is interpreted as
(52.0, 13.0, 55.998000000008687)
which, not accidentally, is the lattitude of Warsaw, Poland.
"""
start_long = math.radians(recalculate_coordinate(start[0], 'deg'))
start_latt = math.radians(recalculate_coordinate(start[1], 'deg'))
end_long = math.radians(recalculate_coordinate(end[0], 'deg'))
end_latt = math.radians(recalculate_coordinate(end[1], 'deg'))
d_latt = end_latt - start_latt
d_long = end_long - start_long
a = math.sin(d_latt/2)**2 + math.cos(start_latt) * math.cos(end_latt) * math.sin(d_long/2)**2
c = 2 * math.asin(math.sqrt(a))
return 6371 * c
if __name__ == '__main__':
warsaw = ((21, 0, 30), (52, 13, 56))
cracow = ((19, 56, 18), (50, 3, 41))
print points2distance(warsaw, cracow)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment