计算具有当前纬度/经度/高度的对象的纬度/经度/高度值 + x、y、z 上的移动 [关闭]

Question

我原以为在 Google 上搜索时会有很多结果，但事实证明，在大多数情况下，就像在 stackoverflow 上一样，问题非常具体，包括 Google 地图、一些 GIS 或现实世界椭圆形。

就我而言，

• 我有一个具有 lat、long 和 alt 值的对象，其中 lat 和 long 是以度为单位的浮点数，alt 是作为十进制值的浮点数，表示距球体中心的距离。所以没有分钟（10°30' 是 10,5）。
• 对象在所有 3 个轴 x、y 和 z 上移动，其中移动相对于其当前位置
• 我需要计算新的 lat、long 和 alt 值

没有多想，我先是这样实现的：

只需将 z 运动添加到 alt，然后计算虚拟球体的周长以获得每单位（米）的度数。有了这个，我首先计算了新的纬度，然后再计算了新的。我知道一个接一个地计算一个值是错误的，但只要我的“对象世界”中的所有计算都以相同的方式完成，整体行为就可以了。我考虑了一些事情，比如当对象绕整个球体旋转时，long 值不会改变，并且绕球体绕半圈在 x 轴上与在 y 轴上不同（x 轴：-180 到 180，y-轴 -90 到 90）和类似的东西，它工作。

但后来我意识到我没有考虑我计算赤道上每米的度数，也没有考虑其他纬度值。那时我知道事情更复杂，我开始在网上搜索。但是我没有找到适合我需要的算法。现在我确信这已经做过很多次了，所以我在这里询问是否有人以前处理过这个话题并且可以指出一个很好的实现:)。

我发现是这样的：

• 将纬度/经度转换为墨卡托投影的算法
• Haversine 公式，用于计算两个纬度/经度值对的距离
• 其他用途的其他配方

这对我没有帮助。

对我帮助最大的是：Calculate latitude and longitude having meters distance from another latitude/longitude point

但我想我还没有完全理解它。

• 什么是弧度的课程？ （关于我有 x 和 y 动作）
• 外推法不考虑海拔/地球半径？

如果有人能帮我解决这个问题，那就太棒了！

（旁注：我在 Erlang 中实现了这个，但这没关系，任何类型的算法都会有所帮助）

---

更新

我实现了上面提到的功能和下面答案中的功能。测试时我得到了错误的值，可能是因为我在实现中犯了错误或计算了错误的测试数据。让我们看看：

实施 1：

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

测试（赤道）正常。

实施 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.

测试： 对象位于 lat/lon/alt (10°,10°,6371000m)。无运动（0m,0m,0m）。预期：

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

但返回：

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

比预期少 360 度...

对象在 (0°,10°,6371000m)，移动 (10m,0m,0m)。预期：

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

不确定是否只是没有显示某些数字。应该是像 10.000089932160591 这样的经度。无论如何 - 返回：

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

虽然我们现在要移动，但经度值仍然错误？虽然我们没有在 Y 轴上移动，但纬度值发生了变化？

同样的位置，现在向东移动 5,000,000 米呢？

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

实施 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}.

类似上面的测试：

物体在 (10°,10°,6371000m)，没有移动

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

在 (0°,10°,6371000m)，移动 (10m,0m,0m)

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

在 (0°,10°,6371000m)，移动 (5000000m,0m,0m)

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

那么：我错过了一些 rad 到 deg 的转换或类似的东西吗？

P.S.：抱歉语法高亮不好，但是 Erlang 似乎不可用，所以我使用了 Shellscript。使其更具可读性。

Answer 1

从你对 cme​​ts 的回答中，我终于有了足够的信息：

您已在地理纬度、经度（等于或类似于 WGS84 纬度/经度）和高度值中给出了对象的位置。高度是从地心开始测量的。

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.

所以这些坐标被定义为球坐标，但要注意顺时针与逆时针（见下文）

您还给出了一个以米为单位的运动矢量 (dx,dy,dz)，它定义了与当前位置的相对运动。
这些矢量被定义为笛卡尔矢量。

您的应用程序既不是导航，也不是飞行控制，更多的是游戏领域。

对于计算，您必须知道 x,y,z Achsis 的相关位置。您应该使用与ECEF 相同的轴对齐方式。

您需要执行以下步骤：

您首先必须知道地理学位不是数学学位：有时重要的是要知道数学学位是逆时针方向的，而地理学位是顺时针测量的。
之后，您的坐标必须转换为弧度。

将球坐标转换为笛卡尔空间。 （球面到笛卡尔变换：见http://en.wikipedia.org/wiki/Spherical_coordinate_system） 但是使用 Math.atan2() 而不是 atan()。

（您可以使用http://www.random-science-tools.com/maths/coordinate-converter.htm 检查您的代码） 检查 Ted Hopp 的回答中是否也使用了这种转换。

转换后，您在笛卡尔 x,y,z 空间中有一个坐标，单位 = 1m。

接下来是笛卡尔运动矢量：确保或转换它，以便您有正确的轴对齐方式，dx 对应于您的 x，-dx 对应于 -x，并且也对应于 y,dy,z,dz .

然后将点添加到dx：

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

然后重新转换回球坐标，如 Wiki 链接或 Ted Hopp 所示。

您在 cmets 中提供的另一个链接（到 stackoverflow）进行其他类型的转换，它们不适用于您的问题。

我推荐一个特殊的测试用例：将球点设置为经度 = 179.9999999 然后加上 100m，结果应该是 -180.0000xx + 逗号后面的东西（两架喷气式战斗机在穿越基准限制时因这个问题坠毁）

Answer 2

假设：

• 赤道在 x-y 平面内
• 纬度 0，经度 0，高度 > 0 在 +x 轴上
• 两极位于经度 0

在这些假设下，从球形（lat、lng、alt）转换为笛卡尔（x、y、z）：

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

从笛卡尔 (x, y, z) 转换为球面 (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

那就听从@Peter的建议吧：

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