Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

As far as I understand the literature there are several ways to convert a set of geoecentric cartesian coordinates to geodetic coordinates

What would be the most efficient way to convert a cartesian coordinates to geodetic coordinates - most efficient meaning fastest (and most direct) when being implemented on a computer?

share|improve this question
1  
Do you need to perform conversions for arbitrary geodetic heights or only at nominal ground surface (zero heights)? –  whuber Jun 28 '12 at 14:00
1  
Ground surface will do for now, but I'm also interested in arbitrary geodetic heights in the long run. –  oschrenk Jun 28 '12 at 14:15
add comment

1 Answer

Depends on your accuracy requirements. One iteration of Heiskanen is efficient, but it takes 3 iterations to match the accuracy of a single iteration of Bowring's 1985 method. Even Lin and Wang's algorithm cannot match the efficiency of Bowring, if the trigonometric optimizations included below are used. Therefore, for all around performance I would recommend Bowring's 1985 algorithm.

   -- Based on B. R. Bowring's 1985 algorithm (single iteration)
   -----------------------------------------------------------------------------
   function ECEF_to_Geo(ECEF  : Vector3;
                        Datum : Datum_Id_Type) return Geographic_Type is
      Spheroid : Spheroid_Type := EARTH.DATUMS.Spheroid(Datum);
      x    : Meters   renames ECEF(1);
      y    : Meters   renames ECEF(2);
      z    : Meters   renames ECEF(3);
      a    : Meters   renames Models(Spheroid).Semi_Major_Axis;
      b    : Meters   renames Models(Spheroid).Semi_Minor_Axis;
      e2   : Unitless renames Models(Spheroid).Eccentricity_Squared;
      eb2  : Unitless renames Models(Spheroid).Second_Eccen_Squared;
      a2   : Scalar := a**2;
      d    : Meters := eb2*b;
      Ome2 : Unitless := 1.0 - e2;
      p2, p, r, tu, tu2, su3, cu3, Phi, Lam, tp, cp, sp, h : Scalar;
      Lat : Radians;
      Lon : Radians;
      Alt : Meters;
      Geo : Geographic_Type;
   begin -- ECEF_to_Geo
      if ((x = 0.0) and (y = 0.0)) then
         r   := abs(z);
         Phi := Scalar'copy_sign(Half_Pi, z);
         Lam := 0.0;
         h   := r - b;
      elsif (z = 0.0) then
         Phi := 0.0;
         Lam := arctan(y, x);
         p   := sqrt(x**2 + y**2);
         h   := p - a;
      else
         p2  := x**2 + y**2;
         p   := sqrt(p2);
         r   := sqrt(p2 + z**2);
         tu  := b*z*(1.0 + d/r)/(a*p);
         tu2 := tu**2;
         cu3 := (1.0/sqrt(1.0 + tu2))**3;
         su3 := cu3*tu2*tu;
         tp  := (z + d*su3)/(p - e2*a*cu3);
         Phi := arctan(tp);
         Lam := arctan(y, x);
         cp  := 1.0/sqrt(1.0 + tp**2);
         sp  := cp*tp;
         h   := p*cp + z*sp - a*sqrt(1.0 - e2*sp**2);
      end if;
      Lat := Radians(Phi);
      Lon := Radians(Lam);
      Alt := h;
      Geo := Make_Geo(Lat, Lon, Alt, Datum);
      return Geo;
   end ECEF_to_Geo;
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.