Calculate lat/long/alt values for object with current lat/long/alt + movement on x, y, z

Question

I thought there would be plenty of results when searching for this on Google, but it turns out that in most cases, just like here on stackoverflow, the questions are very specific and include Google Maps, some GIS or the real world with elliptical shape.

In my case,

• I have an object with lat, long and alt values, where lat and long are floats in degrees and alt is a float as decimal value and representing the distance from the center of the sphere. So there are no minutes (10°30' is 10,5).
• The object moves on all 3 axes x, y and z, where the movement is relative to its current position
• I need to calculate new lat, long and alt values

Without thinking about it a lot, I first implemented it like this:

Just added z movement to alt, then calculated the virtual sphere's circumference to get the degree per unit (meter). With that I calculated the new lat first, and new long afterwards. I knew calculating one value after another is wrong, but as long as all the calculations in my "object world" would be done the same way, the overall behaviour would be ok. I put some thought into things like when the object goes around the entire sphere the long value doesn't change, and that going half around the sphere is different on x axis than on y axis (x axis: -180 to 180, y-axis -90 to 90) and things like that and it worked.

But then I realized that I didn't take into account that I calculated the degree per meter on the equator and didn't take other latitude values into account. At that point I knew things were more complicated and I started searching the web. But I didn't find an algorithm that would fit my needs. Now I'm sure this has been done before a lot of times, so I'm asking here in case someone has dealt with this topic before and can point me to a nice implementation :) .

What I found is this:

• Algorithms converting lat/long to Mercator projection
• Haversine formular for calculating a distance from two lat/long value pairs
• other formulars for other purposes

Which didn't help me.

What helped me the most was this: Calculate latitude and longitude having meters distance from another latitude/longitude point

But I think I didn't fully understand it yet.

• What is the course in radians? (regarding I have x and y movements)
• The extrapolate method doesn't take the altitude / earth's radius into account?

It would be awesome if someone could help me with this!

(Sidenote: I'm implementing this in Erlang, but that doesn't matter, any kind of algorithm would help)

---

UPDATE

I implemented the function mentioned above and the one from the answer below. When testing I got wrong values, maybe because I made mistakes in the implementation or calculated wrong testing data. Let's see:

Implementation 1:

% My own implementation that only works on the equator so far. Simple calculations as mentioned above.

Tests (on equator) ok.

Implementation 2:

% @doc calculates new lat+long+alt values for old values + movement vector
% https://stackoverflow.com/questions/5857523/calculate-latitude-and-longitude-having-meters-distance-from-another-latitude-lo
calc_position2(LastCurrLat, LastCurrLong, LastCurrAlt, MoveX, MoveY, MoveZ) ->
    % first the new altitude
    NewCurrAlt = LastCurrAlt + MoveZ,

    % original algorithm: http://williams.best.vwh.net/avform.htm#LL
    % lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
    % dlon=atan2(sin(tc)*sin(d)*cos(lat1),cos(d)-sin(lat1)*sin(lat))
    % lon=mod(lon1+dlon +pi,2*pi)-pi
    % where:
    % lat1, lon1    - start point in radians
    % d             - distance in radians
    % tc            - course in radians

    % -> for the found implementation to work some value conversions are needed
    CourseDeg = calc_course(MoveX, MoveY),
    CourseRad = deg_to_rad(CourseDeg), % todo: cleanup: in course the calculated values are rad anyway, converting to deg is just an extra calculation
    Distance = calc_distance(MoveX, MoveY),
    DistanceDeg = calc_degrees_per_meter_at_equator(NewCurrAlt) * Distance,
    DistanceRad = deg_to_rad(DistanceDeg),
    Lat1Rad = deg_to_rad(LastCurrLat),
    Lon1Rad = deg_to_rad(LastCurrLong),

    LatRad = math:asin(math:sin(Lat1Rad) * math:cos(DistanceRad) + math:cos(Lat1Rad) * math:sin(DistanceRad) * math:cos(CourseRad)),
    Dlon = math:atan2(math:sin(CourseRad) * math:sin(DistanceRad) * math:cos(Lat1Rad), math:cos(DistanceRad) - math:sin(Lat1Rad) * math:sin(LatRad)),
    LonRad = remainder((Lon1Rad + Dlon + math:pi()), (2 * math:pi())) - math:pi(),

    NewCurrLat = rad_to_deg(LatRad),
    NewCurrLong = rad_to_deg(LonRad),

    {NewCurrLat, NewCurrLong, NewCurrAlt}.

% some trigonometry
% returns angle between adjacent and hypotenuse, with MoveX as adjacent and MoveY as opposite
calc_course(MoveX, MoveY) ->
    case MoveX > 0 of
        true ->
            case MoveY > 0 of
                true ->
                    % tan(alpha) = opposite / adjacent
                    % arc tan to get the alpha
                    % erlang returns radians -> convert to degrees
                    Deg = rad_to_deg(math:atan(MoveY / MoveX));
                false ->
                    Temp = 360 - rad_to_deg(math:atan((MoveY * -1) / MoveX)),
                    case Temp == 360 of
                        true ->
                            Deg = 0.0;
                        false ->
                            Deg = Temp
                    end
            end;
        false ->
            % attention! MoveX not > 0 -> can be 0 -> div by zero
            case MoveX == 0 of
                true ->
                    case MoveY > 0 of
                        true ->
                            Deg = 90.0;
                        false ->
                            case MoveY == 0 of
                                true ->
                                    Deg = 0.0;
                                false ->
                                    Deg = 270.0
                            end
                    end;
                false -> % MoveX < 0
                    case MoveY > 0 of
                        true ->
                            Deg = 180 - rad_to_deg(math:atan(MoveY / (MoveX * -1)));
                        false ->
                            Deg = 180 + rad_to_deg(math:atan((MoveY * -1) / (MoveX * -1)))
                    end
            end
    end,
    Deg.

rad_to_deg(X) ->
    X * 180 / math:pi().
deg_to_rad(X) ->
    X * math:pi() / 180.

% distance = hypetenuse in Pythagorean theorem
calc_distance(MoveX, MoveY) ->
    math:sqrt(math:pow(MoveX,2) + math:pow(MoveY,2)).

calc_degrees_per_meter_at_equator(Alt) ->
    Circumference = 2 * math:pi() * Alt,
    360 / Circumference.

% erlang rem only operates with integers
% https://stackoverflow.com/questions/9297424/stdremainder-in-erlang
remainder(A, B) ->
    A_div_B = A / B,
    N = round(A_div_B),

    case (abs(N - A_div_B) == 0.5) of
    true ->
        A_div_B_Trunc = trunc(A_div_B),

        New_N = case ((abs(A_div_B_Trunc) rem 2) == 0) of
        true -> A_div_B_Trunc;
        false ->
            case (A_div_B >= 0) of
            true -> A_div_B_Trunc + 1;
            false -> A_div_B_Trunc - 1
            end
        end,
        A - New_N * B;
    false ->
        A - N * B
    end.

Tests: Object at lat/lon/alt (10°,10°,6371000m). No movement (0m,0m,0m). Expected:

{1.000000e+001,1.000000e+001,6.371000e+006}

But returns:

{1.000000e+001,-3.500000e+002,6.371000e+006}

360 degree less than expected...

Object at (0°,10°,6371000m), moving (10m,0m,0m). Expected:

{0.000000e+000,1.000009e+001,6.371000e+006}

Not sure if some digits are just not being displayed. Should be something like 10.000089932160591 as longitude. Anyway - returns:

{8.993216e-005,-3.500000e+002,6.371000e+006}

Same wrong longitude value although we're moving now? And a changed latitude value although we didn't move on Y-Axis?

What about same position, now moving 5,000,000m east?

{0.000000e+000,5.496608e+001,6.371000e+006} % expected
{4.496608e+001,-3.500000e+002,6.371000e+006} % returned

Implementation 3:

calc_position3(LastCurrLat, LastCurrLong, LastCurrAlt, MoveX, MoveY, MoveZ) ->
    {CurrX, CurrY, CurrZ} = spherical_to_cartesian(LastCurrLat, LastCurrLong, LastCurrAlt),
    NewX = CurrX + MoveX,
    NewY = CurrY + MoveY,
    NewZ = CurrZ + MoveZ,
    {NewCurrLat, NewCurrLong, NewCurrAlt} = cartesian_to_spherical(NewX, NewY, NewZ),
    {NewCurrLat, NewCurrLong, NewCurrAlt}.

spherical_to_cartesian(Lat, Lon, Alt) ->
    X = Alt * math:cos(Lat) * math:cos(Lon),
    Y = Alt * math:cos(Lat) * math:sin(Lon),
    Z = Alt * math:sin(Lat),
    {X, Y, Z}.

cartesian_to_spherical(X, Y, Z) ->
    R = math:sqrt(math:pow(X,2) + math:pow(Y,2)),
    Alt = math:sqrt(math:pow(X,2) + math:pow(Y,2) + math:pow(Z,2)),
    Lat = math:asin(Z / Alt),
    case R > 0 of
        true ->
            Lon = math:acos(X / R);
        false -> % actually: if R == 0, but it can never be negative (see above)
            Lon = 0
    end,
    {Lat, Lon, Alt}.

Tests like above:

Object at (10°,10°,6371000m), no movement

{1.000000e+001,1.000000e+001,6.371000e+006} % expected
{-5.752220e-001,5.752220e-001,6.371000e+006} % returned

At (0°,10°,6371000m), moving (10m,0m,0m)

{0.000000e+000,1.000009e+001,6.371000e+006} % expected
{0.000000e+000,2.566370e+000,6.370992e+006} % returned

At (0°,10°,6371000m), moving (5000000m,0m,0m)

{0.000000e+000,5.496608e+001,6.371000e+006} % expected
{0.000000e+000,1.670216e+000,3.483159e+006} % returned

So: Did I miss some rad to deg conversions or something similar?

P.S.: Sorry for the bad syntax highlighting, but Erlang doesn't seem to be available, so I took Shellscript. Makes it a bit more readable.

Answer 1

Let's assume:

• the equator is in the x-y plane
• latitude 0, longitude 0, altitude > 0 is on the +x axis
• the poles are at longitude 0

With these assumptions, to convert from spherical (lat, lng, alt) to Cartesian (x, y, z):

x = alt * cos(lat) * cos(lng)
y = alt * cos(lat) * sin(lng)
z = alt * sin(lat)

To convert from Cartesian (x, y, z) to spherical (lat, lng, alt):

r = sqrt(x*x + y*y)
alt = sqrt(x*x + y*y + z*z)
lat = arcsin(z / alt)
       / arccos(x / r) if r > 0
lng = -|
       \ 0 if r == 0

Then follow @Peter's suggestion:

(lat, lng, alt) = toSpherical(movement + toCartesian(lat, lng, alt))

Answer 2

From your answers to the comments I finally have enough information:

You have given a position of an object in geographical latitude, longitude (something equal or similar to WGS84 lat /long) and altitude value. Where altitude is measured from center of earth.

latitude : [-90 , 90] geographical latitude in decimal degrees, 90° is Northpole<br>
longitude: [-180, 180] longitude in decimal degrees, 0° is Greenwhich
altitude: [0 , INF ] in meters from center of earth.

So these coords are defined as Spherical coordinates, but beware of clockwise vs. counter clockwise (see later)

Further you have given a movement vector (dx,dy,dz) measured in meters which defines the relative movement from current position.
These vector is defined as cartesian vector.

Your application is neither for navigation, nor flight controll, its more in the area of games.

For the calculation you must know where the x,y,z Achsis are related. You should use the same Axis alignment that ECEF uses.

You need to perform these steps:

You first must know that geographical degrees are not mathematical degrees: Sometimes its is important to know that in mathematics degrees are counter clockwise, where in geography they are clockwise measured.
After that your coords must be converted to radians.

Convert the Spherical Coordinates to cartesian space. (Spherical to cartesian transformation: see http://en.wikipedia.org/wiki/Spherical_coordinate_system) But use Math.atan2() not atan().

(You might check your code with http://www.random-science-tools.com/maths/coordinate-converter.htm) Check if this conversion is used in the answer of Ted Hopp, too.

After the transformation you have a coordinate in an cartesian x,y,z space with unit = 1m.

Next is the cartesian movement vector: Be sure, or convert it, such that you have the correct axis alignment, that dx corresponds to your x, -dx to -x, and for y,dy, z,dz, too.

Then add the point to dx:

p2.x = p.x + dx;
p2.y = p.y + dy;
p2.z = p.z + dz;

Then reconvert back to spherical coordinates, as shown in the Wiki link or by Ted Hopp.

The other link (to stackoverflow) that you have given in your comments do other types of conversion, they will not work for your problem.

I recomend one special test case: set a spherical point to longitude = 179.9999999 then add 100m, the result should be -180.0000xx + something after the comma (Two jet fighters crashed for this problem when traversing the Datum limit)