如何找到给定纬度/经度以北 x 公里的纬度/经度？
Posted
技术标签:
【中文标题】如何找到给定纬度/经度以北 x 公里的纬度/经度？【英文标题】：How do I find the lat/long that is x km north of a given lat/long? 【发布时间】：2010-11-10 15:24:02 【问题描述】：我有一些生成谷歌地图的 C# 代码。此代码查看我需要在地图上绘制的所有点，然后计算出矩形的边界以包含这些点。然后它将此边界传递给 Google Maps API 以适当地设置缩放级别以显示地图上的所有点。
此代码运行良好，但我有一个新要求。
其中一个点可能具有与之相关的精度。如果是这种情况，那么我在该点周围画一个圆，半径设置为精度值。这再次工作正常，但是我的边界检查现在没有做我想要做的事情。我想让边界框包含完整的圆圈。
这需要一种算法来获取点 x 并计算点 y，该点位于 x 以北 z 米处和 x 以南 z 米处。
有没有人有这个算法，最好是在 C# 中。我确实找到了一个通用算法here，但我似乎没有正确实现这一点，因为我得到的答案是漂移了 1000 公里。
这是通用示例
Lat/lon given radial and distance
A point lat,lon is a distance d out on the tc radial from point 1 if:
lat=asin(sin(lat1)*cos(d)+cos(lat1)*sin(d)*cos(tc))
IF (cos(lat)=0)
lon=lon1 // endpoint a pole
ELSE
lon=mod(lon1-asin(sin(tc)*sin(d)/cos(lat))+pi,2*pi)-pi
ENDIF
这是我的 C# 翻译。
// Extend a Point North/South by the specified distance
public static Point ExtendPoint(Point _pt, int _distance, int _bearing )
Decimal lat = 0.0;
Decimal lng = 0.0;
lat = Math.Asin(Math.Sin(_pt.Lat) * Math.Cos(_distance) + Math.Cos(_pt.Lat) *
Math.Sin(_distance) * Math.Cos(_bearing));
if (Math.Cos(lat) == 0)
lng = _pt.Lng; // endpoint a pole
else
lng = (
(_pt.Lng - Math.Asin(Math.Sin(_bearing) * Math.Sin(_distance) / Math.Cos(lat))
+ Math.PI) % (2 * Math.PI)) - Math.PI;
ret = new Point(lat,lng);
return ret;
我调用这个函数时，方位角为 0 来计算新的北纬位置，值 180 来计算新的南纬位置。
谁能看到我做错了什么或者提供一个已知的工作算法？
【问题讨论】：
对于它的价值，我有一个 php 示例，它可以执行 OP 请求的操作。在我的示例中，它在起始纬度/经度坐标周围绘制了一个框，但代码可以很容易地用于获取单个点、X 公里或英里以外的点。 richardpeacock.com/blog/2011/11/… 【参考方案1】：对于想要java版Eirch的代码的人
/**
* move latlng point by rang and bearing
*
* @param latLng point
* @param range range in meters
* @param bearing bearing in degrees
* @return new LatLng
*/
public static LatLng moveLatLng(LatLng latLng, double range, double bearing)
double EarthRadius = 6378137.0;
double DegreesToRadians = Math.PI / 180.0;
double RadiansToDegrees = 180.0 / Math.PI;
final double latA = latLng.latitude * DegreesToRadians;
final double lonA = latLng.longitude * DegreesToRadians;
final double angularDistance = range / EarthRadius;
final double trueCourse = bearing * DegreesToRadians;
final double lat = Math.asin(
Math.sin(latA) * Math.cos(angularDistance) +
Math.cos(latA) * Math.sin(angularDistance) * Math.cos(trueCourse));
final double dlon = Math.atan2(
Math.sin(trueCourse) * Math.sin(angularDistance) * Math.cos(latA),
Math.cos(angularDistance) - Math.sin(latA) * Math.sin(lat));
final double lon = ((lonA + dlon + Math.PI) % (Math.PI * 2)) - Math.PI;
return new LatLng(lat * RadiansToDegrees, lon * RadiansToDegrees);
【讨论】：
【参考方案2】：如果您有给定的纬度和经度，您可以像这样计算 x 公里纬度变化的正确纬度和经度：
new-lat = ((old-km-north + x-km-change)/40,075) * 360)
^ is the ratio of the ^ times the ratio of the circle
earth the change by 360 to get the total ratio
covers. covered in degrees.
这同样适用于经度。如果你有总距离加上变化，你可以用类似的方式计算总度数。
new-long = ((old-km-east + x-km-change)/40,075) * 360)
^ is the ratio of the ^ times the ratio of the circle
earth the change by 360 to get the total ratio
covers. covered in degrees.
同样，这些计算应该有效，但我在这里失去了纯粹的直觉，但逻辑似乎确实成立。
编辑：正如 Skizz 所指出的，需要使用 2.pi.r.cos(lat) 或 40074.cos(lat) 将 40,075 调整到任何给定纬度的地球圆周
【讨论】：
嗯，这看起来容易多了 :-) 谢谢。我现在试试。什么是 40.075 幻数。是 KM 覆盖的秒数（度）吗？ 那是地球的圆周，尽管试一试。这不符合直觉，因此可能也需要一些调整。 KevDog 的意思是 40,075 随着纬度的增加而减少 - 极点以东 1 公里将增加赤道以东 1 公里以上的经度。因此，您需要将 40,075 替换为给定纬度的周长，即：2.pi.rs.cos(lat)，其中 rs 是地球的半径，lat 是纬度。它确实假设地球是球形的。但是，如果您向东 1 公里然后向北 1 公里而不是向北 1 公里然后向东 1 公里，这种方法会产生不同的结果。 如果您只是在原始位置添加度数，为什么不保持经度不变（因为他只是向北，经度不应该改变）并且只需将 km-change/40,075 添加到 lat ?这个方程对我来说似乎很脆弱并且缺乏大量数据。当有非常详尽的方程式描述如何计算时，为什么还要凭直觉？ 谁能告诉我old-km-north/ old-km-east指的是什么？【参考方案3】：对于懒惰的人，（像我一样；））复制粘贴解决方案，Erich Mirabal 的版本，改动很小：
using System.Device.Location; // add reference to System.Device.dll
public static class GeoUtils
/// <summary>
/// Calculates the end-point from a given source at a given range (meters) and bearing (degrees).
/// This methods uses simple geometry equations to calculate the end-point.
/// </summary>
/// <param name="source">Point of origin</param>
/// <param name="range">Range in meters</param>
/// <param name="bearing">Bearing in degrees</param>
/// <returns>End-point from the source given the desired range and bearing.</returns>
public static GeoCoordinate CalculateDerivedPosition(this GeoCoordinate source, double range, double bearing)
var latA = source.Latitude * DegreesToRadians;
var lonA = source.Longitude * DegreesToRadians;
var angularDistance = range / EarthRadius;
var trueCourse = bearing * DegreesToRadians;
var lat = Math.Asin(
Math.Sin(latA) * Math.Cos(angularDistance) +
Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse));
var dlon = Math.Atan2(
Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA),
Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat));
var lon = ((lonA + dlon + Math.PI) % (Math.PI*2)) - Math.PI;
return new GeoCoordinate(
lat * RadiansToDegrees,
lon * RadiansToDegrees,
source.Altitude);
private const double DegreesToRadians = Math.PI/180.0;
private const double RadiansToDegrees = 180.0/ Math.PI;
private const double EarthRadius = 6378137.0;
用法：
[TestClass]
public class CalculateDerivedPositionUnitTest
[TestMethod]
public void OneDegreeSquareAtEquator()
var center = new GeoCoordinate(0, 0);
var radius = 111320;
var southBound = center.CalculateDerivedPosition(radius, -180);
var westBound = center.CalculateDerivedPosition(radius, -90);
var eastBound = center.CalculateDerivedPosition(radius, 90);
var northBound = center.CalculateDerivedPosition(radius, 0);
Console.Write($"leftBottom: southBound.Latitude , westBound.Longitude rightTop: northBound.Latitude , eastBound.Longitude");
【讨论】：
【参考方案4】：Ed William's rather awesome site 上的两个方程有问题...但我没有分析它们以了解原因。
I found here 的第三个等式似乎给出了正确的结果。
这是 php 中的测试用例...第三个等式是正确的，前两个给出的经度值非常不正确。
<?php
$lon1 = -108.553412; $lat1 = 35.467155; $linDistance = .5; $bearing = 170;
$lon1 = deg2rad($lon1); $lat1 = deg2rad($lat1);
$distance = $linDistance/6371; // convert dist to angular distance in radians
$bearing = deg2rad($bearing);
echo "lon1: " . rad2deg($lon1) . " lat1: " . rad2deg($lat1) . "<br>\n";
// doesn't work
$lat2 = asin(sin($lat1) * cos($distance) + cos($lat1) * sin($distance) * cos($bearing) );
$dlon = atan2(sin($bearing) * sin($distance) * cos($lat1), cos($distance) - sin($lat1) * sin($lat2));
$lon2 = (($lon1 - $dlon + M_PI) % (2 * M_PI)) - M_PI; // normalise to -180...+180
echo "lon2: " . rad2deg($lon2) . " lat2: " . rad2deg($lat2) . "<br>\n";
// same results as above
$lat3 = asin( (sin($lat1) * cos($distance)) + (cos($lat1) * sin($distance) * cos($bearing)));
$lon3 = (($lon1 - (asin(sin($bearing) * sin($distance) / cos($lat3))) + M_PI) % (2 * M_PI)) - M_PI;
echo "lon3: " . rad2deg($lon3) . " lat3: " . rad2deg($lat3) . "<br>\n";
// gives correct answer... go figure
$lat4 = asin(sin($lat1) * cos($linDistance/6371) + cos($lat1) * sin($linDistance/6371) * cos($bearing) );
$lon4 = $lon1 + atan2( (sin($bearing) * sin($linDistance/6371) * cos($lat1) ), (cos($linDistance/6371) - sin($lat1) * sin($lat2)));
echo "lon4: " . rad2deg($lon4) . " lat4: " . rad2deg($lat4) . "<br>\n";
?>
注意，我收到了前两个方程的作者 (Ed Williams) 的电子邮件：
来自我的“实施笔记”：
注意 mod 功能。这似乎在 不同的语言，对符号是否有不同的约定 结果遵循除数或被除数的符号。 （我们想要标志 跟随除数或成为欧几里得。 C 的 fmod 和 Java 的 % 不起作用。） 在本文档中，Mod(y,x) 是 y 除以 x 的余数，并且始终 位于 0
如果您有一个 floor 函数（Excel 中的 int），则返回 floor(x)= “小于或等于 x 的最大整数”，例如地板（-2.3）=-3 和 地板(2.3) =2
mod(y,x) = y - x*floor(y/x)
以下应该在没有地板功能的情况下工作 - 不管 "int" 是否向下截断或舍入：
mod=y - x * int(y/x)
if ( mod < 0) mod = mod + x
php 就像 C 中的 fmod，但出于我的目的，它是“错误的”。
【讨论】：
【参考方案5】：我不确定我是否在这里遗漏了什么，但我认为这个问题可以改写为，“我有一个纬度/经度点，我想找到北 x 米和南 x 米的点那个点。”
如果这是问题，那么您不需要找到新的经度（这使事情变得更简单），您只需要一个新的纬度。一个纬度大约是地球上任何地方的 60 海里长，而一海里是 1,852 米。所以，对于新的纬度，南北 x 米：
north_lat = lat + x / (1852 * 60)
north_lat = min(north_lat, 90)
south_lat = lat - x / (1852 * 60)
south_lat = max(south_lat, -90)
这并不完全准确，因为地球不是一个完美的球体，每个纬度之间正好有 60 海里。但是，其他答案假设纬度线是等距的，所以我假设您不关心这一点。如果您对可能引入多少错误感兴趣，***上有一个很好的表格，在此链接上显示了不同纬度的“纬度每 1° 变化的表面距离”：
http://en.wikipedia.org/wiki/Latitude#Degree_length
【讨论】：
【参考方案6】：我有一段非常相似的代码。与另一个实现相比，它让我得到了非常接近的结果。
我认为你的问题是你使用“距离”作为米的线性距离，而不是弧度的角距离。
/// <summary>
/// Calculates the end-point from a given source at a given range (meters) and bearing (degrees).
/// This methods uses simple geometry equations to calculate the end-point.
/// </summary>
/// <param name="source">Point of origin</param>
/// <param name="range">Range in meters</param>
/// <param name="bearing">Bearing in degrees</param>
/// <returns>End-point from the source given the desired range and bearing.</returns>
public static LatLonAlt CalculateDerivedPosition(LatLonAlt source, double range, double bearing)
double latA = source.Latitude * UnitConstants.DegreesToRadians;
double lonA = source.Longitude * UnitConstants.DegreesToRadians;
double angularDistance = range / GeospatialConstants.EarthRadius;
double trueCourse = bearing * UnitConstants.DegreesToRadians;
double lat = Math.Asin(
Math.Sin(latA) * Math.Cos(angularDistance) +
Math.Cos(latA) * Math.Sin(angularDistance) * Math.Cos(trueCourse));
double dlon = Math.Atan2(
Math.Sin(trueCourse) * Math.Sin(angularDistance) * Math.Cos(latA),
Math.Cos(angularDistance) - Math.Sin(latA) * Math.Sin(lat));
double lon = ((lonA + dlon + Math.PI) % UnitConstants.TwoPi) - Math.PI;
return new LatLonAlt(
lat * UnitConstants.RadiansToDegrees,
lon * UnitConstants.RadiansToDegrees,
source.Altitude);
在哪里
public const double EarthRadius = 6378137.0; // WGS-84 ellipsoid parameters
并且 LatLonAlt 以度/米为单位（转换发生在内部）。 根据需要进行调整。
我假设您可以弄清楚UnitConstants.DegreesToRadians 的值是什么:)
【讨论】：
我在 java 中尝试这种方法，但它不起作用。我得到 dlon = 0 的值。这是我的 dlon 代码：double dlon = Math.atan2(Math.sin(brng)*Math.sin(dist)*Math.cos(lat1), Math.cos(dist) -Math.sin(lat1)*Math.sin(lat2));我已经将纬度和经度的值转换为弧度。每个代码都按照您说的编写。 我无法计算出 UnitConstants.DegreesToRadians、GeospatialConstants.EarthRadius、UnitConstants.TwoPi 和 UnitConstants.RadiansToDegrees 是什么？ :( 有人可以帮忙吗？ EarthRadius = 6378137.0, TwoPi = Math.PI * 2, DegreesToRadians = 0.0174532925, RadiansToDegrees = 57.2957795。【参考方案7】：如果先重投影到UTM再看距离会更准确。
希望对你有帮助
【讨论】：
以上是关于如何找到给定纬度/经度以北 x 公里的纬度/经度？的主要内容，如果未能解决你的问题，请参考以下文章
给定 GPS 的地理位置，如何在给定的 x 米半径内找到纬度和经度