14

これはこの質問のフォローアップです。

私はこれにこだわっているようです。基本的に、標準的な度数システムの座標を参照するように、または国際日付変更線に沿って南極から北への距離を測定し、次に日付のその時点から始まる東の距離を測定するために、前後に変換できる必要があります。ライン。これを行うために (およびいくつかのより一般的な距離測定のものと同様に)、2 つの緯度/経度ポイント間の距離を決定する 1 つの方法と、緯度/経度ポイント、方位、および距離を取得して返す別の方法があります。そのコースの終わりの緯度/経度ポイント。

私が定義した2つの静的メソッドは次のとおりです。

/* Takes two lon/lat pairs and returns the distance between them in kilometers.
*/
public static double distance (double lat1, double lon1, double lat2, double lon2) {
    double theta = toRadians(lon1-lon2);
    lat1 = toRadians(lat1);
    lon1 = toRadians(lon1);
    lat2 = toRadians(lat2);
    lon2 = toRadians(lon2);

    double dist = sin(lat1)*sin(lat2) + cos(lat1)*cos(lat2)*cos(theta);
    dist = toDegrees(acos(dist)) * 60 * 1.1515 * 1.609344 * 1000;

    return dist;
}

/* endOfCourse takes a lat/lon pair, a heading (in degrees clockwise from north), and a distance (in kilometers), and returns
 * the lat/lon pair that would be reached by traveling that distance in that direction from the given point.
 */
public static double[] endOfCourse (double lat1, double lon1, double tc, double dist) {
    double pi = Math.PI;
    lat1 = toRadians(lat1);
    lon1 = toRadians(lon1);
    tc = toRadians(tc);
    double dist_radians = toRadians(dist / (60 * 1.1515 * 1.609344 * 1000));
    double lat = asin(sin(lat1) * cos(dist_radians) + cos(lat1) * sin(dist_radians) * cos(tc));
    double dlon = atan2(sin(tc) * sin(dist_radians) * cos(lat1), cos(dist_radians) - sin(lat1) * sin(lat));
    double lon = ((lon1-dlon + pi) % (2*pi)) - pi;
    double[] endPoint = new double[2];
    endPoint[0] = lat; endPoint[1] = lon;
    return endPoint;
}

そして、これが私がそれをテストするために使用している関数です:

public static void main(String args[]) throws java.io.IOException, java.io.FileNotFoundException {
    double distNorth = distance(0.0, 0.0, 72.0, 0.0);
    double distEast = distance(72.0, 0.0, 72.0, 31.5);
    double lat1 = endOfCourse(0.0, 0.0, 0.0, distNorth)[0];
    double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[1];
    System.out.println("end at: " + lat1 + " / " + lon1);
    return;
}

「終了」値は appx である必要があります。72.0 / 31.5。しかし、代わりに約 1.25 / 0.021 を得ています。

どこかで単位を変換するのを忘れたり、何か愚かなことを見逃しているに違いないと思います...どんな助けも大歓迎です!

更新 1:

メートルを返す距離関数を (正しく) 書きましたが、コメントに間違ってキロメートルを書いてしまいました... もちろん、今日戻ってきたときは混乱しました。とにかく、これで修正されました。endOfCourse メソッドの因数分解エラーを修正しました。また、そのメソッドでもラジアンから度数に戻すのを忘れていたことに気付きました。とにかく: 正しい緯度番号 (71.99...) を取得しているように見えますが、経度番号はかなりずれています (11.5 ではなく 3.54 になります)。

更新 2: 以下で説明するように、テストでタイプミスがありました。現在はコードで修正されています。ただし、経度の数値はまだ間違っています。現在、11.5 ではなく -11.34 になっています。これらの行に何か問題があるに違いないと思います:

double dlon = atan2(sin(tc) * sin(dist_radians) * cos(lat1), cos(dist_radians) - sin(lat1) * sin(lat));
double lon = ((lon1-dlon + pi) % (2*pi)) - pi;
4

5 に答える 5

61

コード内のマジック ナンバーの重大なケースがあります。表現:

 (60 * 1.1515 * 1.609344 * 1000)

2回出てきますが、あまり説明がありません。いくつかの助けを借りて:1.609344は1マイルのキロメートル数です。60 は 1 度の分数です。1000 は 1 キロメートルのメートル数です。1.1515 は海里の法定マイル数です (ありがとう、DanM)。1 海里は、赤道における 1 分間の緯度の長さです。

球状の地球ではなく、球状の地球モデルを使用していると思いますか? 代数は、回転楕円体になるほど複雑ではありません。

最初の公式 - 2 つの緯度と経度のペア間の変換 - は奇数です。答えを整理するには、デルタ緯度 (Δλ) とデルタ経度 (Δφ) の両方が必要です。さらに、ペア間の距離:

(60° N, 30° W), (60° N, 60° W)
(60° N, 60° W), (60° N, 90° W)

同じはずですが、コードが異なる答えを生成すると確信しています。

ですから、球面三角法の参考資料に戻って、何が間違っているのかを確認する必要があると思います。(この件に関する私の本を見つけるのにしばらく時間がかかります。どの箱に入っているかに関係なく、開梱する必要があります。)

[ ...時間が経ち...開梱完了... ]

頂点と辺abcで角度ABCを持つ球面三角形(つまり、辺aBからCなど) を考えると、コサイン公式は次のようになります。

cos a = cos b . cos c + sin b . sin c . cos A

これを問題に当てはめると、与えられた 2 点をBCと呼び、 Aで直角をもつ直角球面三角形を作成します。

