Jak Garmin počítá přímou vzdálenost mezi body

Diskuze o software, který má mnoho společného s touto značkou (MapSource, BaseCamp, POI Loader, WebUpdater a ostatní).

Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod Path » 10.03.09 14:17 (Út)

Drtivá většina navigačních přístrojů (nejen značky Garmin) včetně mapového software (například nám dobře známý MapSource) používá k určení vzdálenosti mezi dvěma body matematický algoritmus, který je sice velmi rychlý, ale není naprosto přesný. Algoritmus známý jako Great Circle vychází z předpokladu, že Země je koule o daném poloměru. Počítaná vzdálenost je tedy nejkratší spojnice dvou bodů na kulové ploše (ortodroma), která leží na kružnici se středem v jejím těžišti (střed koule).

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
Nemáte oprávnění prohlížet přílohy.
Naposledy editoval Path dne 28.02.11 10:41 (Po). Počet editací: 5x.
Důvod: Vzorec opraven
Obrázek uživatele
Path
Site Admin
 
Příspěvků: 3406
Registrace: 15.05.07 23:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod zhrd » 10.03.09 15:52 (Út)

Duvod, proc Vam vychazi rozdil mezi Great Circle (76 499,73630 m) a Haversine (76 498,44083 m) je ten, ze mate v excelu chybu ve vzorecku - pred COS(RADIANS(lat_1)) ma byt minus (1-sin-cos). Po oprave vychazi Haversine 76499,7363, tedy stejne jako Great Circle. Ono je to vlastne pochopitelne, oba vzorecky jsou matematicky identicke a vzhledem k tomu, ze Excel pocita na 15 mist, by ani k tak velke diferenci nemelo dochazet.
zhrd
 
Příspěvků: 105
Registrace: 08.04.08 7:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod Path » 10.03.09 17:10 (Út)

Máte pravdu a díky za upozornění, tohle mi opravdu uteklo: (1−a)

Opraveno
Obrázek uživatele
Path
Site Admin
 
Příspěvků: 3406
Registrace: 15.05.07 23:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod zhrd » 10.03.09 22:36 (Út)

Kdyz tak o tom premyslim, tak se mi ten vysledek podle Vincenty nezda. Zarazi mne totiz, ze vysel VETSI nez vzdalenost na kouli. Pri vypoctu na kouli totiz uvazujete jako polomer koule vetsi poloosu elipsoidu a tak bych rekl, ze vzdalenost na elipsoidu by mela vyjit MENSI nez na kouli, ktera jej obklopuje. Ale nemam na to, abych zkontroloval vzorecky. Sice jsem kdysi pocital vzdalenost na elipsoidu pro navigacni system Omega, ale nejenze jsem jiz davno vse zapomnel, ale pravdepodobne i zhloupl.
zhrd
 
Příspěvků: 105
Registrace: 08.04.08 7:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod Path » 11.03.09 13:07 (St)

Teoreticky ano, ale výpočet vzdálenosti na kouli je zatížen chybou. Ani já jsem neměl chuť a čas prozkoumávat algoritmus Vincenty a naprosto jsem se spoléhal na věrohodnost zdroje. Je možné, že bych někde udělal chybu v převodu do VBA, ale zkoušel jsem ověřit výsledky s několika dalšími kalkulátory a vždy mi vyšly totožné.
Obrázek uživatele
Path
Site Admin
 
Příspěvků: 3406
Registrace: 15.05.07 23:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod zhrd » 11.03.09 14:29 (St)

Nezpochybnuji ani spravnost algoritmu, ani Vase programatorske umeni :-). Jen si nejsem jist korektnosti vstupnich hodnot. Jedna a tataz zemepisna sirka na kouli a na elipsoidu jsou totiz v principu ruzna mista. Sirka na kouli je uhel mezi rovnikem a spojnici bodu se stredem (geocentricka sirka). Sirka na elipsoidu je uhel mezi rovnikem a normalou k povrchu elipsoidu (geodeticka sirka). A proto si myslim (ale klidne pripustim, ze se mylim), ze ve skutecnosti porovnavate vzdalenosti mezi ruznymi body. Jinak receno, jestlize spocitam vzdalenost na kouli mezi dvema body s danou geocentrickou sirkou, pak musim prepocitat tyto sirky na geodeticke sirky a ty pak pouzit do algoritmu na vypocet vzdalenosti na elipsoidu.
zhrd
 
