[pukiwiki]
Ａ地点から～B地点まで～の距離を計算する関数 geodistanceを書いてみますた。
2地点の緯度経度から距離を概算します

誤差が500メートルぐらいありますが、（今必要な値としては）1km単位ぐらいまで許容範囲だったので、これで行くことに。　メモ
-[[測量計算:http://vldb.gsi.go.jp/sokuchi/surveycalc/]]
-[[距離と方位角の計算:http://vldb.gsi.go.jp/sokuchi/surveycalc/algorithm/]]
こちらに、ちゃんとしたアルゴリズムがあります。（多分、ミリ単位の精度）
でも、難しくてよくわからなかったので（涙

[/pukiwiki]

以下　関数　＆　テストコード

# -*- coding: cp932 -*-
import math
from math import cos,sin,acos,sqrt
u"""
#ウィキペディア地球
#http://ja.wikipedia.org/wiki/%E5%9C%B0%E7%90%83
#赤道半径が 6378.137 km、極半径が 6356.752 km
#40,075 km 	地球の赤道の長さ
#緯度 latのときの緯線１度あたりの距離
内積による　2点間の距離概算

半径は2点の緯度で線形補間？
地球の中心から　北緯 n度　　東経 e度を指し示すベクトルを正規化すると

東経　0度がX軸プラスス
東経 90度がY軸プラス
北極方向がZ軸プラス　と定める

z=sin n
t=cos n
x=t*cos e - 0*sin e
y=t*sin e + 0*cos e

"""
r1=6378.137
r2=6356.752
r=(r1+r2)/2.
dr=r1-r2
degree=math.pi/180.
ratio=r2/r1
def geodistance(lat0,lon0,lat1,lon1):
    x0,y0,z0=nvec(lat0,lon0)
    x1,y1,z1=nvec(lat1,lon1)
    dot=x0*x1+y0*y1+z0*z1
    #dot= |A||B|cosθ
    t=acos(dot)
    r=(radius(lat0)+radius(lat1))/2.

    dist=t*r
    return dist
def radius(lat):  
    r=r1-dr*abs(lat/90.)

    return r
def nvec(lat,lon):
    lat2,lon2=lat*degree,lon*degree
    z=ratio*sin(lat2)
    t=cos(lat2)
    x=t*cos(lon2) #- 0*sin e
    y=t*sin(lon2) #+ 0*cos e
    l=sqrt(x**2.+y**2.+z**2.)
    x/=l
    y/=l
    z/=l
    return x,y,z
  
if __name__=="__main__":
    ep=1/1000.
    assert abs(geodistance(0,0,0,1.)*360-40075.01668557849)<ep
    assert abs(geodistance(0,0,0,90.)*4-40075.01668557849)<ep
    assert geodistance(0,0,1,0)<geodistance(89,0,90,0)
    assert 111<geodistance(89,0,90,0)<112

Mathematicaでテストデータを作って、子午線の長さと、緯度毎の緯線一度の長さをチェック。
斜め方向のチェックは　また今度。

子午線の長さで350メートル
緯線の長さで200メートルほどの誤差が出ます
補正は　また　時間があるときに。

斜めの距離は　また今度チェックします。

以下　テストコード

# -*- coding: cp932 -*-
#from pylab import *
from geodistance import geodistance
import random

def test1():
    ep=500.
    e=[]
    for t,res in data1:
        lat0,lat1=t
        lon=360.*(random.random()-0.5)
        ret=geodistance(lat0,lon,lat1,lon)*1000.
        err=ret-res
        e.append(err)
        #print lat0,lat1,lon ,res ,ret,res-ret
        assert abs(err)<ep
    #plot(e)
    #show()
def test2():
    ep=200.
    for lat,res in data2:
        lon0=360.*(random.random()-0.5)
        lon1=lon0+1.
        ret=geodistance(lat,lon0,lat,lon1)*1000.
        #print lat,lon0,lon1 ,res ,ret,res-ret
        assert abs(ret-res)<ep


        
#GeoDistance[(#, 136),(#,137)]&/@Table[x,(x,0, 90)]
data1=(((0, 1), 110574.38855415277), 
((1, 2), 110575.06481069342), 
((2, 3), 110576.41652051672), 
((3, 4), 110578.4420779865), 
((4, 5), 110581.13907684629), 
((5, 6), 110584.50431285569), 
((6, 7), 110588.53378730224), 
((7, 8), 110593.22271138613), 
((8, 9), 110598.56551147692), 
((9, 10), 110604.5558352348), 
((10, 11), 110611.18655859602), 
((11, 12), 110618.44979361467), 
((12, 13), 110626.33689715691), 
((13, 14), 110634.83848044065), 
((14, 15), 110643.94441941284), 
((15, 16), 110653.64386596132), 
((16, 17), 110663.92525994252), 
((17, 18), 110674.77634203236), 
((18, 19), 110686.18416737043), 
((19, 20), 110698.13512000655), 
((20, 21), 110710.61492812574), 
((21, 22), 110723.60868003967), 
((22, 23), 110737.1008409439), 
((23, 24), 110751.07527041009), 
((24, 25), 110765.51524061138), 
((25, 26), 110780.40345526232), 
((26, 27), 110795.72206925572), 
((27, 28), 110811.45270897615), 
((28, 29), 110827.57649328338), 
((29, 30), 110844.07405513317), 
((30, 31), 110860.92556382921), 
((31, 32), 110878.11074786755), 
((32, 33), 110895.60891837666), 
((33, 34), 110913.39899311108), 
((34, 35), 110931.45952097507), 
((35, 36), 110949.76870707203), 
((36, 37), 110968.30443823044), 
((37, 38), 110987.04430899528), 
((38, 39), 111005.9656480526), 
((39, 40), 111025.04554505984), 
((40, 41), 111044.26087785786), 
((41, 42), 111063.58834002037), 
((42, 43), 111083.00446873484), 
((43, 44), 111102.48567295846), 
((44, 45), 111122.00826183661), 
((45, 46), 111141.5484733303), 
((46, 47), 111161.08250305083), 
((47, 48), 111180.58653322102), 
((48, 49), 111200.03676178747), 
((49, 50), 111219.40943158424), 
((50, 51), 111238.68085955613), 
((51, 52), 111257.82746600527), 
((52, 53), 111276.82580377902), 
((53, 54), 111295.65258742096), 
((54, 55), 111314.28472220218), 
((55, 56), 111332.69933301568), 
((56, 57), 111350.8737930881), 
((57, 58), 111368.78575247644), 
((58, 59), 111386.41316628916), 
((59, 60), 111403.73432264249), 
((60, 61), 111420.72787023295), 
((61, 62), 111437.37284557625), 
((62, 63), 111453.64869980412), 
((63, 64), 111469.5353250161), 
((64, 65), 111485.01308013621), 
((65, 66), 111500.06281623749), 
((66, 67), 111514.66590130705), 
((67, 68), 111528.80424437404), 
((68, 69), 111542.46031905425), 
((69, 70), 111555.61718633534), 
((70, 71), 111568.2585167364), 
((71, 72), 111580.36861163845), 
((72, 73), 111591.93242389344), 
((73, 74), 111602.9355775882), 
((74, 75), 111613.36438697355), 
((75, 76), 111623.20587451702), 
((76, 77), 111632.4477880502), 
((77, 78), 111641.07861699462), 
((78, 79), 111649.08760759898), 
((79, 80), 111656.46477725181), 
((80, 81), 111663.20092771783), 
((81, 82), 111669.28765740496), 
((82, 83), 111674.71737255243), 
((83, 84), 111679.48329736333), 
((84, 85), 111683.57948305104), 
((85, 86), 111687.00081579809), 
((86, 87), 111689.74302357297), 
((87, 88), 111691.80268188249), 
((88, 89), 111693.17721832152), 
((89, 90), 111693.86491604138))


#緯度ごとの緯線一度あたりの長さ 赤道では赤道長とイコール
#GeoDistance[{#, 136},{#,137}]&/@Table[x,{x,0, 90}]

data2=(
  (0, 111319.49079327358), 
  (1, 111302.64933943351), 
  (2, 111252.12980021282), 
  (3, 111167.94664186916), 
  (4, 111050.12397269311),
  (5, 110898.69553979335), 
  (6, 110713.7047245823), 
  (7, 110495.20453694796), 
  (8, 110243.25760809594), 
  (9, 109957.93618203873), 
  (10, 109639.32210570878), 
  (11, 109287.50681766603), 
  (12, 108902.59133536776), 
  (13, 108484.68624096578), 
  (14, 108033.91166559128), 
  (15, 107550.39727208568), 
  (16, 107034.28223613203), 
  (17, 106485.71522573847), 
  (18, 105904.85437902318), 
  (19, 105291.86728024643), 
  (20, 104646.93093403339), 
  (21, 103970.23173772973), 
  (22, 103261.96545182857), 
  (23, 102522.33716840531), 
  (24, 101751.56127749757), 
  (25, 100949.86143136208), 
  (26, 100117.47050654185), 
  (27, 99254.63056367503), 
  (28, 98361.59280497493), 
  (29, 97438.61752931186), 
  (30, 96485.97408482508), 
  (31, 95503.9408189941), 
  (32, 94492.80502609798), 
  (33, 93452.8628919924), 
  (34, 92384.41943613302), 
  (35, 91287.78845077672), 
  (36, 90163.29243729111), 
  (37, 89011.26253950551), 
  (38, 87832.03847403749), 
  (39, 86625.96845753022), 
  (40, 85393.40913073995), 
  (41, 84134.72547941183), 
  (42, 82850.29075188849), 
  (43, 81540.48637339566), 
  (44, 80205.70185695408), 
  (45, 78846.33471086824), 
  (46, 77462.79034274847), 
  (47, 76055.48196002372), 
  (48, 74624.83046690872), 
  (49, 73171.26435779205), 
  (50, 71695.2196070158), 
  (51, 70197.13955502294), 
  (52, 68677.47479085259), 
  (53, 67136.68303096831), 
  (54, 65575.22899440922), 
  (55, 63993.584274260145), 
  (56, 62392.227205440766), 
  (57, 60771.64272882016), 
  (58, 59132.32225166838), 
  (59, 57474.7635044625), 
  (60, 55799.47039407), 
  (61, 54106.95285333881), 
  (62, 52397.72668712815), 
  (63, 50672.31341482132), 
  (64, 48931.24010936693), 
  (65, 47175.03923290086), 
  (66, 45404.24846900778), 
  (67, 43619.41055168637), 
  (68, 41821.073091088445), 
  (69, 40009.78839610813), 
  (70, 38186.113293902934), 
  (71, 36350.608946433815), 
  (72, 34503.84066411754), 
  (73, 32646.37771668947), 
  (74, 30778.793141380764), 
  (75, 28901.663548517958), 
  (76, 27015.568924659736), 
  (77, 25121.09243338852), 
  (78, 23218.82021388006), 
  (79, 21309.34117737897), 
  (80, 19393.246801710808), 
  (81, 17471.130923967106), 
  (82, 15543.589531501515), 
  (83, 13611.220551380227), 
  (84, 11674.623638431713), 
  (85, 9734.39996204431), 
  (86, 7791.151991862593), 
  (87, 5845.48328253557), 
  (88, 3897.99825767181), 
  (89, 1949.301993158165), 
  (90, 0))

test1()
test2()