Version 5 of geodesy

Updated 2003-07-10 10:52:16

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 betwen two point on the earth giving 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 mile

And then a conversion factor between degress (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 -> of 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)}]  
 }

Thats it!


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

  • spherical: 63735 metres
  • WGS84 (GPS): 63993 metres

The difference is 258 metres. It depends on the application if this difference is significant enough.