2016-03-29 11 views
3

90°0′0″N 0°0′0″Eという形式の2つの座標ペアを文字列として取得し、半径R = 6371kmの球上のそれらの点間の距離を計算したい。地球儀上の2つの座標間の距離を計算する

「haversine」と「余弦の球法則」の2つの式がインターネットで見つかったが、それらは機能していないようだ。 2*pi*R/4を返す90°の角度の場合、haversineは正しく動作しますが、コサインは失敗し、0を返します。ランダムな座標が異なる別のポイントは、両方のアルゴリズムでfalse値を返します:haversineが高すぎ、コサインが低すぎます。

私の実装は間違っていますか、または誤ったアルゴリズムを選択しましたか?

代わりに、これらの計算をどのように行う必要がありますか?

(そして、はい、私はN/SおよびE/Wをチェックしていないよ、まだ、しかし、テストした座標はすべて北東半球にあることを知っている。)

ここに私はPython 3のコードです:ここで

import math, re 
R = 6371 
PAT = r'(\d+)°(\d+)′(\d+)″([NSEW])' 

def distance(first, second): 
    def coords_to_rads(s): 
     return [math.radians(int(d) +int(m)/60 +int(s)/3600) \ 
       for d, m, s, nswe in re.findall(PAT, s)] 

    y1, x1 = coords_to_rads(first) 
    y2, x2 = coords_to_rads(second) 
    dx = x1 - x2 
    dy = y1 - y2 

    print("coord string:", first, "|", second) 
    print("coord radians:", y1, x1, "|", y2, x2) 
    print("x/y-distances:", dy, dx) 

    a = math.sin(dx/2)**2 + math.cos(x1) * math.cos(x2) * math.sin(dy/2)**2 
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a)) 
    haversine = R * c 

    law_of_cosines = math.acos(math.sin(x1) * math.sin(x2) + \ 
           math.cos(x1) * math.cos(x2)) * R 

    print("HS:", round(haversine, 2), "LOC:", round(law_of_cosines, 2)) 

    return haversine 
    #return law_of_cosines 

if __name__ == '__main__': 
    def test(result, correct): 
     print("result: ", result) 
     print("correct:", correct) 

    test(distance("90°0′0″N 0°0′0″E", "0°0′0″N, 0°0′0″E"), 10007.5) 
    test(distance("51°28′48″N 0°0′0″E", "46°12′0″N, 6°9′0″E"), 739.2) 
    test(distance("90°0′0″N 0°0′0″E", "90°0′0″S, 0°0′0″W"), 20015.1) 
    test(distance("33°51′31″S, 151°12′51″E", "40°46′22″N 73°59′3″W"), 15990.2) 

は、いくつかの出力です:

coord string: 90°0′0″N 0°0′0″E | 0°0′0″N, 0°0′0″E 
coord radians: 1.5707963267948966 0.0 | 0.0 0.0 
x/y-distances: 1.5707963267948966 0.0 
HS: 10007.54 LOC: 0.0 
result: 10007.543398010286 
correct: 10007.5 

coord string: 51°28′48″N 0°0′0″E | 46°12′0″N, 6°9′0″E 
coord radians: 0.8984954989266809 0.0 | 0.8063421144213803 0.10733774899765128 
x/y-distances: 0.09215338450530064 -0.10733774899765128 
HS: 900.57 LOC: 683.85 
result: 900.5669567853056 
correct: 739.2 
+0

良い質問を、それが[codereview.SE]に行くべきかしら?私はあなたがすでに作業している、またはほとんど働いているコードサンプルをここで話題にしていると考えられている多くの質問があることを知っていますが、サイト間の正確な「カットオフ」についてはわかりません。 –

+0

@DavidZ私が理解しているように、CRは最適化を探す作業用のコードスニペットを必要としています。私はここにそれを掲載しました。しかし、それはSOからCRへ、または戻ってマイグレーションされる私の最初のものではないでしょう... –

+0

私には意味があります。とにかく、[Wikipedia](https://en.wikipedia.org/wiki/Haversine_formula)にリンクしたサイトをチェックしたところ、haversine式の定義に矛盾があるようです不一致です)。それを調べましたか? –

答えて

1

あなたがaのあなたの計算にxyを混同ように見えます。あなたは経度()ではなく、緯度余弦(y)を取ることになっています。

私は(すなわちRで増殖しない)angular_distanceにごdistanceを変更し、いくつかの追加試験追加することにより、これを発見:

test(angular_distance("90°0′0″N 0°0′0″E", "89°0′0″N, 0°0′0″E"), math.radians(1)) 
test(angular_distance("90°0′0″N 0°0′0″E", "80°0′0″N, 0°0′0″E"), math.radians(10)) 
test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 0°0′0″E"), math.radians(40)) 
test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 20°0′0″E"), math.radians(40)) 
+0

ありがとうございます。私は単純に 'y1、x1 = coords_to_rads(最初) 'と' y2、x2 = coords_to_rads(秒)'という行を変更し、 'x'と' y'を入れ替えました。 –

関連する問題