最悪の ASCII アート:

                  + C
                 /|
                / |
            a  /  | b
              /   |
             /    |
            /     |
         B +------+ A
              c

cは経度の差に等しくなります。辺bは緯度の差に等しい。角度Aは 90° なので、cos A = 0 です。したがって、aの式は次のようになります。

cos a = cos Δλ . cos Δφ + sin Δλ . sin Δφ . cos 90°

a = arccos (cos Δλ . cos Δφ)

ラジアン単位の角度aは、地球の半径を掛けて距離に変換されます。別の方法として、aの度数 (および度数の分数) を指定すると、1 度は 60 海里であるため、60 * 1.1515 法定マイル、1 度は 60 * 1.1515 * 1.609344 キロメートルになります。メートル単位の距離が必要でない限り、係数 1000 は必要ないと思います。

Paul Tomblinは、方程式のソースとしてAviation Formulary v1.44を指摘しています。実際、位置の違いが小さい場合の数値的に安定したバージョンとともに、そこにあります。

基本的な三角法に行くと、次のこともわかっています。

cos (A - B) = cos A . cos B + sin A . sin B

私が与えた式にそれを 2 回適用すると、Aviation Formulary の式になる可能性があります。

(私の参考文献: AE Roy and D Clarke (2003) による「Astronomy: Principles and Practice, Fourth Edition」。私のコピーは、1977 年の初版、Adam Hilger、ISBN 0-85274-346-7 です。)


NBチェックアウト (Google) 'define:"海里"'; 海里は現在、定義により 1852 m (1.852 km) になっているようです。乗数 1.1515 は、海里の古い定義である約 6080 フィートに対応します。bcスケール 10 を使用すると、次のようになります。

(1852/(3*0.3048))/1760
1.1507794480

どの要因が効果的かは、基礎が何であるかによって異なります。


第 1 原理からの 2 番目の問題を見ると、少し異なる設定があり、「もう 1 つの」球面三角法である正弦公式が必要です。

sin A   sin B   sin C
----- = ----- = -----
sin a   sin b   sin c

前の図を適応させる:

                  + C
                 /|
                / |
            a  /  | b
           |  /   |
           |X/    |
           |/     |
         B +------+ A
              c

開始点B、角度X = 90° - B、長さ (角度) a、角度A = 90° が与えられます。あなたが求めているのは、 b (緯度のデルタ) とc (経度のデルタ) です。

したがって、次のようになります。

sin a   sin b
----- = ----
sin A   sin B

または

        sin a . sin B
sin b = -------------
            sin A

または、A = 90°であるため、sin A = 1、sin B = sin (90° - X) = cos X:

sin b = sin a . cos X

つまり、移動距離を角度aに変換し、その正弦を取り、コース方向の余弦を掛けて、結果の逆正弦をとります。

ab (計算されたばかり)、およびABが与えられると、コサイン式を適用してcを取得できます。Cの値がないため、単純にサイン公式を再適用してcを取得することはできないことに注意してください。また、球面三角法で遊んでいるため、C = 90° - B (和球面三角形の角の角度は 180° よりも大きくなる可能性があります; すべての角度が 90° に等しい正球面三角形を考えてみてください。これは完全に実現可能です)。


于 2008-12-23T15:58:58.137 に答える
6

http://www.movable-type.co.uk/scripts/latlong.htmlをご覧ください

そのサイトには、役立つはずのさまざまな式と Javascript コードがたくさんあります。私はそれを C# と SQL Server UDF の両方に正常に変換し、あらゆる場所で使用しています。

たとえば、距離を計算する Javascript の場合:

var R = 6371; // km
var φ1 = lat1.toRadians();
var φ2 = lat2.toRadians();
var Δφ = (lat2-lat1).toRadians();
var Δλ = (lon2-lon1).toRadians();

var a = Math.sin(Δφ/2) * Math.sin(Δφ/2) +
        Math.cos(φ1) * Math.cos(φ2) *
        Math.sin(Δλ/2) * Math.sin(Δλ/2);
var c = 2 * Math.atan2(Math.sqrt(a), Math.sqrt(1-a));

var d = R * c; 

楽しみ!

于 2008-12-23T15:59:47.807 に答える
2

km とラジアンの変換が間違っています。海里は 1/60 度なので、1.15... がマイルから海里への変換であり、1.6... が km から法定マイルへの変換であると仮定すると、

   nm = km /  (1.1515 * 1.609344);
   deg = nm / 60;
   rad = toRadians(deg);

つまり、1000 倍もずれていると思います。

于 2008-12-23T16:07:03.533 に答える
1

あなたの更新された質問について:すべきではありません

double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[0];

なれ

double lon1 = endOfCourse(lat1, 0.0, 90.0, distEast)[1];
于 2008-12-23T17:33:54.810 に答える
0

他の回答と更新で言及されている実装エラーは別として、これらの式の大きな問題を理解しました。

大きな問題はこれでした: Distance メソッド (2 点間の距離を計算するため) は、大円距離を計算していました。もちろん、これは理にかなっています。これが 2 点間の最短経路です。 ただし、同じ平行線 (緯度線) 上にある 2 点間の大圏距離は、赤道にいる場合を除き、緯度線に沿って直接移動した場合の 2 点間の距離と同じではありません。

したがって、機能は正しく機能しています。ただし、元の質問で私が提案した代替座標系では、IDL に沿った北の距離と、結果として得られる緯度での平行線に沿った東の距離のみを見る必要があります。また、特定の緯線に沿った距離の計算は、大円に沿った距離の計算とはまったく異なります。

とにかく、そこにあります。

于 2008-12-26T20:55:59.580 に答える