Geodesy is the "science concerned with the study of the shape and size of the earth in the geometrical sense and with the study of certain physical phenomena, such as gravity, in seeking explanation of fine irregularities in the earth's shape." (D.H. Maling: Coordinate Systems and Map Projections, 2nd edition, Pergamon press 1992).
So what has Tcl to do with that? All geographic applications like GIS need geodesy for their calculations. One of the important calculations is the measurement of the distance between two points on the earth given the two geographical positions. This distance can be calculated using T. Vincenty's algorithm as published in Survey Review [L1 ]:
Vincenty T., (1975) "Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations", Survey Review 23, Number 176, p 88-93
The only thing we have to know beyond the two positions are the parameters of the ellipsoid in which the positions are expressed. This is important since the result will vary according to the used ellipsoid. So we define:
# all known Ellipsoids:
# each entry is a list with: ("major semi-axis","minor semi-axis","name")
# axes are in metres
set App(ellipsoid,wgs84) [list 6378137 6356750.321 "WGS84"]
set App(ellipsoid,grs80) [list 6378137 6356752.314 "GRS80"]
set App(ellipsoid,clarke) [list 6378206.4 6356583.8 "Clarke (1866) / NAD27"]
set App(ellipsoid,int) [list 6378388 6356911.946 "International (1909/1924)"]
set App(ellipsoid,krassovsky) [list 6378245 6356863.019 "Krassovsky (1940)"]
set App(ellipsoid,bessel) [list 6377397.155 6356078.963 "Bessel (1841)"]
set App(ellipsoid,wgs72) [list 6378135 6356750.52 "WGS72"]
set App(ellipsoid,wgs66) [list 6378145 6356759.769 "WGS66"]
set App(ellipsoid,airy) [list 6377563.396 6356256.909 "Airy 1830"]Often you will want to use the WGS84 or GRS80 ellipsoid since GPS data are expressed with reference to these. Now we need to define some constant values:
set App(PI) [expr {atan(1)*4}] ;# this is pi
set App(nm) 1852.0 ;# this is the number of meters that go into a nautical mileAnd then a conversion factor between degrees (ranging from 0 to 360) and radians (ranging from 0 to 2*pi):
set App(d2r) [expr {2*$App(PI)/360.0}]Now we are ready to go along with the algorithmic solution. It uses a loop to calculate the correct value iteratively until the result is good enough:
# method -> one of the elements of the above ellipsoid array, e.g., 'ellipsoid,bessel'
# x1 ... y2 -> position of the points 1 and 2 in decimal degrees
# returns the distance between the points expressed in metres
proc DistEllipsoid {method x1 y1 x2 y2} {
global App
# set our precision
set eps 0.00000000005
# convert to radians:
set lon1 [expr {$x1 * $App(d2r)}]
set lon2 [expr {$x2 * $App(d2r)}]
set lat1 [expr {$y1 * $App(d2r)}]
set lat2 [expr {$y2 * $App(d2r)}]
# get ellipsoid data:
foreach {a0 b0 name} $App($method) break
# flatness:
set flat [expr {($a0 - $b0) / $a0}]
set r [expr {1 - $flat}]
#
set tanu1 [expr {$r * tan($lat1)}]
set tanu2 [expr {$r * tan($lat2)}]
#
set x [expr {atan($tanu1)}]
set cosu1 [expr {cos($x)}]
set sinu1 [expr {sin($x)}]
#
set x [expr {atan($tanu2)}]
set cosu2 [expr {cos($x)}]
set sinu2 [expr {sin($x)}]
#
set omega [expr {$lon2 - $lon1}]
set lambda $omega
# loop untile the result is good enough:
while {1} {
set testlambda $lambda
set s2s [expr {pow($cosu2 * sin($lambda),2) + \
pow($cosu1 * $sinu2 - $sinu1 * $cosu2 * cos($lambda),2)}]
set ss [expr {sqrt($s2s)}]
set cs [expr {$sinu1 * $sinu2 + $cosu1 * $cosu2 * cos($lambda)}]
set tansigma [expr {$ss / $cs}]
set sinalpha [expr {$cosu1 * $cosu2 * sin($lambda) / $ss}]
set x [expr {asin($sinalpha)}]
set cosalpha [expr {cos($x)}]
set c2sm [expr {$cs - (2 * $sinu1 * $sinu2 / pow($cosalpha,2))}]
set c [expr {$flat / 16.0 * pow($cosalpha,2) * \
(4 + $flat * (4 - 3 * pow($cosalpha,2)))}]
set lambda [expr {$omega + (1 - $c) * \
$flat * $sinalpha * (asin($ss) + $c * $ss * ($c2sm + \
$c * $cs * (-1 + 2 * pow($c2sm,2))))}]
# result good enough?
if {abs($testlambda - $lambda) <= 0.00000000000005} {break}
}
set u2 [expr {pow($cosalpha,2) * (pow($a0,2) - pow($b0,2)) / pow($b0,2)}]
set a [expr {1 + ($u2 / 16384.0) * (4096 + $u2 * \
(-768 + $u2 * (320 - 175 * $u2)))}]
set b [expr {($u2 / 1024.0) * (256 + $u2 * (-128 + $u2 * (74 - 47 * $u2)))}]
set dsigma [expr {$b * $ss * ($c2sm + ($b / 4.0) * \
($cs * (-1 + 2 * pow($c2sm,2)) - ($b / 6.0) * $c2sm * \
(-3 + 4 * pow($ss,2)) * (-3 + 4 * pow($c2sm,2))))}]
return [expr {$b0 * $a * (asin($ss) - $dsigma)}]
}That's it!
KBK There are a large number of calculations for map projections available in the mapproj package in tcllib already - http://tmml.sourceforge.net/doc/tcllib/mapproj.html . If you're dealing with maps at the scale of a continent, or of most countries or US States, you won't notice the errors between the ellipsoid and the sphere. If you are dealing with larger scale maps, I'd certainly welcome code that fits a reference geoid into the framework! (In any case, even large scale maps are usually projected on a polyconic or transverse-Mercator projection, with geodetic latitude lifted from the geoid onto a sphere (and the map scale adjusted accordingly).
Related pages:
[ CL needs to make time to demonstrate how much simpler the calculation is for a sphere, and how little inaccuracy that introduces.] [Nice work above, though, Torsten.]
TR - you are right, calculation for a sphere is much simpler:
set d [expr {acos(sin($y1)*sin($y2)+cos($y1)*cos($y2)*cos($x1-$x2))}]
set distance_in_metres [expr {180*60.0/3.1415926535*$d*1852}]See Canvas geometrical calculations for details. Under the assumption that 1 nautical mile is the same as 1 Minute in latitude, I can give an example:
Two points: x1=10 , y1=54 , x1=10.5 , y2=54
(Btw: I live at 10°8'16''E 54°22'14''N in WGS84 coordinates)
The difference is 131 metres. It depends on the application if this difference is significant enough.
We can, however, come a little closer to a good result by using the sperical formula with the parameters of the WGS84 ellipsoid. The distance $d above is the arc distance. We just had to multiply by the earth's radius to get the distance in metres. So, if we multiply by 6378137 (the major semi-axis of the WGS84 ellipsoid, we get:
Now we are "only" 72 metres away from the correct distance. This may still be too much for finding and digging up a treasure chest in a desert given a pair of coordinates ...
DG -- Here's two flavors of a sherical one I did:
proc geographical_distance {lat1 long1 lat2 long2} {
global pi
# Radius of the earth. Assume a perfect sphere.
set r_smi 3963.1
set r_km 6378
# convert to radians
set y1 [expr {$lat1*($pi/180)}]
set x1 [expr {$long1*($pi/180)}]
set y2 [expr {$lat2*($pi/180)}]
set x2 [expr {$long2*($pi/180)}]
if {0} {
# found @ http://jan.ucc.nau.edu/~cvm/latlon_formula.html
set c [expr {acos(cos($y1)*cos($x1)*cos($y2)*cos($x2) + \
cos($y1)*sin($x1)*cos($y2)*sin($x2) + sin($y1)*sin($y2))}]
} else {
# The "Haversine Formula"
# found @ http://mathforum.org/library/drmath/view/51879.html
set dlon [expr {$x2 - $x1}]
set dlat [expr {$y2 - $y1}]
set a [expr {pow((sin($dlat/2)),2) + \
cos($y1)*cos($y2)*pow((sin($dlon/2)),2)}]
set c [expr {2 * asin(sqrt($a))}]
}
# statute_miles kilometers
return [list [expr {$r_smi * $c}] [expr {$r_km * $c}]]
}Is there anyway to modify the ellipsoid equation at the top to return arc distance instead of meters? It looks so optimized, I can't understand the innards at all.
Lars H: Isn't there a problem that "arc distance" isn't all that well-defined on an ellipsoid? The angle between the two vectors from the midpoint to the two given points is one possible definition (and not at all hard to compute), but the problem is that this angle corresponds to the route from one point to the other that lies in the plane determined by the two given points and the midpoint, and this is in general not the shortest route!
Of course, if you only want the distance "in degrees" you can just translate the length to an angle as for a circle or something. The numerically easiest conversion is probably that 1 meter = 9e-6 degrees, based on the original definition of the meter as 1/10000000 of the distance from the north pole to the equator (along the meridian of Paris (however that didn't come out quite right, see [L2 ])). If that is good enough depends on the application.
AMG: At 2010-04-15 08:21:13, [email protected] edited the Haversine Formula to change from "set c [expr {2 * atan2(sqrt($a),sqrt(1-$a))}]" to "set c [expr {2 * asin(sqrt($a))}]". I'm afraid I don't understand this change. Could someone please explain it? If you think it's wrong (I suspect it might be), and you have more confidence in your assessment than I do, go ahead and revert the change back to page version 21.
Lars H: atan2($y,$x), by definition, returns a solution theta to tan($theta) == $y/$x. asin($y) by definition returns a solution theta to sin($theta) == $y. For y being sqrt($a) and x being sqrt(1-$a) both formulae compute the same angle theta, since x^2 + y^2 = 1. I must say those pow's and sqrt's in the code look a bit odd, though: anyone with a classical training in numerical analysis would have used the trigonometric identities for cosine of the double angle to get rid of them. Or on second though, maybe not: doing it with the sqrt's and pow's avoids cancellation, which for short distances could be very noticable. Still, the following rewrite using hypot might be a better alternative
set c [expr {2*asin( hypot( sin($dlat/2), sqrt(cos($y1)*cos($y2)) * sin($dlon/2) ) )}]
AMG: Oops, when I saw the change I did a quick check for a random value of $a and got different results, hence my idea that the identity did not hold. Now I redid the test and got the same result. I must have made a typo. Sorry guys! :^)
Longitude and latitude can be expressed in a number of ways. LongLat Converter converts between the different formats.
dcd: For the mathematical sluts among us who still think the earth is locally pretty darn flat, have forgotten the difference between the Gilead and the Geodesy, and who need to get a useful calculation done in a short amount of time, may I present the following two alternate distance calculations, cribbed from Wikipedia. It's very informative to compare these calculations with the Haversine / Great Circle version above. For distances under a few hundred kms, the flat earth variant performs as well as the Haversine formula and the FCC algorithm differs in ways that don't seem to justify it's use if you're in a hurry. Anyway, here they are:
# assume all points are {lat lon} where lat and lon
# are expressed in decimal degrees
# spherical projection on a flat earth
# surprisingly accurate for distances under a few hundred kilometers
# and requires only one cosine and five multiplications.
# (no divisions if divide by 2 == a right shift in fixed point :-)
# The ellipsoidal FCC flat earth projection requires 5 cosines
# and nine multiplications
proc fedist {pt1 pt2} {
set dlat [expr {[lat $pt1] - [lat $pt2]}]
set dlon [expr {[lon $pt1] - [lon $pt2]}]
set dlon [expr {$dlon * cos([d2r [expr {([lat $pt1] + [lat $pt2])/2.0}]])}]
set dlat [deg2km $dlat]
set dlon [deg2km $dlon]
expr {sqrt($dlat * $dlat + $dlon * $dlon)}
}
# for details on this algorithm see FCC 47 CFR 73.208
# it's an ellipsoidal (Clarke 1866) projection onto a flat earth
# the FCC says not to use this for distances over 475km
# probably isn't as good in Australia as in N. Amer.
proc fccdist {pt1 pt2} {
set K11 111.13209
set K12 0.56605
set K13 0.0012
set K21 111.41513
set K22 0.09455
set K23 0.00012
set dlat [expr {[lat $pt1] - [lat $pt2]}]
set dlon [expr {[lon $pt1] - [lon $pt2]}]
set midpt [d2r [expr {([lat $pt1] + [lat $pt2])/2.0}]]
set mphi $midpt
set K1 $K11
set K2 [expr {$K21 * cos($mphi)}]; # cos(phi)
set mphi [expr {$mphi + $midpt}]
set K1 [expr {$K1 - $K12 * cos($mphi)}]; # cos(2phi)
set mphi [expr {$mphi + $midpt}]
set K2 [expr {$K2 - $K22 * cos($mphi)}]; # cos(3phi)
set mphi [expr {$mphi + $midpt}]
set K1 [expr {$K1 + $K13 * cos($mphi)}]; # cos(4phi)
set mphi [expr {$mphi + $midpt}]
set K2 [expr {$K2 + $K23 * cos($mphi)}]; # cos(5phi)
expr {sqrt(pow($K1*$dlat,2) + pow($K2*$dlon,2))}
}
# some helper procs
proc deg2km d {expr {$d * 111.19492664455873}}
package require math::constants
math::constants::constants {degtorad}
proc lat p {lindex $p 0}
proc lon p {lindex $p 1}
proc d2r d {expr {$d * $::degtorad}}