Later, you discover the ellipsoid you are using is not the ellipsoid used by the rest of the world. Most of the world uses WGS 84 (including GPS) and your coordinates are GPS based, so you need to use the right formulas to do the conversion.
There is a lot of literature on Internet about this, but as soon you start to read you discover a big problem. Conversion from longitude, latitude and altitude is easy and straightforward, for example in Octave code:
Note: I wrote the functions in Octave/Matlab code without advanced functions to make the task of porting to other languages easy. For example, many sqrt instructions can be written in a more clear way using Octave's norm function.
Surprisingly, WGS 84 has no inverse model: you can calculate Cartesian coordinates easily from WGS 84 latitude, longitude and altitude. But going from Cartesian coordinates to WGS 84 latitude, longitude and altitude is very tricky.
There are two methods to make the conversion. The first one involves making iterations, so conversion precision is function of the number of iterations and convergence criteria. I want a simple and fast conversion, so I discarded this option.
The other method is a simple algorithm:
Well, I just said there is no inverse model, and now I show you a simple algorithm to make the conversion. What is the catch? The catch is this algorithm only works for points situated at altitude = 0, or really close to 0. If you play with it, you will find out it works fine finding latitude and longitude, but fails with altitude.
The solution I found to solve this problem is:
- Get longitude, latitude and altitude with xyz2lla using desired x, y, z values.
- Discard altitude and get ground coordinates of calculated latitude and longitude using lla2xyz(lon,lat,0)
- Call this new ground Cartesian coordinates gx, gy, gz
- Get the module of vectors [x y z] and [gx gy gz].
- Altitude is the distance between points x,y,z and gx,gy,gz. Subtract length of vector [gx gy gz] to length of vector [x y z] to get altitude.
In Octave code can be implemented in this way:
To test this function, I wrote this small Octave script to create 100000 conversions. It generates points at random latitude, longitude and altitude, converts them to X-Y-Z values, and then, it uses the X-Y-Z values to calculate latitude, longitude and altitude using my modified xyz2lla function.. At the end, a small report is presented:
And the output of this little script in my system is:
As you can see, the RMS error in altitude is a little bit under 1.6 meters, not bad for such simple function.
Some interesting readings about this:
https://gist.github.com/klucar/1536194 the functions I ported to Octave.
I hope this helps to someone who need it.
Miguel A. Vallejo, EA4EOZ