问题描述
我认为在Google上搜索时会得到很多结果，但事实证明，在大多数情况下，就像此处的stackoverflow一样，问题非常具体，包括Google Maps，一些GIS或带有椭圆形。
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.
在我的情况下，
- 我有一个具有lat，long和alt值的对象，其中lat和long是以度为单位的浮点数，而alt是以十进制值表示的浮点数，代表与球体中心的距离。因此没有分钟（10°30'为10,5）。
- 对象在x，y和z的所有3个轴上移动相对于其当前位置
- 我需要计算新的纬度，经度和alt值
- 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:
只需将z移动添加到alt，然后计算虚拟球体的圆周得到单位（米）的度数。这样，我首先计算了新纬度，然后又计算了很长时间。我知道一个接一个地计算值是错误的，但是只要我对象世界中的所有计算都以相同的方式进行，则总体行为就可以了。我考虑了一些问题，例如当对象绕整个球体移动时，长值不变，并且绕球体的一半在x轴上的方向与在y轴上的方向不同（x轴：-180至180，y-
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 :) .
我发现的是：
- 将经纬度转换为墨卡托投影的算法
- 用于计算与两个纬度/经度值对之间的距离的Haversine公式
- 其他用于其他目的的公式
没有帮助我。
对我最有用的是：
但是我想我还不太了解。
But I think I didn't fully understand it yet.
- 弧度的方向是什么？ （关于我有x和y运动）
- 外推方法是否不考虑海拔高度/地球半径？
如果有人可以帮助我，那就太好了！
（旁注：我正在实现此功能在Erlang中，但这没关系，任何一种算法都可以帮助）
(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:
实现1：
% My own implementation that only works on the equator so far. Simple calculations as mentioned above.
（在赤道上）测试正常。
Tests (on equator) ok.
实现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）。预期：
Tests:Object at lat/lon/alt (10°,10°,6371000m). No movement (0m,0m,0m). Expected:
{1.000000e+001,1.000000e+001,6.371000e+006}
但返回：
{1.000000e+001,-3.500000e+002,6.371000e+006}
比预期的少360度。 ..
360 degree less than expected...
对象位于（0°，10°，6371000m），正在移动（10m，0m，0m）。预期：
Object at (0°,10°,6371000m), moving (10m,0m,0m). Expected:
{0.000000e+000,1.000009e+001,6.371000e+006}
不确定是否仅显示了一些数字。经度应类似于10.000089932160591。无论如何-返回值：
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}
相同的经度值，尽管我们现在要移动？尽管我们没有沿Y轴移动，但纬度值却发生了变化？
Same wrong longitude value although we're moving now? And a changed latitude value although we didn't move on Y-Axis?
相同位置怎么办，现在向东移动5,000,000m？
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
实现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），没有运动
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
在（0°，10°，6371000m）处移动（10m，0m，0m）
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
在（0°，10°，6371000m）处移动（5000000m，0m，0m）
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
所以：我是否错过了rad到deg转换或类似的东西？
So: Did I miss some rad to deg conversions or something similar?
PS：抱歉，语法突出显示错误，但是Erlang似乎不可用，所以我选择了Shellscript。
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.
推荐答案
从您对评论的回答中，我终于有了足够的信息：
From your answers to the comments I finally have enough information:
已根据地理纬度，经度（等于或类似于WGS84纬度/经度）和海拔高度给出了对象的位置。
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)
您还给出了以米为单位的运动矢量（dx，dy，dz），该运动矢量定义了从当前位置开始的相对运动。 >这些矢量被定义为笛卡尔矢量。
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.
要进行计算，您必须知道x，y，z Achsis与哪里相关。您应该使用使用的Axis对齐方式。
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 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.
转换笛卡尔坐标系的球坐标。 （从球形到笛卡尔的转换：请参见）
，但请使用Math.atan2（）而不是atan（）。
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().
（您可以使用）
检查Ted Hopp的答案中是否也使用了这种转换。
(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.
转换后，您在x，y，z直角坐标系中的坐标为unit = 1m。
After the transformation you have a coordinate in an cartesian x,y,z space with unit = 1m.
下一个是笛卡尔运动矢量：请确保或对其进行转换，以使您具有正确的轴对齐，并且dx对应于x，-dx对应于-x，对于y，dy，z，dz也是如此。
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.
然后将点添加到dx：
p2.x = p.x + dx;
p2.y = p.y + dy;
p2.z = p.z + dz;
然后重新转换回球形坐标，如Wiki链接或Ted Hopp所示。
Then reconvert back to spherical coordinates, as shown in the Wiki link or by Ted Hopp.
您在注释中给出的另一个链接（指向stackoverflow）进行其他类型的转换，它们将无法解决您的问题。
The other link (to stackoverflow) that you have given in your comments do other types of conversion, they will not work for your problem.
我推荐一个特殊的测试用例：将球形点设置为经度= 179.9999999，然后添加100m，结果应为-180.0000xx +逗号后的值（两个喷气式战斗机在越过基准极限时因该问题而坠毁）
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)
这篇关于使用当前lat / long / alt +在x，y，z上移动来计算对象的lat / long / alt值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！