od Výpočet vzdálenosti (spojnice) dvou bodů je v tomto případě opravdu velmi rychlý, ale při malých vzdálenostech je zatížený chybou danou zaokrouhlením. Chyba výpočtu je v běžném použití navigačních přístrojů naprosto zanedbatelná a menší zátěž procesoru GPS přijímače navíc pozitivně ovlivňuje spotřebu energie. Algoritmus Great Circle můžete najít v souboru Coordinates.xls, který jsem v minulosti poskytl ke stažení zde.
Great Circle (vzdálenost mezi A > B):
R = poloměr Země
d = R * arccos(sin(lat2)*sin(lat1)+cos(lat2)*cos(lat1)*cos(lon2 – lon1))
Garmin používá algoritmus Great Circle nejen v přístrojích, ale také v MapSource. Pokud chcete využít další algoritmy pro výpočet vzdálenosti, např. při analýze tracklogu, můžete také využít o něco přesnější algoritmus Haversine. Tento algoritmus je zatížen chybou ± 3 m na 1 km.
Haversine (vzdálenost mezi A > B):
R = poloměr Země
Δlat = lat2− lat1
Δlong = long2− long1
a = sin²(Δlat/2) + cos(lat1).cos(lat2).sin²(Δlong/2)
c = 2.atan2(√a, √(1−a))
d = R.c
Předchozí dva algoritmy vycházejí z předpokladu, že Země je koule. Zde bylo v minulosti zmíněno, že fyzikální tvar Země není koule, ale geoid, který lze matematicky jen velmi obtížně definovat. Pro náhradu geoidu se jako nejvhodnější těleso hodí rotační elipsoid. A právě rotačního elipsoidu pro vlastní definici výpočtu vzdálenosti využil už v roce 1975 americký geodet Thaddeus Vincenty, jehož algoritmus dokáže spočítat vzdálenost mezi dvěma body s přesností až na 0,5 mm v rámci celého elipsoidu. Přestože tento algoritmus není na první pohled úplně srozumitelný, je dostatečně popsán zde včetně jeho přepracování do javasriptu.
Pro zájemce jsem jsem algoritmus Vincenty převedl z javascriptu do VBA:
Vincenty (vzdálenost mezi A > B):
- vba Code: Vybrat vše
'funkce pro výpočet dvou zeměpisných bodů pomocí algoritmu Thaddeus Vincenty
'přesnost výpočtu na teoretickém elipsoidu činí 0,5 mm nebo 0.000015"
' http://en.wikipedia.org/wiki/Thaddeus_Vincenty
' http://www.movable-type.co.uk/scripts/l ... centy.html
Public Function fceVzalenostVinc(lat1, lon1, lat2, lon2)
Dim a, b, f, L, Pi
Dim U1, U2, sinU1, sinU2, cosU1, cosU2
Dim lambda, lambdaP
Dim iIterace
If Len(lat2) = 0 Or Len(lon2) = 0 Then
fceVzalenostVinc = 0
Exit Function
End If
'=== elipsoid WGS-84 =========================================
a = 6378137 ' velká poloosa
b = 6356752.314245 ' malá poloosa
f = 1 / 298.257223563 ' reciproká hodnota zloštění
L = ToRad(lon2 - lon1) ' rozdíl zeměpisných délek v radiánech
'==============================================================
Pi = Application.WorksheetFunction.Pi
U1 = Atn((1 - f) * Tan(ToRad(lat1)))
U2 = Atn((1 - f) * Tan(ToRad(lat2)))
sinU1 = Sin(U1)
cosU1 = Cos(U1)
sinU2 = Sin(U2)
cosU2 = Cos(U2)
lambda = L
lambdaP = 2 * Pi
iIterace = 20
Do
iIterace = iIterace - 1
sinLambda = Sin(lambda)
cosLambda = Cos(lambda)
sinSigma = Sqr((cosU2 * sinLambda) * (cosU2 * sinLambda) + _
(cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * (cosU1 * sinU2 - sinU1 * cosU2 * cosLambda))
If sinSigma = 0 Then
fceVzalenostVinc= 0 'stávající bod
Exit Function
End If
cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda
sigma = WorksheetFunction.Atan2(cosSigma, sinSigma)
sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma
cosSqAlpha = 1 - sinAlpha * sinAlpha
If cosSqAlpha = 0 Then
cos2SigmaM = 0 'pro rovník je cosSqAlpha = 0
Else
cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha
End If
C = f / 16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha))
lambdaP = lambda
lambda = L + (1 - C) * f * sinAlpha * _
(sigma + C * sinSigma * (cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)))
Loop While (Abs(lambda - lambdaP) > 10 ^ -12) And iIterace > 0
If iIterace = 0 Then
fceVzalenostVinc = CVErr(xlErrNA) 'selhání konvergence
Exit Function
End If
uSq = cosSqAlpha * (a * a - b * b) / (b * b)
AA = 1 + uSq / 16384 * (4096 + uSq * (-768 + uSq * (320 - 175 * uSq)))
BB = uSq / 1024 * (256 + uSq * (-128 + uSq * (74 - 47 * uSq)))
deltaSigma = BB * sinSigma * (cos2SigmaM + BB / 4 * (cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM) - _
BB / 6 * cos2SigmaM * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM * cos2SigmaM)))
s = Round(b * AA * (sigma - deltaSigma), 3) 'zaokrouhlíme na tři desetinná místa
fceVzalenostVinc = s
End Function
'pomocná funkce
Function ToRad(degrees)
ToRad = WorksheetFunction.Radians(degrees)
End Function
Test
Všechny tři algoritmy jsem použil na výpočet přímé vzdálenosti mezi Prahou a Turnovem s těmito výsledky:
Great Circle: 76 499,73630 m
Haversine: 76 499,73630 m
Vincenty: 76 540,94800 m
Vincenty spočítal téměř o 43 metrů jinou vzdálenost než předchozí dva algoritmy. Uvažujete-li o vlastní podrobnější analýze tracklogů (import do excelu zde), můžete pro výpočet vzdáleností využít některého z výše uvedených algoritmů. Nejvhodnější algoritmus bude určitě Vincenty.
Výše uvedené algoritmy jsou k dispozici zde (MS Excel):
Vzdalenosti_MapSource.rar- (14.67 KB) staženo 613 krát
