''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'' [http://www.surveyreview.org]: 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 mile And 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! ---- Related pages: * [Canvas geometrical calculations] * [GIS] ---- [[ [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) * spherical (1' = 1nm): 32657 metres * WGS84 (GPS): 32788 metres 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: * spherical (pseudo-WGS84): 32716 metres 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 * atan2(sqrt($a),sqrt(1-$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 [http://www-groups.dcs.st-and.ac.uk/~history/Mathematicians/Mechain.html])). If that is good enough depends on the application. ---- Longitude and latitude can be expressed in a number of ways. [LongLat Converter] converts between the different formats. ---- [Category Geography]