Lat/Lng 点到 Minor Arc 段的距离
Distance from Lat/Lng point to Minor Arc segment
我需要计算从 lat/lng GPS 点 P 到另外 2 个 lat/lng GPS 点 A 和 B 描述的线段的最短距离。
'Cross-track distance'帮我算出P到A、B所描述的大圆的最短距离
然而，这不是我想要的。我需要 P 和 A-B 的 线段 之间的距离，而不是整个大圆。
我使用了 http://www.movable-type.co.uk/scripts/latlong.html
中的以下实现
Formula: dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R
where:
δ13 is (angular) distance from start point to third point
θ13 is (initial) bearing from start point to third point
θ12 is (initial) bearing from start point to end point
R is the earth’s radius
以下图片有望展示我正在尝试解决的问题：
在第一张图片中，绿色线表示的交叉轨道距离是正确的，确实是到线段 AB 的最短距离。
在第二张图中显示了交叉轨道距离的问题，在这种情况下，我希望最短距离是简单距离 AP，但是交叉轨道距离给出了 [=34 指示的距离=]红色行。
如何更改我的算法以考虑到这一点，或检查点 X 是否在 AB 内。有可能通过计算来做到这一点吗？还是迭代是唯一可能的（昂贵的）解决方案？ （沿AB取N个点，计算P点到所有这些点的最小距离）
为了简单起见，图像中的所有线条都是直的。实际上，这些是大圆上的小圆弧
对于 100 - 1000m 的球形问题，很容易只转换为
笛卡尔坐标 space，使用等角投影。
然后继续学校数学：
使用函数 "distance from line segment" 很容易找到现成的实现。
此函数使用（有时 returns）直线 A、B 上的投影点 X 的相对 forward/backward 位置。值为
- 在区间 [0,1] 中，如果投影点在线段内。
- 如果X在A之前在外面，则为负数，
- 如果在B之后在外面则>1。
如果相对位置在 0,1 之间，则采用正常距离，如果在起点和线终点 A、B 的较短距离之外。
这种/或非常相似的笛卡尔实现的一个例子是Shortest distance between a point and a line segment
首先，一些术语：
我们的弧线是从 p1 到 p2 绘制的。
我们的第三点是p3。
与大圆相交的虚点是 p4。
p1 由 lat1,lon1 定义； p2 通过 lat2,lon2；等等
dis12 是 p1 到 p2 的距离；等等
bear12 是从 p1 到 p2 的方位角；等等
dxt 是交叉轨道距离。
dxa是跨弧距离，我们的目标！
注意交叉轨道公式依赖于相对方位，bear13-bear12
我们有3个案例要处理。
情况一：相对方位为钝角。所以，dxa=dis13.
案例2.1：相对方位是锐角，p4落在我们的弧上。
所以，dxa=dxt.
案例2.2：相对方位是锐角，p4超出了我们的范围。
所以，dxa=dis23
算法：
第一步： 如果相对方位角是钝角，dxa=dis13
完毕！
第 2 步： 如果相对方位角是锐角：
2.1: 查找 dxt。
2.3:寻找dis12。
2.4:寻找dis14。
2.4: 如果 dis14>dis12，dxa=dis23。
完毕！
2.5: 如果我们到达这里，dxa=abs(dxt)
MATLAB代码：
function [ dxa ] = crossarc( lat1,lon1,lat2,lon2,lat3,lon3 )
%// CROSSARC Calculates the shortest distance in meters
%// between an arc (defined by p1 and p2) and a third point, p3.
%// Input lat1,lon1,lat2,lon2,lat3,lon3 in degrees.
lat1=deg2rad(lat1); lat2=deg2rad(lat2); lat3=deg2rad(lat3);
lon1=deg2rad(lon1); lon2=deg2rad(lon2); lon3=deg2rad(lon3);
R=6371000; %// Earth's radius in meters
%// Prerequisites for the formulas
bear12 = bear(lat1,lon1,lat2,lon2);
bear13 = bear(lat1,lon1,lat3,lon3);
dis13 = dis(lat1,lon1,lat3,lon3);
diff = abs(bear13-bear12);
if diff > pi
diff = 2 * pi - diff;
end
%// Is relative bearing obtuse?
if diff>(pi/2)
dxa=dis13;
else
%// Find the cross-track distance.
dxt = asin( sin(dis13/R)* sin(bear13 - bear12) ) * R;
%// Is p4 beyond the arc?
dis12 = dis(lat1,lon1,lat2,lon2);
dis14 = acos( cos(dis13/R) / cos(dxt/R) ) * R;
if dis14>dis12
dxa=dis(lat2,lon2,lat3,lon3);
else
dxa=abs(dxt);
end
end
end
function [ d ] = dis( latA, lonA, latB, lonB )
%DIS Finds the distance between two lat/lon points.
R=6371000;
d = acos( sin(latA)*sin(latB) + cos(latA)*cos(latB)*cos(lonB-lonA) ) * R;
end
function [ b ] = bear( latA,lonA,latB,lonB )
%BEAR Finds the bearing from one lat/lon point to another.
b=atan2( sin(lonB-lonA)*cos(latB) , ...
cos(latA)*sin(latB) - sin(latA)*cos(latB)*cos(lonB-lonA) );
end
示例输出： 演示所有案例。请参阅下面的地图。
>> crossarc(-10.1,-55.5,-15.2,-45.1,-10.5,-62.5)
ans =
7.6709e+05
>> crossarc(40.5,60.5,50.5,80.5,51,69)
ans =
4.7961e+05
>> crossarc(21.72,35.61,23.65,40.7,25,42)
ans =
1.9971e+05
地图上那些相同的输出！：
演示案例 1：
演示案例 2.1：
演示案例 2.2：
归功于：http://www.movable-type.co.uk/scripts/latlong.html
对于公式
和：http://www.darrinward.com/lat-long/?id=1788764
用于生成地图图像。
/**
* Calculates the euclidean distance from a point to a line segment.
*
* @param v the point
* @param a start of line segment
* @param b end of line segment
* @return an array of 2 doubles:
* [0] distance from v to the closest point of line segment [a,b],
* [1] segment coeficient of the closest point of the segment.
* Coeficient values < 0 mean the closest point is a.
* Coeficient values > 1 mean the closest point is b.
* Coeficient values between 0 and 1 mean how far along the segment the closest point is.
*
* @author Afonso Santos
*/
public static
double[]
distanceToSegment( final R3 v, final R3 a, final R3 b )
{
double[] results = new double[2] ;
final R3 ab_ = b.sub( a ) ;
final double ab = ab_.modulus( ) ;
final R3 av_ = v.sub( a ) ;
final double av = av_.modulus( ) ;
if (ab == 0.0) // a and b coincide
{
results[0] = av ; // Distance
results[1] = 0.0 ; // Segment coeficient.
}
else
{
final double avScaProjAb = av_.dot(ab_) / ab ;
final double abCoeficient = results[1] = avScaProjAb / ab ;
if (abCoeficient <= 0.0) // Point is before start of the segment ?
results[0] = av ; // Use distance to start of segment.
else if (abCoeficient >= 1.0) // Point is past the end of the segment ?
results[0] = v.sub( b ).modulus() ; // Use distance to end of segment.
else // Point is within the segment's start/end perpendicular boundaries.
{
if (avScaProjAb >= av) // Test to avoid machine float representation epsilon rounding errors that would result in expection on sqrt.
results[0] = 0.0 ; // a, b and v are colinear.
else
results[0] = Math.sqrt( av * av - avScaProjAb * avScaProjAb ) ; // Perpendicular distance from point to segment.
}
}
return results ;
}
上述方法需要笛卡尔 3D space 个参数，而您要求使用 lat/lon 个参数。要进行转换，请使用
/**
* Calculate 3D vector (from center of earth).
*
* @param latDeg latitude (degrees)
* @param lonDeg longitude (degrees)
* @param eleMtr elevation (meters)
* @return 3D cartesian vector (from center of earth).
*
* @author Afonso Santos
*/
public static
R3
cartesian( final double latDeg, final double lonDeg, final double eleMtr )
{
return versor( latDeg, lonDeg ).scalar( EARTHMEANRADIUS_MTR + eleMtr ) ;
}
对于 3D/R3 代码的其余部分或如何计算到 path/route/track 检查的距离
https://sourceforge.net/projects/geokarambola/
将 Java 版本添加到 wdickerson 答案：
public static double pointToLineDistance(double lon1, double lat1, double lon2, double lat2, double lon3, double lat3) {
lat1 = Math.toRadians(lat1);
lat2 = Math.toRadians(lat2);
lat3 = Math.toRadians(lat3);
lon1 = Math.toRadians(lon1);
lon2 = Math.toRadians(lon2);
lon3 = Math.toRadians(lon3);
// Earth's radius in meters
double R = 6371000;
// Prerequisites for the formulas
double bear12 = bear(lat1, lon1, lat2, lon2);
double bear13 = bear(lat1, lon1, lat3, lon3);
double dis13 = dis(lat1, lon1, lat3, lon3);
// Is relative bearing obtuse?
if (Math.abs(bear13 - bear12) > (Math.PI / 2))
return dis13;
// Find the cross-track distance.
double dxt = Math.asin(Math.sin(dis13 / R) * Math.sin(bear13 - bear12)) * R;
// Is p4 beyond the arc?
double dis12 = dis(lat1, lon1, lat2, lon2);
double dis14 = Math.acos(Math.cos(dis13 / R) / Math.cos(dxt / R)) * R;
if (dis14 > dis12)
return dis(lat2, lon2, lat3, lon3);
return Math.abs(dxt);
}
private static double dis(double latA, double lonA, double latB, double lonB) {
double R = 6371000;
return Math.acos(Math.sin(latA) * Math.sin(latB) + Math.cos(latA) * Math.cos(latB) * Math.cos(lonB - lonA)) * R;
}
private static double bear(double latA, double lonA, double latB, double lonB) {
// BEAR Finds the bearing from one lat / lon point to another.
return Math.atan2(Math.sin(lonB - lonA) * Math.cos(latB), Math.cos(latA) * Math.sin(latB) - Math.sin(latA) * Math.cos(latB) * Math.cos(lonB - lonA));
}
并添加 python Sga 实现的翻译：
def bear(latA, lonA, latB, lonB):
# BEAR Finds the bearing from one lat / lon point to another.
return math.atan2(
math.sin(lonB - lonA) * math.cos(latB),
math.cos(latA) * math.sin(latB) - math.sin(latA) * math.cos(latB) * math.cos(lonB - lonA)
)
def pointToLineDistance(lon1, lat1, lon2, lat2, lon3, lat3):
lat1 = math.radians(lat1)
lat2 = math.radians(lat2)
lat3 = math.radians(lat3)
lon1 = math.radians(lon1)
lon2 = math.radians(lon2)
lon3 = math.radians(lon3)
R = 6378137
bear12 = bear(lat1, lon1, lat2, lon2)
bear13 = bear(lat1, lon1, lat3, lon3)
dis13 = distance( (lat1, lon1), (lat3, lon3)).meters
# Is relative bearing obtuse?
if math.fabs(bear13 - bear12) > (math.pi / 2):
return dis13
# Find the cross-track distance.
dxt = math.asin(math.sin(dis13 / R) * math.sin(bear13 - bear12)) * R
# Is p4 beyond the arc?
dis12 = distance((lat1, lon1), (lat2, lon2)).meters
dis14 = math.acos(math.cos(dis13 / R) / math.cos(dxt / R)) * R
if dis14 > dis12:
return distance((lat2, lon2), (lat3, lon3)).meters
return math.fabs(dxt)
添加 wdickerson 实现的 ObjectiveC 翻译：
#define DEGREES_RADIANS(angle) ((angle) / 180.0 * M_PI)
#define RADIANS_DEGREES(angle) ((angle) / M_PI * 180)
- (double)crossArcFromCoord:(CLLocationCoordinate2D)fromCoord usingArcFromCoord:(CLLocationCoordinate2D)arcCoord1 toArcCoord:(CLLocationCoordinate2D)arcCoord2 {
fromCoord.latitude = DEGREES_RADIANS(fromCoord.latitude);
fromCoord.longitude = DEGREES_RADIANS(fromCoord.longitude);
arcCoord1.latitude = DEGREES_RADIANS(arcCoord1.latitude);
arcCoord1.longitude = DEGREES_RADIANS(arcCoord1.longitude);
arcCoord2.latitude = DEGREES_RADIANS(arcCoord2.latitude);
arcCoord2.longitude = DEGREES_RADIANS(arcCoord2.longitude);
double R = 6371000; // Earth's radius in meters
// Prerequisites for the formulas
double bear12 = [self bearFromCoord:arcCoord1 toCoord:arcCoord2];
double bear13 = [self bearFromCoord:arcCoord1 toCoord:fromCoord];
double dis13 = [self distFromCoord:arcCoord1 toCoord:fromCoord];
double diff = fabs(bear13 - bear12);
if (diff > M_PI) {
diff = 2 * M_PI - diff;
}
// Is relative bearing obtuse?
if (diff > (M_PI/2)) {
return dis13;
}
// Find the cross-track distance
double dxt = asin(sin(dis13 / R) * sin(bear13 - bear12)) * R;
// Is p4 beyond the arc?
double dis12 = [self distFromCoord:arcCoord1 toCoord:arcCoord2];
double dis14 = acos(cos(dis13 / R) / cos(dxt / R)) * R;
if (dis14 > dis12) {
return [self distFromCoord:arcCoord2 toCoord:fromCoord];
}
return fabs(dxt);
}
- (double)distFromCoord:(CLLocationCoordinate2D)coord1 toCoord:(CLLocationCoordinate2D)coord2 {
double R = 6371000;
return acos(sin(coord1.latitude) * sin(coord2.latitude) + cos(coord1.latitude) * cos(coord2.latitude) * cos(coord2.longitude - coord2.longitude)) * R;
}
- (double)bearFromCoord:(CLLocationCoordinate2D)fromCoord toCoord:(CLLocationCoordinate2D)toCoord {
return atan2(sin(toCoord.longitude - fromCoord.longitude) * cos(toCoord.latitude),
cos(fromCoord.latitude) * sin(toCoord.latitude) - (sin(fromCoord.latitude) * cos(toCoord.latitude) * cos(toCoord.longitude - fromCoord.longitude)));
}
添加 python+numpy 实现（现在您可以将经度和纬度作为数组传递并同时计算所有距离而无需循环）。
def _angularSeparation(long1, lat1, long2, lat2):
"""All radians
"""
t1 = np.sin(lat2/2.0 - lat1/2.0)**2
t2 = np.cos(lat1)*np.cos(lat2)*np.sin(long2/2.0 - long1/2.0)**2
_sum = t1 + t2
if np.size(_sum) == 1:
if _sum < 0.0:
_sum = 0.0
else:
_sum = np.where(_sum < 0.0, 0.0, _sum)
return 2.0*np.arcsin(np.sqrt(_sum))
def bear(latA, lonA, latB, lonB):
"""All radians
"""
# BEAR Finds the bearing from one lat / lon point to another.
result = np.arctan2(np.sin(lonB - lonA) * np.cos(latB),
np.cos(latA) * np.sin(latB) - np.sin(latA) * np.cos(latB) * np.cos(lonB - lonA)
)
return result
def pointToLineDistance(lon1, lat1, lon2, lat2, lon3, lat3):
"""All radians
points 1 and 2 define an arc segment,
this finds the distance of point 3 to the arc segment.
"""
result = lon1*0
needed = np.ones(result.size, dtype=bool)
bear12 = bear(lat1, lon1, lat2, lon2)
bear13 = bear(lat1, lon1, lat3, lon3)
dis13 = _angularSeparation(lon1, lat1, lon3, lat3)
# Is relative bearing obtuse?
diff = np.abs(bear13 - bear12)
if np.size(diff) == 1:
if diff > np.pi:
diff = 2*np.pi - diff
if diff > (np.pi / 2):
return dis13
else:
solved = np.where(diff > (np.pi / 2))[0]
result[solved] = dis13[solved]
needed[solved] = 0
# Find the cross-track distance.
dxt = np.arcsin(np.sin(dis13) * np.sin(bear13 - bear12))
# Is p4 beyond the arc?
dis12 = _angularSeparation(lon1, lat1, lon2, lat2)
dis14 = np.arccos(np.cos(dis13) / np.cos(dxt))
if np.size(dis14) == 1:
if dis14 > dis12:
return _angularSeparation(lon2, lat2, lon3, lat3)
else:
solved = np.where(dis14 > dis12)[0]
result[solved] = _angularSeparation(lon2[solved], lat2[solved], lon3[solved], lat3[solved])
if np.size(lon1) == 1:
return np.abs(dxt)
else:
result[needed] = np.abs(dxt[needed])
return result
我需要计算从 lat/lng GPS 点 P 到另外 2 个 lat/lng GPS 点 A 和 B 描述的线段的最短距离。
'Cross-track distance'帮我算出P到A、B所描述的大圆的最短距离
然而，这不是我想要的。我需要 P 和 A-B 的 线段 之间的距离，而不是整个大圆。
我使用了 http://www.movable-type.co.uk/scripts/latlong.html
中的以下实现Formula: dxt = asin( sin(δ13) ⋅ sin(θ13−θ12) ) ⋅ R
where:
δ13 is (angular) distance from start point to third point
θ13 is (initial) bearing from start point to third point
θ12 is (initial) bearing from start point to end point
R is the earth’s radius
以下图片有望展示我正在尝试解决的问题：
在第一张图片中，绿色线表示的交叉轨道距离是正确的，确实是到线段 AB 的最短距离。
在第二张图中显示了交叉轨道距离的问题，在这种情况下，我希望最短距离是简单距离 AP，但是交叉轨道距离给出了 [=34 指示的距离=]红色行。
如何更改我的算法以考虑到这一点，或检查点 X 是否在 AB 内。有可能通过计算来做到这一点吗？还是迭代是唯一可能的（昂贵的）解决方案？ （沿AB取N个点，计算P点到所有这些点的最小距离）
为了简单起见，图像中的所有线条都是直的。实际上，这些是大圆上的小圆弧
对于 100 - 1000m 的球形问题，很容易只转换为
笛卡尔坐标 space，使用等角投影。
然后继续学校数学：
使用函数 "distance from line segment" 很容易找到现成的实现。
此函数使用（有时 returns）直线 A、B 上的投影点 X 的相对 forward/backward 位置。值为
- 在区间 [0,1] 中，如果投影点在线段内。
- 如果X在A之前在外面，则为负数，
- 如果在B之后在外面则>1。
如果相对位置在 0,1 之间，则采用正常距离，如果在起点和线终点 A、B 的较短距离之外。
这种/或非常相似的笛卡尔实现的一个例子是Shortest distance between a point and a line segment
首先，一些术语：
我们的弧线是从 p1 到 p2 绘制的。
我们的第三点是p3。
与大圆相交的虚点是 p4。
p1 由 lat1,lon1 定义； p2 通过 lat2,lon2；等等
dis12 是 p1 到 p2 的距离；等等
bear12 是从 p1 到 p2 的方位角；等等
dxt 是交叉轨道距离。
dxa是跨弧距离，我们的目标！
注意交叉轨道公式依赖于相对方位，bear13-bear12
我们有3个案例要处理。
情况一：相对方位为钝角。所以，dxa=dis13.
案例2.1：相对方位是锐角，p4落在我们的弧上。 所以，dxa=dxt.
案例2.2：相对方位是锐角，p4超出了我们的范围。 所以，dxa=dis23
算法：
第一步： 如果相对方位角是钝角，dxa=dis13
完毕！
第 2 步： 如果相对方位角是锐角：
2.1: 查找 dxt。
2.3:寻找dis12。
2.4:寻找dis14。
2.4: 如果 dis14>dis12，dxa=dis23。
完毕！
2.5: 如果我们到达这里，dxa=abs(dxt)
MATLAB代码：
function [ dxa ] = crossarc( lat1,lon1,lat2,lon2,lat3,lon3 )
%// CROSSARC Calculates the shortest distance in meters
%// between an arc (defined by p1 and p2) and a third point, p3.
%// Input lat1,lon1,lat2,lon2,lat3,lon3 in degrees.
lat1=deg2rad(lat1); lat2=deg2rad(lat2); lat3=deg2rad(lat3);
lon1=deg2rad(lon1); lon2=deg2rad(lon2); lon3=deg2rad(lon3);
R=6371000; %// Earth's radius in meters
%// Prerequisites for the formulas
bear12 = bear(lat1,lon1,lat2,lon2);
bear13 = bear(lat1,lon1,lat3,lon3);
dis13 = dis(lat1,lon1,lat3,lon3);
diff = abs(bear13-bear12);
if diff > pi
diff = 2 * pi - diff;
end
%// Is relative bearing obtuse?
if diff>(pi/2)
dxa=dis13;
else
%// Find the cross-track distance.
dxt = asin( sin(dis13/R)* sin(bear13 - bear12) ) * R;
%// Is p4 beyond the arc?
dis12 = dis(lat1,lon1,lat2,lon2);
dis14 = acos( cos(dis13/R) / cos(dxt/R) ) * R;
if dis14>dis12
dxa=dis(lat2,lon2,lat3,lon3);
else
dxa=abs(dxt);
end
end
end
function [ d ] = dis( latA, lonA, latB, lonB )
%DIS Finds the distance between two lat/lon points.
R=6371000;
d = acos( sin(latA)*sin(latB) + cos(latA)*cos(latB)*cos(lonB-lonA) ) * R;
end
function [ b ] = bear( latA,lonA,latB,lonB )
%BEAR Finds the bearing from one lat/lon point to another.
b=atan2( sin(lonB-lonA)*cos(latB) , ...
cos(latA)*sin(latB) - sin(latA)*cos(latB)*cos(lonB-lonA) );
end
示例输出： 演示所有案例。请参阅下面的地图。
>> crossarc(-10.1,-55.5,-15.2,-45.1,-10.5,-62.5)
ans =
7.6709e+05
>> crossarc(40.5,60.5,50.5,80.5,51,69)
ans =
4.7961e+05
>> crossarc(21.72,35.61,23.65,40.7,25,42)
ans =
1.9971e+05
地图上那些相同的输出！：
演示案例 1：
演示案例 2.1：
演示案例 2.2：
归功于：http://www.movable-type.co.uk/scripts/latlong.html
对于公式
和：http://www.darrinward.com/lat-long/?id=1788764
用于生成地图图像。
/**
* Calculates the euclidean distance from a point to a line segment.
*
* @param v the point
* @param a start of line segment
* @param b end of line segment
* @return an array of 2 doubles:
* [0] distance from v to the closest point of line segment [a,b],
* [1] segment coeficient of the closest point of the segment.
* Coeficient values < 0 mean the closest point is a.
* Coeficient values > 1 mean the closest point is b.
* Coeficient values between 0 and 1 mean how far along the segment the closest point is.
*
* @author Afonso Santos
*/
public static
double[]
distanceToSegment( final R3 v, final R3 a, final R3 b )
{
double[] results = new double[2] ;
final R3 ab_ = b.sub( a ) ;
final double ab = ab_.modulus( ) ;
final R3 av_ = v.sub( a ) ;
final double av = av_.modulus( ) ;
if (ab == 0.0) // a and b coincide
{
results[0] = av ; // Distance
results[1] = 0.0 ; // Segment coeficient.
}
else
{
final double avScaProjAb = av_.dot(ab_) / ab ;
final double abCoeficient = results[1] = avScaProjAb / ab ;
if (abCoeficient <= 0.0) // Point is before start of the segment ?
results[0] = av ; // Use distance to start of segment.
else if (abCoeficient >= 1.0) // Point is past the end of the segment ?
results[0] = v.sub( b ).modulus() ; // Use distance to end of segment.
else // Point is within the segment's start/end perpendicular boundaries.
{
if (avScaProjAb >= av) // Test to avoid machine float representation epsilon rounding errors that would result in expection on sqrt.
results[0] = 0.0 ; // a, b and v are colinear.
else
results[0] = Math.sqrt( av * av - avScaProjAb * avScaProjAb ) ; // Perpendicular distance from point to segment.
}
}
return results ;
}
上述方法需要笛卡尔 3D space 个参数，而您要求使用 lat/lon 个参数。要进行转换，请使用
/**
* Calculate 3D vector (from center of earth).
*
* @param latDeg latitude (degrees)
* @param lonDeg longitude (degrees)
* @param eleMtr elevation (meters)
* @return 3D cartesian vector (from center of earth).
*
* @author Afonso Santos
*/
public static
R3
cartesian( final double latDeg, final double lonDeg, final double eleMtr )
{
return versor( latDeg, lonDeg ).scalar( EARTHMEANRADIUS_MTR + eleMtr ) ;
}
对于 3D/R3 代码的其余部分或如何计算到 path/route/track 检查的距离 https://sourceforge.net/projects/geokarambola/
将 Java 版本添加到 wdickerson 答案：
public static double pointToLineDistance(double lon1, double lat1, double lon2, double lat2, double lon3, double lat3) {
lat1 = Math.toRadians(lat1);
lat2 = Math.toRadians(lat2);
lat3 = Math.toRadians(lat3);
lon1 = Math.toRadians(lon1);
lon2 = Math.toRadians(lon2);
lon3 = Math.toRadians(lon3);
// Earth's radius in meters
double R = 6371000;
// Prerequisites for the formulas
double bear12 = bear(lat1, lon1, lat2, lon2);
double bear13 = bear(lat1, lon1, lat3, lon3);
double dis13 = dis(lat1, lon1, lat3, lon3);
// Is relative bearing obtuse?
if (Math.abs(bear13 - bear12) > (Math.PI / 2))
return dis13;
// Find the cross-track distance.
double dxt = Math.asin(Math.sin(dis13 / R) * Math.sin(bear13 - bear12)) * R;
// Is p4 beyond the arc?
double dis12 = dis(lat1, lon1, lat2, lon2);
double dis14 = Math.acos(Math.cos(dis13 / R) / Math.cos(dxt / R)) * R;
if (dis14 > dis12)
return dis(lat2, lon2, lat3, lon3);
return Math.abs(dxt);
}
private static double dis(double latA, double lonA, double latB, double lonB) {
double R = 6371000;
return Math.acos(Math.sin(latA) * Math.sin(latB) + Math.cos(latA) * Math.cos(latB) * Math.cos(lonB - lonA)) * R;
}
private static double bear(double latA, double lonA, double latB, double lonB) {
// BEAR Finds the bearing from one lat / lon point to another.
return Math.atan2(Math.sin(lonB - lonA) * Math.cos(latB), Math.cos(latA) * Math.sin(latB) - Math.sin(latA) * Math.cos(latB) * Math.cos(lonB - lonA));
}
并添加 python Sga 实现的翻译：
def bear(latA, lonA, latB, lonB):
# BEAR Finds the bearing from one lat / lon point to another.
return math.atan2(
math.sin(lonB - lonA) * math.cos(latB),
math.cos(latA) * math.sin(latB) - math.sin(latA) * math.cos(latB) * math.cos(lonB - lonA)
)
def pointToLineDistance(lon1, lat1, lon2, lat2, lon3, lat3):
lat1 = math.radians(lat1)
lat2 = math.radians(lat2)
lat3 = math.radians(lat3)
lon1 = math.radians(lon1)
lon2 = math.radians(lon2)
lon3 = math.radians(lon3)
R = 6378137
bear12 = bear(lat1, lon1, lat2, lon2)
bear13 = bear(lat1, lon1, lat3, lon3)
dis13 = distance( (lat1, lon1), (lat3, lon3)).meters
# Is relative bearing obtuse?
if math.fabs(bear13 - bear12) > (math.pi / 2):
return dis13
# Find the cross-track distance.
dxt = math.asin(math.sin(dis13 / R) * math.sin(bear13 - bear12)) * R
# Is p4 beyond the arc?
dis12 = distance((lat1, lon1), (lat2, lon2)).meters
dis14 = math.acos(math.cos(dis13 / R) / math.cos(dxt / R)) * R
if dis14 > dis12:
return distance((lat2, lon2), (lat3, lon3)).meters
return math.fabs(dxt)
添加 wdickerson 实现的 ObjectiveC 翻译：
#define DEGREES_RADIANS(angle) ((angle) / 180.0 * M_PI)
#define RADIANS_DEGREES(angle) ((angle) / M_PI * 180)
- (double)crossArcFromCoord:(CLLocationCoordinate2D)fromCoord usingArcFromCoord:(CLLocationCoordinate2D)arcCoord1 toArcCoord:(CLLocationCoordinate2D)arcCoord2 {
fromCoord.latitude = DEGREES_RADIANS(fromCoord.latitude);
fromCoord.longitude = DEGREES_RADIANS(fromCoord.longitude);
arcCoord1.latitude = DEGREES_RADIANS(arcCoord1.latitude);
arcCoord1.longitude = DEGREES_RADIANS(arcCoord1.longitude);
arcCoord2.latitude = DEGREES_RADIANS(arcCoord2.latitude);
arcCoord2.longitude = DEGREES_RADIANS(arcCoord2.longitude);
double R = 6371000; // Earth's radius in meters
// Prerequisites for the formulas
double bear12 = [self bearFromCoord:arcCoord1 toCoord:arcCoord2];
double bear13 = [self bearFromCoord:arcCoord1 toCoord:fromCoord];
double dis13 = [self distFromCoord:arcCoord1 toCoord:fromCoord];
double diff = fabs(bear13 - bear12);
if (diff > M_PI) {
diff = 2 * M_PI - diff;
}
// Is relative bearing obtuse?
if (diff > (M_PI/2)) {
return dis13;
}
// Find the cross-track distance
double dxt = asin(sin(dis13 / R) * sin(bear13 - bear12)) * R;
// Is p4 beyond the arc?
double dis12 = [self distFromCoord:arcCoord1 toCoord:arcCoord2];
double dis14 = acos(cos(dis13 / R) / cos(dxt / R)) * R;
if (dis14 > dis12) {
return [self distFromCoord:arcCoord2 toCoord:fromCoord];
}
return fabs(dxt);
}
- (double)distFromCoord:(CLLocationCoordinate2D)coord1 toCoord:(CLLocationCoordinate2D)coord2 {
double R = 6371000;
return acos(sin(coord1.latitude) * sin(coord2.latitude) + cos(coord1.latitude) * cos(coord2.latitude) * cos(coord2.longitude - coord2.longitude)) * R;
}
- (double)bearFromCoord:(CLLocationCoordinate2D)fromCoord toCoord:(CLLocationCoordinate2D)toCoord {
return atan2(sin(toCoord.longitude - fromCoord.longitude) * cos(toCoord.latitude),
cos(fromCoord.latitude) * sin(toCoord.latitude) - (sin(fromCoord.latitude) * cos(toCoord.latitude) * cos(toCoord.longitude - fromCoord.longitude)));
}
添加 python+numpy 实现（现在您可以将经度和纬度作为数组传递并同时计算所有距离而无需循环）。
def _angularSeparation(long1, lat1, long2, lat2):
"""All radians
"""
t1 = np.sin(lat2/2.0 - lat1/2.0)**2
t2 = np.cos(lat1)*np.cos(lat2)*np.sin(long2/2.0 - long1/2.0)**2
_sum = t1 + t2
if np.size(_sum) == 1:
if _sum < 0.0:
_sum = 0.0
else:
_sum = np.where(_sum < 0.0, 0.0, _sum)
return 2.0*np.arcsin(np.sqrt(_sum))
def bear(latA, lonA, latB, lonB):
"""All radians
"""
# BEAR Finds the bearing from one lat / lon point to another.
result = np.arctan2(np.sin(lonB - lonA) * np.cos(latB),
np.cos(latA) * np.sin(latB) - np.sin(latA) * np.cos(latB) * np.cos(lonB - lonA)
)
return result
def pointToLineDistance(lon1, lat1, lon2, lat2, lon3, lat3):
"""All radians
points 1 and 2 define an arc segment,
this finds the distance of point 3 to the arc segment.
"""
result = lon1*0
needed = np.ones(result.size, dtype=bool)
bear12 = bear(lat1, lon1, lat2, lon2)
bear13 = bear(lat1, lon1, lat3, lon3)
dis13 = _angularSeparation(lon1, lat1, lon3, lat3)
# Is relative bearing obtuse?
diff = np.abs(bear13 - bear12)
if np.size(diff) == 1:
if diff > np.pi:
diff = 2*np.pi - diff
if diff > (np.pi / 2):
return dis13
else:
solved = np.where(diff > (np.pi / 2))[0]
result[solved] = dis13[solved]
needed[solved] = 0
# Find the cross-track distance.
dxt = np.arcsin(np.sin(dis13) * np.sin(bear13 - bear12))
# Is p4 beyond the arc?
dis12 = _angularSeparation(lon1, lat1, lon2, lat2)
dis14 = np.arccos(np.cos(dis13) / np.cos(dxt))
if np.size(dis14) == 1:
if dis14 > dis12:
return _angularSeparation(lon2, lat2, lon3, lat3)
else:
solved = np.where(dis14 > dis12)[0]
result[solved] = _angularSeparation(lon2[solved], lat2[solved], lon3[solved], lat3[solved])
if np.size(lon1) == 1:
return np.abs(dxt)
else:
result[needed] = np.abs(dxt[needed])
return result