Příspěvků: 105
Registrace: 08.04.08 7:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod Path » 11.03.09 15:25 (St)

Vaše připomínka se mi zdá poměrně logická. Osobně se necítím být v dané problematice odborníkem, a už vůbec ne geodet nebo zapálený matematik. Proto jsem také uváděl zdroj, kde je algoritmus uveden a každý si to může překontrolovat. Je dobře, že někdo reaguje. Přiznám se ale, že nerozumím tomu, proč by vstupní hodnoty (šířky a délky) pro výpočet vzdálenosti dle Vincentyho měly být korigovány, když vycházejí ze stejného referenčního elipsoidu, který je v samotném algoritmu definován. To už bych spíš pochopil nutnou transformaci pro výpočet ortodromy? Nahlodal jste mě a celkem mě to zajímá. Pokud budete mít čas a chuť, můžete to více rozebrat?
Obrázek uživatele
Path
Site Admin
 
Příspěvků: 3406
Registrace: 15.05.07 23:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod zhrd » 11.03.09 16:21 (St)

Ono je to jednoduche. Vezmu si jako priklad sirku Prahy tak, jak ji mate ve svem prikladu, tj. 50,08667. Prvni vec, kterou si musim ujasnit je to, v jakem souradnicovem systemu to je. Ve WGS-84? V S-42? JTSK? Nebo je to geocentricka sirka?
Predpokladejme, ze ona sirka je geocentricka (kdyz uz byla pouzita pro vypocet vzdalenosti na kouli). A ted se ptejme, jake je geodeticka sirka tohoto bodu (vysvetleni tech pojmu je napr. http://www.natur.cuni.cz/~bayertom/Mmk/mk0.pdf). Prepocet je tg(geodeticka)=(a/b)^2*tg(geocentricka). Takze na elipsoidu ma tentyz bod (geodetickou) sirku 50,27595566. Stejne bych mel nalozit se sirkou Turnova (50,58655 -> 50,77520333). Pokud tyto sirky pouziji pro vypocet vzdalenosti na elipsoidu, dostanu 76347,646, tedy mensi nez vzdalenost na velke kouli a proto se mi ten vysledek libi vic (nerikam ze je spravny).
zhrd
 
Příspěvků: 105
Registrace: 08.04.08 7:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod Path » 11.03.09 16:52 (St)

Předpokladem pro všechny uvedené výpočty je, že polohy bodu vycházejí z geodetického geocentrického systému WGS-84, který je standardem v systému GPS a implicitně defaultní systém pro navigační přístroje i SW Garmin. Všechny polohy jsou tedy uváděny ve WGS-84 (pokud si uživatel nenastaví něco jiného). Ve algoritmu je rovněž uvedeno o jaký systém se jedná a o jiném systému než WGS-84 jsem v této souvislosti ani neuvažoval.
Obrázek uživatele
Path
Site Admin
 
Příspěvků: 3406
Registrace: 15.05.07 23:15 (Út)

Re: Jak Garmin počítá přímou vzdálenost mezi body

Odeslatod zhrd » 11.03.09 16:59 (St)

Asi si nerozumime. Nemohu pro vsechny vypocty (tj. jak na kouli, tak na elipsoidu) pouzit tytez souradnice. Na kouli musim pouzit geocentricke, na elipsoidu geodeticke.
Pojem "geodetického geocentrického " je protimluv.
Mezitim jsem si udelal kontrolni vypocet. Pokud vezmu vypocet na kouli, ale misto vetsi poloosy pouziji vzdalenost od rovniku, jaky je na elipsoidu v nasi sirce (tj. 6365416 m), tak mi vzdalenost na (ted jiz tecne) kouli vyjde 76347,15991, tedy prakticky stejna jako elipsoidu. Tak si myslim, ze mam pravdu.
zhrd
 
Příspěvků: 105
Registrace: 08.04.08 7:15 (Út)

Další

Zpět na Garmin: Software

Kdo je online

Online uživatelé v tomto fóru: CommonCrawl [Bot], Yandex [Bot] a 2 návštevníci.