我以为在 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。使它更具可读性。