地球上に線分(大円部分)があります。線分は、その両端の座標によって定義されます。明らかに、2 つの点は 2 つの線分を定義するので、短い線分に関心があるとします。
3 番目の点が与えられ、線と点の間の (最短) 距離を探しています。
すべての座標は経度\緯度 (WGS 84) で指定されます。
距離の計算方法を教えてください。
合理的なプログラミング言語でのソリューションで十分です。
地球上に線分(大円部分)があります。線分は、その両端の座標によって定義されます。明らかに、2 つの点は 2 つの線分を定義するので、短い線分に関心があるとします。
3 番目の点が与えられ、線と点の間の (最短) 距離を探しています。
すべての座標は経度\緯度 (WGS 84) で指定されます。
距離の計算方法を教えてください。
合理的なプログラミング言語でのソリューションで十分です。
ask Dr. Mathのアイデアに基づいた、私自身のソリューションを次に示します。フィードバックをいただければ幸いです。
最初に免責事項。このソリューションは、球に対して正しいです。地球は球体ではなく、座標系 (WGS 84) は地球が球体であるとは想定していません。したがって、これは単なる概算であり、誤差を実際に見積もることはできません。また、距離が非常に短い場合は、すべてが同一平面上にあると仮定することで、適切な近似を取得することもおそらく可能です。繰り返しになりますが、距離をどれだけ「小さく」する必要があるかわかりません。
さぁ、ビジネスへ。線分 A、B の端点、および 3 番目の点 C を呼び出します。基本的に、アルゴリズムは次のようになります。
次の 3 つのベクトル積を使用して、C に最も近い直線 AB 上の点 T を計算します。
G = A×B
F = C×G
T = G×F
T を正規化し、地球の半径を掛けます。
C と、A と B によって定義される大円との間の距離を探している場合、これらの手順で十分です。私のように、C と短い方の線分との間の距離に関心がある場合は、それを確認する追加の手順を実行する必要があります。 T は確かにこのセグメントにあります。そうでない場合は、必然的に最も近い点は端 A または B のいずれかになります。最も簡単な方法は、どちらかを確認することです。
一般的に、3 つのベクトル積の背後にある考え方は次のとおりです。最初のもの (G) は、A と B の大円の平面 (つまり、A、B と原点を含む平面) を示します。2 番目の (F) は、C を通り、G に垂直な大円を示します。次に、T は、F と G によって定義される大円の交点であり、正規化と R の乗算によって正しい長さにされます。
これを行うための Java コードの一部を次に示します。
大円上の最も近い点を見つける. 入力と出力は長さ 2 の配列です。中間配列の長さは 3 です。
double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
normalize(t);
multiplyByScalar(t, R_EARTH);
return fromCartsian(t);
}
セグメント上の最も近い点を見つける:
double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (distance(a,c) < distance(b,c)) ? a : c;
}
これは、A および B と同じ大円上にあることがわかっている点 T が、この大円の短いセグメント上にあるかどうかをテストする簡単な方法です。ただし、それを行うためのより効率的な方法があります。
boolean onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.abs(distance(a,b)-distance(a,t)-distance(b,t)) < PRECISION;
}
Ask Dr. Mathから、点から大円までの距離を試してください。経度/緯度を球座標に変換し、地球の半径に合わせてスケーリングする必要がありますが、これは良い方向のようです。
誰かがそれを必要とする場合、これは c# に移植された loleksy の回答です。
private static double _eQuatorialEarthRadius = 6378.1370D;
private static double _d2r = (Math.PI / 180D);
private static double PRECISION = 0.1;
// Haversine Algorithm
// source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates
private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
return (1000D * HaversineInKM(lat1, long1, lat2, long2));
}
private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.Pow(Math.Sin(dlat / 2D), 2D) + Math.Cos(lat1 * _d2r) * Math.Cos(lat2 * _d2r)
* Math.Pow(Math.Sin(dlong / 2D), 2D);
double c = 2D * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;
return d;
}
// Distance between a point and a line
static double pointLineDistanceGEO(double[] a, double[] b, double[] c)
{
double[] nearestNode = nearestPointGreatCircle(a, b, c);
double result = HaversineInKM(c[0], c[1], nearestNode[0], nearestNode[1]);
return result;
}
// source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
private static double[] nearestPointGreatCircle(double[] a, double[] b, double [] c)
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}
private static bool onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.Abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
}
// source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
private static double[] toCartsian(double[] coord) {
double[] result = new double[3];
result[0] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Cos(deg2rad(coord[1]));
result[1] = _eQuatorialEarthRadius * Math.Cos(deg2rad(coord[0])) * Math.Sin(deg2rad(coord[1]));
result[2] = _eQuatorialEarthRadius * Math.Sin(deg2rad(coord[0]));
return result;
}
private static double[] fromCartsian(double[] coord){
double[] result = new double[2];
result[0] = rad2deg(Math.Asin(coord[2] / _eQuatorialEarthRadius));
result[1] = rad2deg(Math.Atan2(coord[1], coord[0]));
return result;
}
// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
double[] result = new double[3];
result[0] = a[1] * b[2] - a[2] * b[1];
result[1] = a[2] * b[0] - a[0] * b[2];
result[2] = a[0] * b[1] - a[1] * b[0];
return result;
}
private static double[] normalize(double[] t) {
double length = Math.Sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
double[] result = new double[3];
result[0] = t[0]/length;
result[1] = t[1]/length;
result[2] = t[2]/length;
return result;
}
private static double[] multiplyByScalar(double[] normalize, double k) {
double[] result = new double[3];
result[0] = normalize[0]*k;
result[1] = normalize[1]*k;
result[2] = normalize[2]*k;
return result;
}
私は基本的に今同じことを探していますが、厳密に言えば、大円のセグメントを持つことは気にせず、完全な円上の任意の点までの距離が必要なだけです。
現在調査中の 2 つのリンク:
このページでは「クロストラック距離」について言及していますが、これは基本的にあなたが探しているものと思われます。
また、PostGIS メーリング リストの次のスレッドでは、(1) 2D 平面上の線距離に使用されるのと同じ式 (PostGIS の line_locate_point を使用) を使用して大円上の最も近い点を決定し、次に(2) その点と回転楕円体上の 3 番目の点の間の距離を計算します。数学的にステップ (1) が正しいかどうかはわかりませんが、驚くでしょう。
http://postgis.refractions.net/pipermail/postgis-users/2009-July/023903.html
最後に、「関連」の下に次のリンクが表示されていることを確認しました。
これは、イデオンフィドルとして受け入れられた回答の完全なコードです(ここにあります):
import java.util.*;
import java.lang.*;
import java.io.*;
/* Name of the class has to be "Main" only if the class is public. */
class Ideone
{
private static final double _eQuatorialEarthRadius = 6378.1370D;
private static final double _d2r = (Math.PI / 180D);
private static double PRECISION = 0.1;
// Haversine Algorithm
// source: http://stackoverflow.com/questions/365826/calculate-distance-between-2-gps-coordinates
private static double HaversineInM(double lat1, double long1, double lat2, double long2) {
return (1000D * HaversineInKM(lat1, long1, lat2, long2));
}
private static double HaversineInKM(double lat1, double long1, double lat2, double long2) {
double dlong = (long2 - long1) * _d2r;
double dlat = (lat2 - lat1) * _d2r;
double a = Math.pow(Math.sin(dlat / 2D), 2D) + Math.cos(lat1 * _d2r) * Math.cos(lat2 * _d2r)
* Math.pow(Math.sin(dlong / 2D), 2D);
double c = 2D * Math.atan2(Math.sqrt(a), Math.sqrt(1D - a));
double d = _eQuatorialEarthRadius * c;
return d;
}
// Distance between a point and a line
public static void pointLineDistanceTest() {
//line
//double [] a = {50.174315,19.054743};
//double [] b = {50.176019,19.065042};
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
//double [] c = {50.184373,19.054657};
double [] c = {52.008308, 17.542927};
double[] nearestNode = nearestPointGreatCircle(a, b, c);
System.out.println("nearest node: " + Double.toString(nearestNode[0]) + "," + Double.toString(nearestNode[1]));
double result = HaversineInM(c[0], c[1], nearestNode[0], nearestNode[1]);
System.out.println("result: " + Double.toString(result));
}
// source: http://stackoverflow.com/questions/1299567/how-to-calculate-distance-from-a-point-to-a-line-segment-on-a-sphere
private static double[] nearestPointGreatCircle(double[] a, double[] b, double c[])
{
double[] a_ = toCartsian(a);
double[] b_ = toCartsian(b);
double[] c_ = toCartsian(c);
double[] G = vectorProduct(a_, b_);
double[] F = vectorProduct(c_, G);
double[] t = vectorProduct(G, F);
return fromCartsian(multiplyByScalar(normalize(t), _eQuatorialEarthRadius));
}
@SuppressWarnings("unused")
private static double[] nearestPointSegment (double[] a, double[] b, double[] c)
{
double[] t= nearestPointGreatCircle(a,b,c);
if (onSegment(a,b,t))
return t;
return (HaversineInKM(a[0], a[1], c[0], c[1]) < HaversineInKM(b[0], b[1], c[0], c[1])) ? a : b;
}
private static boolean onSegment (double[] a, double[] b, double[] t)
{
// should be return distance(a,t)+distance(b,t)==distance(a,b),
// but due to rounding errors, we use:
return Math.abs(HaversineInKM(a[0], a[1], b[0], b[1])-HaversineInKM(a[0], a[1], t[0], t[1])-HaversineInKM(b[0], b[1], t[0], t[1])) < PRECISION;
}
// source: http://stackoverflow.com/questions/1185408/converting-from-longitude-latitude-to-cartesian-coordinates
private static double[] toCartsian(double[] coord) {
double[] result = new double[3];
result[0] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.cos(Math.toRadians(coord[1]));
result[1] = _eQuatorialEarthRadius * Math.cos(Math.toRadians(coord[0])) * Math.sin(Math.toRadians(coord[1]));
result[2] = _eQuatorialEarthRadius * Math.sin(Math.toRadians(coord[0]));
return result;
}
private static double[] fromCartsian(double[] coord){
double[] result = new double[2];
result[0] = Math.toDegrees(Math.asin(coord[2] / _eQuatorialEarthRadius));
result[1] = Math.toDegrees(Math.atan2(coord[1], coord[0]));
return result;
}
// Basic functions
private static double[] vectorProduct (double[] a, double[] b){
double[] result = new double[3];
result[0] = a[1] * b[2] - a[2] * b[1];
result[1] = a[2] * b[0] - a[0] * b[2];
result[2] = a[0] * b[1] - a[1] * b[0];
return result;
}
private static double[] normalize(double[] t) {
double length = Math.sqrt((t[0] * t[0]) + (t[1] * t[1]) + (t[2] * t[2]));
double[] result = new double[3];
result[0] = t[0]/length;
result[1] = t[1]/length;
result[2] = t[2]/length;
return result;
}
private static double[] multiplyByScalar(double[] normalize, double k) {
double[] result = new double[3];
result[0] = normalize[0]*k;
result[1] = normalize[1]*k;
result[2] = normalize[2]*k;
return result;
}
public static void main(String []args){
System.out.println("Hello World");
Ideone.pointLineDistanceTest();
}
}
コメント付きのデータでは問題なく機能します。
//line
double [] a = {50.174315,19.054743};
double [] b = {50.176019,19.065042};
//point
double [] c = {50.184373,19.054657};
最寄りのノード: 50.17493121381319,19.05846668493702
しかし、私はこのデータに問題があります:
double [] a = {52.00118, 17.53933};
double [] b = {52.00278, 17.54008};
//point
double [] c = {52.008308, 17.542927};
最も近いノードは次のとおりです: 52.00834987257176,17.542691313436357 これは間違っています。
2点で指定された線分は閉じた線分ではないと思います。
数千メートルまでの距離については、球から平面への問題を単純化します。次に、簡単な三角形の計算を使用できるため、問題は非常に単純です。
点 A と点 B があり、直線 AB までの距離 X を探します。それで:
Location a;
Location b;
Location x;
double ax = a.distanceTo(x);
double alfa = (Math.abs(a.bearingTo(b) - a.bearingTo(x))) / 180
* Math.PI;
double distance = Math.sin(alfa) * ax;
球上の2点間の最短距離は、2点を通過する大円の小さい方の側です。あなたはすでにこれを知っていると確信しています。ここhttp://www.physicsforums.com/archive/index.php/t-178252.htmlにも同様の質問があり、数学的にモデル化するのに役立つ場合があります。
正直なところ、これのコード化された例を取得する可能性がどれほどあるかはわかりません。