Author Topic: Dead reckoning function... (geodesic measurements)  (Read 6120 times)

Planetary Researcher

  • Ra (Power Member)
  • *****
  • Posts: 134
Dead reckoning function... (geodesic measurements)
« on: March 05, 2005, 04:11:06 PM »
Does anyone out there have a mathematical function which will take:


initial latitude (or co-latitude)
initial longitude
direction of travel (azimuthal angle east from north)
distance of travel
radius of the sphere whose surface you're on


and output the new latitude and longitude you'd find yourself at? Needs to be accurate for small distances, relative to the radius of the sphere.

thare

  • GIS Support Team
  • Isis (Extreme Power Member)
  • *****
  • Posts: 1217
    • http://webgis.wr.usgs.gov
    • Email
Dead reckoning function... (geodesic measurements)
« Reply #1 on: March 05, 2005, 04:12:57 PM »
It appears you are searching for some of Vincenty's formulas. Searching the web, I found a version written in Python. See the next post. Othewise, I have done this done this within ArcView 3.x calling the ReturnPointAlongGeodesic function. In ArcView 3.x and its' Avenue programming language, you can return the resultant X, Y point basically in one call:

aEndPoint = aSpheroid.ReturnPointAlongGeodesic(aStartPoint, aDistance, anAzimuth)

This states, given defined Spheroid call the routine ReturnPointAlongGeodesic sending it the Start Point (X, Y), the Distance, and the Azimuth and return the End Point. Gotta love that easy to read syntax. To bad it is basically a dead language. My extension Geodetic Profiler has some sample code you can use which calls this routine:
http://webgis.wr.usgs.gov/pgis_discussion/viewtopic.php?t=49

This can also be done in ArcMap using the geodcoord function. Search for that routine at http://arcobjectsonline.esri.com/ for some sample code that calls that routine.

There is also a downloadable Excel spreadsheet and online forms to calculate Vincenty Direct and Vincenty Inverse formulas:
http://www.ga.gov.au/nmd/geodesy/datums/calcs.jsp

Trent

thare

  • GIS Support Team
  • Isis (Extreme Power Member)
  • *****
  • Posts: 1217
    • http://webgis.wr.usgs.gov
    • Email
Dead reckoning function... (geodesic measurements)
« Reply #2 on: March 05, 2005, 04:13:46 PM »
The original python code for Vincenty's direct formula was found at (comments removed to fit here):
http://wegener.mechanik.tu-darmstadt.de/GMT-Help/Archiv/att-8710/Geodetic_py
Code: [Select]
# Vincenty's Direct formulae
def  vinc_pt( f, a, phi1, lembda1, alpha12, s ) :
   """
   Returns the lat and long of projected point and reverse azimuth
   given a reference point and a distance and azimuth to project.
   lats, longs and azimuths are passed in decimal degrees
   Returns ( phi2,  lambda2,  alpha21 ) as a tuple
   """
   piD4 = math.atan( 1.0 )
   two_pi = piD4 * 8.0
   phi1    = phi1    * piD4 / 45.0
   lembda1 = lembda1 * piD4 / 45.0
   alpha12 = alpha12 * piD4 / 45.0
   if ( alpha12 < 0.0 ) :
      alpha12 = alpha12 + two_pi
   if ( alpha12 > two_pi ) :
      alpha12 = alpha12 - two_pi

   b = a * (1.0 - f)
   TanU1 = (1-f) * math.tan(phi1)
   U1 = math.atan( TanU1 )
   sigma1 = math.atan2( TanU1, math.cos(alpha12) )
   Sinalpha = math.cos(U1) * math.sin(alpha12)
   cosalpha_sq = 1.0 - Sinalpha * Sinalpha
   
   u2 = cosalpha_sq * (a * a - b * b ) / (b * b)
   A = 1.0 + (u2 / 16384) * (4096 + u2 * (-768 + u2 * \
      (320 - 175 * u2) ) )
   B = (u2 / 1024) * (256 + u2 * (-128 + u2 * (74 - 47 * u2) ) )
   
   # Starting with the approx
   sigma = (s / (b * A))
   last_sigma = 2.0 * sigma + 2.0   # something impossible
   
   # Iterate the following 3 eqs unitl no sig change in sigma
   # two_sigma_m , delta_sigma
   while ( abs( (last_sigma - sigma) / sigma) > 1.0e-9 ) :

      two_sigma_m = 2 * sigma1 + sigma
      delta_sigma = B * math.sin(sigma) * ( math.cos(two_sigma_m) \
            + (B/4) * (math.cos(sigma) * \
            (-1 + 2 * math.pow( math.cos(two_sigma_m), 2 ) -  \
            (B/6) * math.cos(two_sigma_m) * \
            (-3 + 4 * math.pow(math.sin(sigma), 2 )) *  \
            (-3 + 4 * math.pow( math.cos (two_sigma_m), 2 ))))) \
     
      last_sigma = sigma
      sigma = (s / (b * A)) + delta_sigma
   
   phi2 = math.atan2 ( (math.sin(U1) * math.cos(sigma) + math.cos(U1) * math.sin(sigma) * math.cos(alpha12) ), \
      ((1-f) * math.sqrt( math.pow(Sinalpha, 2) +  \
      pow(math.sin(U1) * math.sin(sigma) - math.cos(U1) * math.cos(sigma) * math.cos(alpha12), 2))))
   
   lembda = math.atan2( (math.sin(sigma) * math.sin(alpha12 )), (math.cos(U1) * math.cos(sigma) -  \
           math.sin(U1) *  math.sin(sigma) * math.cos(alpha12)))
   
   C = (f/16) * cosalpha_sq * (4 + f * (4 - 3 * cosalpha_sq ))
   omega = lembda - (1-C) * f * Sinalpha *  \
      (sigma + C * math.sin(sigma) * (math.cos(two_sigma_m) + \
      C * math.cos(sigma) * (-1 + 2 * math.pow(math.cos(two_sigma_m),2) )))

   lembda2 = lembda1 + omega
   alpha21 = math.atan2 ( Sinalpha, (-math.sin(U1) * math.sin(sigma) +  \
      math.cos(U1) * math.cos(sigma) * math.cos(alpha12)))

   alpha21 = alpha21 + two_pi / 2.0
   if ( alpha21 < 0.0 ) :
      alpha21 = alpha21 + two_pi
   if ( alpha21 > two_pi ) :
      alpha21 = alpha21 - two_pi

   phi2       = phi2       * 45.0 / piD4
   lembda2    = lembda2    * 45.0 / piD4
   alpha21    = alpha21    * 45.0 / piD4
   return phi2,  lembda2,  alpha21
   # END