-
Haversine formula
/wiki/Haversine_formula
Calculate distance between two points
on a globe
/wiki/Calculate_distance_bet
ween_two_points_on_a_globe#Haversi
ne_Fo
rmula
Calculate distance
between two points on a globe
Contents
?
?
1 Haversine Formula
2 Implementations
o
2.1 APEX
o
2.2
Java
o
2.3 Objective C
o
2.4
Perl
o
2.5 PostgreSQL PL/PERLU
o
2.6
JavaScript
o
2.7 C++
o
2.8
Ruby
o
2.9 Python
o
2.10
MySQL
o
2.11 MSSQL
o
2.12
C#
o
2.13 Excel
o
2.14
PHP
o
2.15 Smalltalk
o
2.16
Tcl
o
2.17 NCL
Haversine Formula
The
haversine formula
is an
equation important in navigation, giving
great-circle distances between two
points on a sphere from their
longitudes and latitudes. It is a
special case of a more general formula
in spherical trigonometry, the law of
haversines, relating the sides
and
angles of spherical
Implementations
APEX
public
class LocatorUtil {
private
static Double EARTH_RADIUS = 6371.00; // Radius in
Kilometers
default
public
static
Double
calculateDistance(Double
lat1,
Double
lon1,
Double
lat2, Double lon2){
Double Radius = _RADIUS; //6371.00;
Double dLat = ans(lat2-lat1);
Double dLon = ans(lon2-lon1);
Double a = (dLat/2) * (dLat/2) +
(ans(lat1)) *
(ans(lat2))
*
(dLon/2) * (dLon/2);
Double c = 2 * ((a));
return Radius *
c;
}
public static
Double toRadians(Double degree){
//
Value degree * Pi/180
Double res =
degree * 3.1415926 / 180;
return res;
}
}
Java
view plaincopy to
clipboardprint?
1.
import
nt;
2.
3.
public
class
DistanceCalculator {
4.
5.
private
double
Radius;
6.
7.
// R = earth's radius (mean radius =
6,371km)
8.
// Constructor
9.
DistanceCalculator(
double
R)
{
10.
Radius = R;
11.
}
12.
13.
public
double
CalculationByDistance(GeoPoint StartP, GeoPoint
EndP) {
14.
double
lat1
= itudeE6()/1E6;
15.
double
lat2 = itudeE6()/1E6;
16.
double
lon1 =
gitudeE6()/1E6;
17.
double
lon2 =
gitudeE6()/1E6;
18.
double
dLat =
ans(lat2-lat1);
19.
double
dLon =
ans(lon2-lon1);
20.
double
a =
(dLat/
2
) *
(dLat/
2
) +
21.
(ans(lat1)) * (ans(lat2)) *
22.
(dLon/
2
) *
(dLon/
2
);
23.
double
c =
2
* ((a));
24.
return
Radius * c;
25.
}
26.
}
Objective C
// adapted from
C++ code on this page
+ (CGFloat)direct
MetersFromCoordinate:(CLLocationCoordinate2D)from
toCoordinate:(CLLocationCoordinate2D)to
{
static const double
DEG_TO_RAD = 0.943295769236907684886;
static const double
EARTH_RADIUS_IN_METERS = 6372797.560856;
double latitudeArc = (de -
de) * DEG_TO_RAD;
double longitudeArc =
(ude - ude) * DEG_TO_RAD;
double
latitudeH = sin(latitudeArc * 0.5);
latitudeH *= latitudeH;
double lontitudeH = sin(longitudeArc *
0.5);
lontitudeH *= lontitudeH;
double tmp = cos(de*DEG_TO_RAD) *
cos(de*DEG_TO_RAD);
return
EARTH_RADIUS_IN_METERS * 2.0 * asin(sqrt(latitudeH
+
tmp*lontitudeH));
}
Perl
# $$lat1 and $$lon1 are
the coordinates of the first point in radians
# $$lat2 and $$lon2 are the coordinates
of the second point in radians
my $$a =
sin(($$lat2 - $$lat1)/2.0);
my $$b =
sin(($$lon2 - $$lon1)/2.0);
my $$h =
($$a*$$a) + cos($$lat1) * cos($$lat2) * ($$b*$$b);
my $$theta = 2 * asin(sqrt($$h)); #
distance in radians
#
in
order
to
find
the
distance,
multiply
$$theta
by
the
radius
of
the
earth,
e.g.
# $$theta * 6,372.7976 =
distance in kilometres (value from
/wiki/Earth_radius)
PostgreSQL PL/PERLU
view
plaincopy to clipboardprint?
1.
create
or
replace
function
distance_in_km (
numeric
,
numeric
,
numeric
,
nu
meric
)
returns
numeric
as
2.
$$body$$
3.
use strict;
4.
use Math::Trig
qw(great_circle_distance deg2rad);
5.
6.
my $$lat1 =
shift(@_);
7.
my $$lon1 = shift(@_);
8.
my $$lat2 =
shift(@_);
9.
my $$lon2 = shift(@_);
10.
11.
# Notice the
90 - latitude: phi zero
is
at
the North Pole.
12.
sub NESW { deg2rad($$_[0]), deg2rad(90 -
$$_[1]) }
13.
my @L = NESW( $$lat1, $$lon1 );
14.
my @T = NESW( $$lat2, $$lon2 );
15.
my $$km = great_circle_distance(@L, @T,
6378);
16.
17.
return
$$km;
18.
$$body$$
19.
language plperlu strict security
definer;
Query that
returns nearest 10 locations within 25 km. Assumes
the
location is stored as two numeric
fields, lat and lng.
view plaincopy to
clipboardprint?
1.
SELECT
*, (
2.
6371 *
acos(cos(radians(search_lat)) * cos(radians(lat) )
*
3.
cos(radians(lng) -
radians(search_lng)) + sin(radians(search_lat)) *
sin
(radians(lat)))
4.
)
AS
distance
5.
FROM
table
6.
WHERE
lat != search_lat
AND
lng != search_lng
AND
distance < 25
7.
ORDER
BY
distance
8.
FETCH
10
ONLY
JavaScript
view plaincopy to
clipboardprint?
1.
var
R = 6371;
// km
2.
var
dLat = (lat2-lat1)*/180;
3.
var
dLon = (lon2-lon1)*/180;
4.
var
a = (dLat/2) * (dLat/2) +
5.
(lat1*/180) * (lat2*/180) *
6.
(dLon/2) * (dLon/2);
7.
var
c = 2 * ((a));
8.
var
d = R * c;
?
Original Source
C++
view plaincopy to
clipboardprint?
1.
/// @brief The
usual PI/180 constant
2.
static
const
double
DEG_TO_RAD =
0.943295769236907684886;
3.
/// @brief
Earth's quatratic mean radius for
WGS-84
4.
static
const
double
EARTH_RADIUS_IN_METERS = 6372797.560856;
5.
6.
/** @brief Computes the arc, in radian,
between two WGS-84 positions.
7.
*
8.
* The result is equal to Distan
ce(from,to)/EARTH_RADIUS_IN_METERS
od
e>
9.
* =
2*asin(sqrt(h(d/EARTH_RADIUS_IN_METERS
)))
10.
*
11.
* where:
12.
*
d
is the distance in meters between 'from' and 'to'
positions.
i>
13.
*
h
is the haversine function:
h(x)=sin?(x/2)
14.
*
15.
*
16.
* The haversine formula
gives:
17.
*
h(d/R) =
h()+h()+cos()*co
s()
18.
*
19.
* @sa
/wiki/Law_of_haversines
20.
*/
21.
double
ArcInRadians(
const
Position&
from,
const
Position& to) {
22.
double
latitudeArc = ( - )
* DEG_TO_RAD;
23.
double
longitudeArc = ( - )
* DEG_TO_RAD;
24.
double
latitudeH =
sin(latitudeArc * 0.5);
25.
latitudeH
*= latitudeH;
26.
double
lontitudeH =
sin(longitudeArc * 0.5);
27.
lontitudeH *= lontitudeH;
28.
double
tmp =
cos(*DEG_TO_RAD) * cos(*DEG_TO_RAD);
29.
return
2.0 *
asin(sqrt(latitudeH + tmp*lontitudeH));
30.
}
31.
32.
/** @brief Computes the distance, in
meters, between two WGS-84 positions.
33.
*
34.
* The result is equal to EARTH_
RADIUS_IN_METERS*ArcInRadians(from,to
)<
/code>
35.
*
36.
* @sa ArcInRadians
37.
*/
38.
double
DistanceInMeters(
const
Position& from,
const
Position& to) {
39.
return
EARTH_RADIUS_IN_METERS*ArcInRadians(from, to);
40.
}
?
Original Source
Ruby
view plaincopy to
clipboardprint?
1.
#
2.
#
3.
# haversine
formula to compute the great circle distance
between two points
given their latitude
and longitudes
4.
#
5.
#
Copyright (C) 2008, 360VL, Inc
6.
#
Copyright (C) 2008, Landon Cox
7.
#
8.
# (Landon
Cox)
9.
#
contact:
10.
#
/blog/businesscard/
11.
#
12.
#
LICENSE: GNU Affero GPL v3
13.
# The ruby
implementation of the Haversine formula is free
software: you can
redistribute it
and/or modify
14.
# it under
the terms of the GNU Affero General Public License
version 3 as p
ublished by the Free
Software Foundation.
15.
#
16.
#
This program is distributed in the hope that it
will be useful, but WITHOU
T ANY
WARRANTY; without even the
17.
# implied
warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.
See the GNU Affero
General Public
18.
# License
version 3 for more details. /licenses/
19.
#
20.
# Landon Cox
- 9/25/08
21.
#
22.
#
Notes:
23.
#
24.
#
translated into Ruby based on information
contained in:
25.
#
/library/drmath/view/ Doctors Rick and
Pe
terson - 4/20/99
26.
#
/scripts/
27.
#
/wiki/Haversine_formula
28.
#
29.
#
This formula can compute accurate distances
between two points given latit
ude and
longitude, even for
30.
# short
distances.
31.
32.
#
PI = 3.1415926535
33.
RAD_PER_DEG =
0.017453293
# PI/180
34.
35.
#
the great circle distance d will be in whatever
units R is in
36.
37.
Rmiles = 3956
# radius of the great circle in
miles
38.
Rkm = 6371
# radius in kilometers...some
algorithms use 6367
39.
Rfeet =
Rmiles * 5282
# radius in
feet
40.
Rmeters = Rkm
* 1000
# radius in meters
41.
42.
@distances
=
Hash.
new
# this
is global because if computing lots of track
point distances, it didn't
make
43.
# sense to new a Hash each time over
potentially 100
's of thousands of
points
44.
45.
=
begin
rdoc
46.
given two lat/lon points, compute the distance
between the two points usin
g the
haversine formula
47.
the result
will be a Hash of distances which are
key
'd by '
mi
','
km
','
ft
'
,
and '
m'
48.
=
end
49.
50.
def
haversine_distance(
lat1, lon1, lat2, lon2 )
51.
52.
dlon = lon2
- lon1
53.
dlat = lat2 - lat1
54.
55.
dlon_rad = dlon * RAD_PER_DEG
56.
dlat_rad =
dlat * RAD_PER_DEG
57.
58.
lat1_rad = lat1 * RAD_PER_DEG
59.
lon1_rad =
lon1 * RAD_PER_DEG
60.
61.
lat2_rad = lat2 * RAD_PER_DEG
62.
lon2_rad =
lon2 * RAD_PER_DEG
63.
64.
# puts
lat_rad}
65.
66.
a = (dlat_rad/2)**2 + (lat1_rad) * (lat2_rad) *
Ma
(dlon_rad/2)**2
67.
c = 2 * (
(a))
68.
69.
dMi =
Rmiles * c
# delta between the
two points in miles
70.
dKm = Rkm *
c
# delta in
kilometers
71.
dFeet =
Rfeet * c
# delta in
feet
72.
dMeters =
Rmeters * c
# delta in
meters
73.
74.
@distances
[
] = dMi
75.
@distances
[
] =
dKm
76.
@distances
[
] = dFeet
77.
@distances
[
] =
dMeters
78.
end
79.
80.
def
test_haversine
81.
82.
lon1 = -104.88544
83.
lat1 =
39.06546
84.
85.
lon2 =
-104.80
86.
lat2 = lat1
87.
88.
haversine_distance( lat1, lon1, lat2, lon2 )
89.
90.
puts
91.
puts
92.
puts
93.
puts
94.
puts
95.
96.
if
@distances
[
'km'
].to_s =~
/7.376*/
97.
puts
98.
else
99.
puts
100.
end
101.
102.
end
103.
104.
test_haversine
?
Original
Source
Python
view plaincopy to
clipboardprint?
1.
#coding:UTF-8
2.
3.
Python implementation of Haversine
formula
4.
Copyright
(C) <2009> Bartek G??rny
5.
6.
This program is free software: you can
redistribute it and/or modify
7.
it under the terms of the GNU General Public
License as published by
8.
the Free
Software Foundation, either version 3 of the
License, or
9.
(at your
option) any later version.
10.
11.
This
program is distributed in the hope that it will be
useful,
12.
but WITHOUT
ANY WARRANTY; without even the implied warranty
of
13.
MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the
14.
GNU General
Public License for more details.
15.
16.
You should have received a copy of the GNU General
Public License
17.
along with
this program. If not, see
.
18.
19.
20.
import
math
21.
22.
def
recalculate_coordinate(val, _as=None):
23.
24.
Accepts a
coordinate as a tuple (degree, minutes,
seconds)
25.
You can
give only one of them (e.g. only minutes as a
floating point num
ber) and it will be
duly
26.
recalculated into degrees, minutes and
seconds.
27.
Return
value can be specified as 'deg', 'min' or 'sec';
default return v
alue is a proper
coordinate tuple.
28.
29.
deg, min, sec = val
30.
# pass outstanding values from right to
left
31.
min = (min
or
0) + int(sec) / 60
32.
sec = sec % 60
33.
deg = (deg
or
0) + int(min) / 60
34.
min = min % 60
35.
# pass decimal part from left to
right
36.
dfrac,
dint = (deg)
37.
min = min +
dfrac * 60
38.
deg = dint
39.
mfrac,
mint = (min)
40.
sec = sec +
mfrac * 60
41.
min = mint
42.
if
_as:
43.
sec = sec
+ min * 60 + deg * 3600
44.
if
_as ==
'sec'
:
return
sec
45.
if
_as ==
'min'
:
return
sec / 60
46.
if
_as ==
'deg'
:
return
sec / 3600
47.
return
deg, min, sec
48.
49.
50.
def
points2distance(start,
end):
51.
52.
Calculate
distance (in kilometers) between two points given
as (long, la
tt) pairs
53.
based on Haversine formula
(/wiki/Haversine_formu
la).
54.
Implementation inspired by JavaScript
implementation from
/scripts/
55.
Accepts
coordinates as tuples (deg, min, sec), but
coordinates can be gi
ven in any form -
e.g.
56.
can
specify only minutes:
57.
(0,
3133.9333, 0)
58.
is
interpreted as
59.
(52.0,
13.0, 55.998)
60.
which,
not accidentally, is the lattitude of Warsaw,
Poland.
61.
62.
start_long = s(recalculate_coordinate(start[0],
'deg'
))
63.
start_latt
= s(recalculate_coordinate(start[1],
'deg'
))
64.
end_long =
s(recalculate_coordinate(end[0],
'deg'
))
65.
end_latt =
s(recalculate_coordinate(end[1],
'deg'
))
66.
d_latt =
end_latt - start_latt
67.
d_long =
end_long - start_long
68.
a =
(d_latt/2)**2 + (start_latt) * (end_latt) *
ma
(d_long/2)**2
69.
c = 2 *
((a))
70.
return
6371 * c
71.
72.
73.
if
__name__ ==
'__main__'
:
74.
warsaw =
((21, 0, 30), (52, 13, 56))
75.
cracow =
((19, 56, 18), (50, 3, 41))
76.
print
points2distance(warsaw, cracow)
?
Original
Source
MySQL
view
plaincopy to clipboardprint?
1.
DELIMITER //
2.
CREATE
DEFINER=`someone`@`somehost`
FUNCTION
`GetDistance`(
3.
lat1
numeric
(9,6),
4.
lon1
numeric
(9,6),
5.
lat2
numeric
(9,6),
6.
lon2
numeric
(9,6)
7.
)
RETURNS
decimal
(10,5)
8.
READS SQL
DATA
9.
BEGIN
10.
DECLARE
x
decimal
(20,10);
11.
DECLARE
pi
decimal
(21,20);
12.
SET
pi = 3.979323846;
13.
SET
x = sin( lat1 * pi/180
) * sin( lat2 * pi/180 ) + cos(
14.
lat1 *pi/180
) * cos( lat2 * pi/180 ) * cos(
abs
( (lon2 * pi/180) -
15.
(lon1 *pi/180) ) );
16.
SET
x = acos( x );
17.
RETURN
( 1.852 * 60.0 *
((x/pi)*180) ) / 1.609344;
18.
END
//
view plaincopy to
clipboardprint?
1.
Calculate the
distance
between
2 points
with
latitude & longitude
in
MySQL:
2.
3.
DISTANCE
in
Kilometres:
4.
5.
6371 * acos( cos( radians(LATITUDE1) )
* cos( radians( LATITUDE2 ) ) * cos(
radians( LONGITUDE1 ) -
radians(LONGITUDE2) ) + sin( radians(LATITUDE1) )
*
sin( radians( LATITUDE2 ) ) )
AS
DISTANCE
6.
7.
DISTANCE
in
Miles:
8.
9.
3956 * acos(
cos( radians(LATITUDE1) ) * cos( radians(
LATITUDE2 ) ) * cos(
radians(
LONGITUDE1 ) - radians(LONGITUDE2) ) + sin(
radians(LATITUDE1) ) *
sin( radians(
LATITUDE2 ) ) )
AS
DISTANCE
List all points
in table having distance between a designated
point (we
use a random point -
lat:45.20327, long:23.7806) less than 50 KM,
with latitude & longitude, in MySQL
(the table fields are coord_lat and
coord_long):
view plaincopy
to clipboardprint?
1.
List
all
having
DISTANCE<50,
in
Kilometres (considered
Earth radius 6371 KM)
:
2.
3.
SELECT
denumire, (6371 *
acos( cos( radians(45.20327) ) * cos( radians(
coor
d_lat ) ) * cos( radians( 23.7806 )
- radians(coord_long) ) + sin(
radians(4
5.20327) ) * sin(
radians(coord_lat) ) ))
AS
distanta
FROM
obiective
WHERE
c
oord_lat<>
''
AND
coord_long<>
''
having
distanta<50
ORDER
BY
distanta
desc
The above example was tested in MySQL
5.0.95 and 5.5.16 (Linux).
MSSQL
view plaincopy to
clipboardprint?
1.
CREATE
FUNCTION
[dbo].[GetDistance]
2.
(
3.
@lat1
Float
(8),
4.
@long1
Float
(8),
5.
@lat2
Float
(8),
6.
@long2
Float
(8)
7.
)
8.
RETURNS
Float
(8)
9.
AS
10.
BEGIN
11.
DECLARE
@R
Float
(8);
12.
DECLARE
@dLat
Float
(8);
13.
DECLARE
@dLon
Float
(8);
14.
DECLARE
@a
Float
(8);
15.
DECLARE
@c
Float
(8);
16.
DECLARE
@d
Float
(8);
17.
18.
SET
@R = 3960;
19.
SET
@dLat = RADIANS(@lat2 -
@lat1);
20.
SET
@dLon =
RADIANS(@long2 - @long1);
21.
22.
SET
@a = SIN(@dLat / 2) *
SIN(@dLat / 2) + COS(RADIANS(@lat1))
23.
* COS(RADIANS(@lat2)) * SIN(@dLon / 2) *SIN(@dLon
/
2);
24.
SET
@c = 2 *
ASIN(
MIN
(SQRT(@a)));
25.
SET
@d = @R * @c;
26.
RETURN
@d;
27.
END
28.
GO
Pentultimate line of
formula is junk. Is: SET @c = 2 *
ASIN(MIN(SQRT(@a))); Should be: SET @c
= 2 * atn2(sqrt(@a),
sqrt(1-@a)); Not
sure what R is here, use 6371 for approx Km
differences.
Not sure about
that. I believe this returns Statute Miles based
on the
Radius @R value
C#
view plaincopy to
clipboardprint?
1.
using
System;
2.
namespace
HaversineFormula
3.
{
4.
///
5.
/// The distance type to return the
results in.
6.
///
7.
public
enum
DistanceType { Miles,
Kilometers };
8.
///
9.
/// Specifies a Latitude / Longitude
point.
10.
///
11.
public
struct
Position
12.
{
13.
public
double
Latitude;
14.
public
double
Longitude;
15.
}
16.
class
Haversine
17.
{
18.
///
19.
/// Returns the distance in miles or
kilometers of any two
20.
/// latitude / longitude
points.
21.
///
22.
///
name=”pos1
″
>
23.
///
name=”pos2
″
>
24.
///
name=”type”>
25.
///
26.
public
double
Distance(Position
pos1, Position pos2, DistanceType ty
pe)
27.
{
28.
double
R = (type == ) ? 3960
: 6371;
29.
double
dLat =
this
.toRadian(de - de);
30.
double
dLon =
this
.toRadian(ude - ude);
31.
double
a = (dLat / 2) *
(dLat / 2) +
32.
(
this
.toRadian(de)) * (
this
.toRad
ian(de)) *
33.
(dLon / 2) * (dLon / 2);
34.
double
c = 2 * ((1, (a)));
35.
double
d = R * c;
36.
return
d;
37.
}
38.
///
39.
/// Convert to Radians.
40.
///
41.
///
42.
///
43.
private
double
toRadian(
double
val)
44.
{
45.
return
( / 180) * val;
46.
}
47.
}
48.
}
To use this code you will
need 2 latitude/longitude points. I created a
struct to hold the latitude/longitude
points. Once we have the points,
we
create the Haversine class and call the Distance
method, passing
in the points and an
enum specifying whether to return the results in
Kilometers or Miles. We end up with the
following:
view plaincopy to
clipboardprint?
1.
Position pos1
=
new
Position();
2.
de
= 40.7486;
3.
ude = -73.9864;
4.
Position pos2
=
new
Position();
5.
de
= 24.7486;
6.
ude = -72.9864;
7.
Haversine hv =
new
Haversine();
8.
double
result = ce(pos1,
pos2, ters);
?
Original Source
Excel
view plaincopy to
clipboardprint?
1.
Public
Function
getDistance(latitude1, longitude1, latitude2,
longitude2)
2.
earth_radius = 6371
3.
Pi =
3.14159265
4.
deg2rad = Pi / 180
5.
6.
dLat = deg2rad
* (latitude2 - latitude1)
7.
dLon = deg2rad
* (longitude2 - longitude1)
8.
9.
a = Sin(dLat /
2) * Sin(dLat / 2) + Cos(deg2rad * latitude1) *
Cos(deg2rad *
latitude2) * Sin(dLon /
2) * Sin(dLon / 2)
10.
c = 2 *
(Sqr(a))
11.
12.
d = earth_radius * c
13.
14.
getDistance =
d
15.
16.
End
Function
PHP
view plaincopy to
clipboardprint?
1.
function
getDistance(
$$latitude1
,
$$longitude1
,
$$latitude2
,
$$longitude2
) {
2.
$$earth_radius
= 6371;
3.
4.
$$dLat
=
deg2rad(
$$latitude2
-
$$latitude1
);
5.
$$dLon
=
deg2rad(
$$longitude2
-
$$longitude1
);
6.
7.
$$a
=
sin(
$$dLat
/2) *
sin(
$$dLat
/2) +
cos(deg2rad(
$$latitude1
)) * c
os(deg2ra
d(
$$latitude2
)) * sin(
$$dLon
/2) *
sin(
$$dLon
/2);
8.
$$c
= 2 *
asin(sqrt(
$$a
));
9.
$$d
=
$$earth_radius
*
$$c
;
10.
11.
return
$$d
;
12.
}
Smalltalk
lat1Rad := latitude * Float pi / 180.
lon1Rad := longitude * Float pi / 180.
lat2Rad := latitude2 * Float pi / 180.
lon2Rad := longitude2 * Float pi / 180.
earthRadius := 6371.00.
dLat := lat2Rad
- lat1Rad.
dLon := lon2Rad - lon1Rad.
dLatSinSqrd := (dLat / 2) sin squared.
dLonSinSqrd := (dLon / 2) sin squared.
cosLatLat := lat2Rad cos * lat1Rad cos.
a := dLatSinSqrd + (cosLatLat *
dLonSinSqrd).
c := 2 * a sqrt arcSin.
distance := earthRadius * c.
^distance
Tcl
view plaincopy to
clipboardprint?
1.
proc deg2rad
deg {expr
{
$$deg
*atan(1)*8/360}}
2.
3.
proc dist {lat1 lat2 lon1 lon2} {
4.
set R 6371
5.
set dLat [deg2rad [expr
{(
$$lat2
-
$$lat1
)}]]
6.
set dLon
[deg2rad [expr {(
$$lon2
-
$$lon1
)}]]
7.
set sdlat2
[expr {sin(
$$dLat
/2)}]
8.
set sdlon2 [expr
{sin(
$$dLon
/2)}]
9.
set a [expr {
$$sdlat2
*
$$sdlat2
+ cos([deg2rad
$$lat1
])*cos([deg2rad
$$lat2
])*
$$s
dlo
n2
*
$$sdlon2
}]
10.
set d [expr {2*
$$R
*asin(sqrt(
$$a
))}]
11.
}
NCL
view
plaincopy to clipboardprint?
1.
undef
(
)
2.
function
haversine_dist(
3.
lat1[*]:
numeric,
4.
lon1[*]: numeric,
5.
lat2[*]:
numeric,
6.
lon2[*]: numeric
7.
)
8.
local a, c, delPhi, delPsi, R, pi,
deg2rad
9.
begin
10.
; Earth
radius
in
km,
and
other constants
11.
R = 6371.
12.
pi = 4.*atan(1.)
13.
deg2rad =
pi/180.
14.
; check that all arrays are of the same
length
15.
if
(
16.
dimsizes(lat1) .ne. dimsizes(lat2)
.
or
.
17.
dimsizes(lon1) .ne. dimsizes(lon2)
.
or
.
18.
dimsizes(lat1) .ne. dimsizes(lon1)
19.
)
then
20.
print(
ing
)
21.
exit
22.
end
if
23.
; Implement
Haversine distance formula
24.
; see /wiki/C
alculate_Distance_Between_Two_Points_on
_a_Globe
25.
delPhi = (lat1 - lat2) * deg2rad
26.
delPsi = (lon1 - lon2) * deg2rad
27.
a =
sin(delPhi/2.)*sin(delPhi/2.) +
28.
cos(lat1*
deg2rad)*cos(lat2*deg2rad)*sin(delPsi/2.)*sin(delP
si/2.)
29.
c = 2.*atan2(sqrt(a), sqrt(1-a))
30.
dist = R * c
31.
return
(dist)
32.
end
Categories
:
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
?
Geometry
Apex
Navigation
Perl
MSSQL
Java
JavaScript
C++
Ruby
Objective
C
Python
Tcl
C
sharp
Excel
Smalltalk
Objective-C
Arc
and
distance
between
two
points
on
Earth
surface
/2008/10/
Distance using Longitiude and latitude
using c++
/Articles/22488/Distance-using-
Longitiude-and-latitude-using-c
Movable Type
Scripts
Calculate
distance,
bearing
and
more
between Latitude/Longitude
points
This page presents a
variety of calculations for latitude/longitude
points, with the formul?
and code
fragments for implementing them.
All
these formul? are for calculations on the basis of
a spherical earth (ignoring ellipsoidal
effects)
–
which
is accurate enough
*
for most purposes… [In fact, the earth
is very slightly
ellipsoidal; using a
spherical model gives errors typically up to 0.3%
–
see notes for further
details].
窗体顶端
Great-circle distance between two
points
Enter the co-
ordinates into the text boxes to try out the
calculations. A variety of formats are
accepted, principally:
deg-min-sec suffixed with N/S/E/W (e.g.
40°44′55″N, 73 59 11W), or
signed decimal degrees
without compass direction, where negative
indicates
west/south (e.g. 40.7486,
-73.9864):
50 03 59N
Point 1:
,
005 42 53W
58 38 38N
Point 2:
Distance:
,
003 04 12W
968.9
km
(to 4 SF
)
*
Initial bearing:
009°07′11″
Final
bearing:
011°16′31″
Midpoint:
54°21′44″N,
004°31′50″W
And you can see
it on a map (aren‘t those Google guys
wonderful!)
... hide map
窗体底端
Distance
This uses the
?
haversine
‘ formula to
calculate the great
-circle distance
between two points
–
that is, the shortest distance
over the earth‘s surface –
giving an ?as
-the-
crow-
flies‘ distance
between
the points (ignoring any hills, of course!).
a = sin?(Δφ/2) + cos(φ
1
).cos(φ
2
).sin?(Δλ/2)
Haversine
c =
2(√
a
, √(1?a)
)
formula:
d = R.c
where
φ is
latitude, λ is longitude, R is
earth’s radius (mean radius =
6,371km)
note that angles
need to be in radians to pass to trig
functions!
var
R
=6371;// km
var
dLat
=(
lat2
-
lat1
).<
/p>
toRad
();
var
dLon
=(
lon2
-
lon1
).<
/p>
toRad
();
var
lat1
=
lat1
.
toRad
();
var
lat2
=
lat2
.
toRad
();
var
a
=Math.
sin
(
dLat
/2)
*Math.
sin
(
dLat
/2)+
Math.
sin
(
dLon
/2)*Math
.
sin
(
dLon
< br>/2)*Math.
cos
(
l
at1
)*Math.
cos
(
lat2
);
var
c
=2*Math.
atan2
(Math.
sqrt
(
a
),Math.
sq
rt
(1-
a
));
< br>
JavaScript:
var
d
=
R
*
c
;
The
haversine
formula
1
?remains
particularly well
-conditioned for
numerical computation
even at small
distances‘ –
unlike calculations based
on the
spherical law of
cosines
. The
?versed
sine‘
is 1-
cosθ, and the
?half
-versed-
sine‘
(1
-
cosθ)/2 = sin?(θ/2) as
used above. It
was published by Roger
Sinnott in
Sky & Telescope
magazine in 1984 (―Virtues of the
Haversine‖), though known about for
much longer by navigators. (For the curious,
c
is the
angular
distance in radians, and
a
is the square of half the chord length between the
points).
A (surprisingly marginal)
performance improvement can be obtained, of
course, by factoring
out the terms
which get squared.
Spherical Law of
Cosines
In fact, when Sinnott published
the haversine formula, computational precision was
limited.
Nowadays, JavaScript (and most
modern computers & languages) use IEEE 754 64-bit
floating-point numbers, which provide
15 significant figures of precision. With this
precision,
the simple
spherical law of cosines
formula (
cos c = cos a cos b + sin a
sin b cos C
) gives
well-
conditioned results down to distances as small as
around 1 metre.
(
Note that
the geodetic
form of the law of cosines
is rearranged from the canonical one so that the
latitude can be used directly,
rather
than the
colatitude
).
This makes the simpler law of cosines a
reasonable 1-line alternative to the haversine
formula for many purposes. The choice
may be driven by coding context, available trig
functions (in different languages),
etc.
Spherical
law of
cosines:
d = acos( sin(φ
1
).sin(φ
2
) + cos(φ
1
).cos(φ
2
).cos(Δλ) )
.R
var
R
=6371;// km
var
d
=Math.
acos
(Math.
sin
(
lat1
)*Math.
sin
(
lat2
)+
Math.
cos
(
lat1
)*Math.
cos
(
lat2
)*
JavaScript:
Excel:
Math.
cos
(
lon2
-
lon1
))*
R
;
=ACOS(SIN(lat1
)*SIN(lat2)+COS(lat1)*COS(lat2)*COS(lon2-lon1))*63
71
(
Note that here and in
all subsequent code fragments, for simplicity I do
not show conversions from
degrees to
radians; see code below for complete
versions
).
Equirectangular
approximation
If performance is an
issue and accuracy less important, for small
distances
Pythagoras‘
theorem
can be used on an
equirectangular
projection
:
*
Formula
x =
Δλ.cos(φ)
y = Δφ
d =
R
.√
x?
+
y?
var
x
< br>=(
lon2
-
lon1
)*Math.
cos
((
< br>lat1
+
lat2
)/2);
var
y
=(
lat2
-
lat1
);
var
d
=Math.
sqrt
(
< br>x
*
x
+
y
*
y<
/p>
)
*
R
;
JavaScript:
(lat/lon in radians!)
This uses just one trig and one sqrt
function
–
as against half-
a-dozen trig functions for cos
law, and
7 trigs + 2 sqrts for haversine. Accuracy is
somewhat complex: along meridians
there
are no errors, otherwise they depend on distance,
bearing, and latitude, but are small
enough for many
purposes
*
(and often trivial
compared with the spherical approximation
itself).
Bearing
Baghdad to Osaka
–
not
a constant bearing!
In general, your
current heading will vary as you follow a great
circle path (orthodrome); the
final
heading will differ from the initial heading by
varying degrees according to distance and
latitude (if you were to go from say
35°
N,45°
E (Baghdad) to
35°
N,135°
E (Osaka), you
would
start on a heading of
60°
and end up on a heading of
120°
!).
This formula is for
the initial bearing (sometimes referred to as
forward azimuth) which if
followed in a
straight line along a great-circle arc will take
you from the start point to the end
point:
1
Formula:
θ
= atan2(
sin(Δλ).cos(φ
2
), cos(φ
1
).sin(φ
2
) ?
sin(φ
1
).cos(φ
2
).cos(Δλ) )
var
y
=Math.
sin
(
dLon
)*M
ath.
cos
(
lat2
);
var
x
=Math.
cos
(
lat1
)*Math.
sin
(
lat2
)-
Math.
sin
(
lat1
)*Math.
cos
(
lat2
)*Math.
cos
(
dLon
);
JavaScript:
var
brng
=Math.
atan2
(
y
,
x
).
toDeg
();
< br>=ATAN2(COS(lat1)*SIN(lat2)-SIN(lat1)*COS(lat2)* COS(lon2-lon1),
SIN(lon2-lon1)*COS(lat2))
Excel:
*note that Excel reverses the arguments
to ATAN2
–
see notes
below
Since atan2 returns
values in the range -
π ... +π (that is,
-180°
... +180°
),
to normalise the result to a
compass
bearing (in the range 0°
...
360°
, with
?ve values
transformed into the range 180° ... 360°),
convert to degrees and then use
(θ+360)
% 360, where % is modulo.
For final bearing, simply take the
initial
bearing from the
end
point to the
start
point and reverse it
(using
θ = (θ+180) % 360).
Midpoint
This is the half-
way point along a great circle path between the
two points.
1
Formula:
B
x
=
cos(φ
2
).cos(Δλ)
B
y
=
cos(φ
2
).sin(Δλ)
φ
m
=
atan2( sin(φ
1
) +
sin(φ
2
), √((cos(φ
1
)+B
x
)?
+ B
y
?
) )
λ
m
=
λ
1
+
atan2(B
y
, cos(φ
1<
/p>
)+B
x
)
v
arBx=Math.
cos
(
lat
2
)*Math.
cos
(
dLon
);
varB
y=Math.
cos
(
lat2
)*Math.
sin
(
< br>dLon
);
var
lat3
=Math.
atan2
(Math.
sin
(
lat1
)+Math.
sin
(
lat2
),
Math.
sqrt
((
Math.
cos
(
lat1
)+Bx)*(Math.
cos
(
lat1
)+Bx)+By*By));
JavaScript:
var
lon3
=
lon1
+Math
.
atan2
(By,Math.
co
s
(
lat1
)+Bx);
Just as the initial bearing may
vary from the final bearing, the midpoint may not
be located half-way
between
latitudes/longitudes; the midpoint between
35°
N,45°
E and
35°
N,135°
E is around
45°
N,90°
E.
Destination point given
distance and bearing from start point
Given a start point, initial bearing,
and distance, this will calculate the destination
point and
final bearing travelling
along a (shortest distance) great circle arc.
窗体顶端
Destination
point along great-circle given distance and
bearing from start point
53
°19′14″N
Start
point:
,
001
°43′47″W
096
°01′18″
Bearing:
Distance:
km
Destination point:
53°11′18″N, 000°08′00″E
Final bearing:
view map
hide map
097°30′52″
窗体底端
Formula:
φ
2
=
asin( sin(φ
1
)*cos(d/R) +
cos(φ
1
)*sin(d/R)*cos(θ)
)
λ
2
= λ
1
+
atan2(
sin(θ)*sin(d/R)*cos(φ
1
), cos
(d/R)?sin(φ
1
)*sin(φ
2
) )
φ is latitude, λ is
longitude,
θ
is the bearing
(in radians, clockwise from north),
d
is the distance
travelled,
R
is
the earth’s radius (
d/R
is
the angular distance, in
where
radians)
var
lat2
=Math.
asin
(Math.
sin
(
lat1
)*Math.
cos
(
d
/
R
)+
Math.
< br>cos
(
lat1
)*Math
.
sin
(
d
/
R
)*Math.
cos
(
brng
));
var
lon2
=
lon1
+Math.
atan2
(Math.
sin
(
brng
)*Math.
s
in
(
d
/
R
)*Math.
cos
(
lat1
),
JavaScript:
Math.
cos
(
d
/
R
)-Math.
sin
(
lat1
)*Math.
sin
< br>(
lat2
));
lat2: =ASIN(SIN(lat1)*COS(d/R) +
COS(lat1)*SIN(d/R)*COS(brng))
lon2:
=lon1
+
ATAN2(COS(d/R)-SIN(lat1)*SIN(lat2),
SIN(brng)*SIN(d/R)*COS(lat1))
Excel:
* Remember that Excel
reverses the arguments to ATAN2
–
see notes below
For final bearing, simply take the
initial
bearing from the
end
point to the
start
point and reverse it
(using
θ = (θ+180) % 360).
Intersection
of two paths given start points and bearings
This is a rather more complex
calculation than most others on this page, but
I've been asked
for it a number of
times. See below for the JavaScript.
窗体顶端
Intersection
of two great-circle paths
51.885 N
Point 1:
,
0.235 E
Brng 1:
108.63
°
49.008 N
Point 2:
,
2.549 E
Brng 2:
32.72
°
Intersection point:
50°54′06″N, 004°29′39″E
窗体底端
Formula:
d
12
=
( √(sin?(Δφ/2)
+
cos(φ
1
).cos(φ
2
).s
in?(Δλ/2))
)
-
-
-
-
-
-
-
-
-
上一篇：grasshopper中文版运算器名称对照翻译
下一篇：GH命令英文详解