0

我以为在 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)和类似的东西,它工作。

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

我发现的是这样的:

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

这对我没有帮助。

对我帮助最大的是:计算距离另一个纬度/经度点有米远的纬度和经度

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

  • 弧度的课程是什么?(关于我有 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.

测试:对象在纬度/经度/高度(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,000m 呢?

{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 的转换或类似的东西吗?

PS:对不起语法高亮不好,但是Erlang似乎不可用,所以我拿了Shellscript。使它更具可读性。

4

2 回答 2

1

让我们假设:

  • 赤道在xy平面
  • 纬度 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))
于 2012-12-04T21:42:08.920 回答
1

从您对评论的回答中,我终于有了足够的信息:

您已在地理纬度、经度(等于或类似于 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.

所以这些坐标被定义为Spherical coordinate,但要注意顺时针与逆时针(见下文)

此外,您还给出了一个以米为单位的运动矢量 (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 所示。

您在评论中提供的另一个链接(到 stackoverflow)进行其他类型的转换,它们不适用于您的问题。

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

于 2012-12-09T14:02:39.720 回答