3
Doporučeno
1 hlasy
SQL&GG04: Vincentyho metoda
geoget návod statistika
Každý člověk je od přírody líný a rád si ulehčuje práci. Stejně je tomu i v IT branži. A tedy i při práci s databázemi. Jak známo, relačních databází se ptáme na odpověď jazykem SQL. A stejně jako čeština má svá pravidla českého pravopisu, tak má i SQL svůj standard. A tento standard má své revize. A stejně jako se u češtiny občas mění pravidla (asi nejčastěji Velká a malá písmena), tak se v jednotlivých revizích mění i SQL.
Nutno podotknout, že na rozdíl od češtiny se jazyk SQL postupně spíše rozšiřuje, ale již existující pravidla se takřka nemění. A již v roce 1999 byl do SQL jazyka přidán operátor WITH. A tím se dostáváme zpět k té lenosti a já hned ukážu proč.
Představte si, že chcete v databázi pracovat s výsledkem nějakého dotazu. Například se chcete podívat na kešky, které jste našli v roce 2017 a to buď T5 a/nebo v zahraničí a/nebo to jsou prosté 1/1ky:
SELECT * FROM geocache WHERE dtfound LIKE '2017%' AND terrain >= 5 UNION -- použito, aby nedocházelo k duplicitám SELECT * FROM geocache WHERE dtfound LIKE '2017%' AND country <> 'Czech Republic' UNION SELECT * FROM geocache WHERE dtfound LIKE '2017%' AND difficulty = 1 AND terrain = 1;Tak a teď budete chtít to samé ale pro rok 2015. Musíte na 3 místech přepsat filtr roku. Pokud budete znát klauzuli WITH, tak celý dotaz napíšete následujícícm způsobem:
WITH mustr AS ( SELECT * FROM geocache WHERE dtfound LIKE '2017%' ) SELECT * FROM mustr WHERE terrain >= 5 UNION SELECT * FROM mustr WHERE country <> 'Czech Republic' UNION SELECT * FROM mustr WHERE difficulty = 1 AND terrain = 1;Někteří možná podotknete že to není úplně pěkně pedagogické, už první dotaz by šel napsat tak aby se datumový filtr přepisoval na jednom místě i bez WITH klauzule. Ale princip je snad zřejmý.
Jednotlivé části (dotazy) v klauzuli WITH jde řetězit a mohou se odkazovat na předcházející. Takže bude fungovat například následující (rok lze měnit hned v prvním WITH dotazu):
WITH rok AS ( SELECT * FROM geocache WHERE dtfound LIKE '2015%' ), duben AS ( SELECT 'duben' text, name, dtfound FROM rok WHERE dtfound LIKE '____04%' ), rijen AS ( SELECT 'rijen' text, name, dtfound FROM rok WHERE dtfound LIKE '____10%' ), paty_den_v_mesici AS ( SELECT 'paty den v mesici' text, name, dtfound FROM rok WHERE dtfound LIKE '%05' ), sesteho_rijna AS ( SELECT 'sesteho rijna' text, name, dtfound FROM rijen WHERE dtfound LIKE '%06' ) SELECT * FROM duben UNION SELECT * FROM rijen UNION SELECT * FROM paty_den_v_mesici UNION SELECT * FROM sesteho_rijna;Zde si dovolím krátkou vsuvku o operátoru LIKE - porovnává text s maskou a pokud maska odpovídá textu, tak vrací logickou jedničku (pravdu). Je citlivý na VELKÁ a malá písmena a v masce se dají použít tzv. zástupné znaky:
- _ (podtržítko) - nahrazuje právě jeden znak
- % (procento) - nahrazuje jakýkoli počet jakýchkoli znaků (i žádný znak - prázdný řetězec)
WITH byl ale původně do SQL přidán kvůli něčemu jinému. To něco se nazývá rekurze. To můžete chápat například tak, že WITH umožňuje volat ten samý dotaz v rámci sama sebe. Abstraktno dotažené takřka k dokonalosti. Kdo ví, jak funguje v matematice faktoriál, případně číselné řady, nebo snad i fraktál, tak má trochu výhodu pří chápání rekurze. My ostatní, co máme problém s abstraktními pojmy a matematiku jsme přestali chápat, když se mezi čísla začala motat abeceda si snad vystačíme s následujícím. Budeme chtít vypsat dotazem čísla od 1 do 10:
WITH cyklus AS ( SELECT 1 counter/* začneme od jedničky */ UNION ALL SELECT counter + 1/* budeme přičítat jedničku */ FROM cyklus/* voláme stejný dotaz */ WHERE counter < 10/* dokud je counter menší než 10, tak voláme další krok */ ) SELECT counter FROM cyklus;Schválně si ho zkuste spustit. Dotaz, který nesahá do žádné tabulky a stejně vrátí 10 řádků. Rekurze, iterace, otevřené dveře pro použití SQL jazyka pro dalo by se říci programování. Pokud se zamyslíme, tak budeme schopni postavit jednoduchý iterační cyklus s koncovou podmínkou.
Teď si dovolím krok stranou. Takový malý teoretický úvod k řešené úloze.
V rámci geocachingu často řešíme dvě úlohy které jsou k sobě inverzní. Projekci bodu (výpočet cílového bodu, pokud známe počáteční bod, azimut do cílového a vzdálenost mezi body) a výpočet vzdálenosti/azimutu mezi dvěma danými body. V dnešní době už to nikdo nepočítá v ruce nebo na kalkulačce, ale máme na to apky v mobilech, vzorečky do excelů a umí to i GeoGet.
Obecně se pro tyto výpočty používají dvě aproximace zeměkoule. Jednou z nich je koule. Tam je situace celkem jednoduchá, koule je koule, vzorečky jsou jasně dané, jen ta přesnost v realitě je horší. Na kilometrových vzdálenostech se výpočet oproti realitě může lišit o desítky metrů. A to co? Nechceme.
Druhá aproximace tvaru zeměkoule je elipsoid. Elipsoid je jednoduše řešeno splácnutá koule. Ne o moc a jen v jednom směru. Rovník je pravidelný kruh, ale vzdálenost mezi póly je menší než vzdálenost mezi dvěma protilehlými body na rovníku. V geocachingu je velmi oblíbený elipsoid z geodetického standardu WGS84. A protože elipsoid není koule, tak se podstatně liší i výpočet výše zmíněných úloh. Už si nevystačíme se vzorečky pro kouli, ale pokud budeme počítat správně, tak se dostaneme k hodnotám, které budou velmi (spíš VELMI) odpovídat realitě. V dnešní době rychlých počítačů nejsou potřebné kalkulace nijak pomalé. Ale když se na elipsoidech začalo počítat, tak ještě všechny výpočty probíhaly buď ručně za pomoci tabulek goniometrických funkcí, nebo na počítačích s velmi omezeným výkonem.
A právě tehdy si pan Vincenty dal práci a právě pro výpočetní stroje připravil algoritmus, který byl iterační, po každé iteraci kontroloval teoretickou možnou odchylku a skončil, když ta odchylka byla přijatelná. Jinými slovy, když funkce konvergovala dostatečně blízko k očekávané hodnotě. Prostě aplikace fyzikálního zaokrouhlování a zanedbávání v aplikované matematice - znám o tom jeden vtip, chcete ho slyšet? Jako přijatelná odchylka se mu jevila hodnota odpovídající zlomkům milimetru na vzdálenosti desitek kilometrů (a odpovídající v azimutu). Elipsoid je totiž také jakýmsi zprůměrováním něčeho čemu se říká geoid a co ještě více odpovídá skutečnému tvaru zeměkoule.
Vincentyho metoda (VM), jak se způsob výpočtu nazývá je použita v GeoGetu, ve Statoru. Sama o sobě není nějak programátorsky náročná na přípravu, pokud ji děláte v programovací jazyce. Pár podmínek, jeden cyklus, spousta goniometrických funkcí. Vše celkem dobře na internetu popsáno pro různé programovací jazyky. Ne tak ale aplikace do SQL. V něm není jednoduchý cyklus a v SQLite v základu ani goniometrické funkce. Navíc tak hromada vzorečků - jen samotný přepis je chůze po laně nad propastí - jedna chybka a musíte hledat kde jste zapomněli znaménko mínus nebo závorku. Uvidíte dále.
Pokud by vás zajímalo teoretické pozadí VM, tak doporučuju wiki, případně původní článek a nebo programátorské zpracování problému. Pro mě osobně je to trochu španělská vesnice.
Proč se samotnou implementací VM do SQL vůbec zabývat, když v Geogetu i Statoru je už naimplementováno. Jak zaznělo na gc.cz fóru: "až se bude hledat případ nezměrného úsilí vynaloženého za zcela nicotným cílem, tohle hraní s SQL musí skončit na předních místech.". Berte to jako chození do fitka. Mozek, stejně jako svaly nicneděláním krní a atrofuje. Takže následující berte prosím jako praktické mentální cvičení "v rámci boje s tím Němcem, co nám furt schovává věci". A navíc jsem chtěl dostat stejné výsledky jako dává Stator.
Ale zpět k SQL. Výše jsme si ukázali jednoduchý dotaz, který nám vrátí čísla od 1 do 10. Ale co když budeme chtít v jednom kroku provést nějaký složitější výpočet. Není problém. Není problém ani nastavit počáteční parametry. A dotaz může vypadat třeba takto:
WITH parametry AS (/* nastavení počátečních parametrů výpočtu */ SELECT 100 pocatecni_hodnota, 1.143 pocatecni_delitel ), cyklus AS ( SELECT pocatecni_hodnota hodnota, pocatecni_delitel delitel FROM parametry UNION ALL SELECT hodnota / delitel value,/* kalkulace v rámci kroku cyklu */ delitel + 0.04 FROM cyklus WHERE value > 1e-10/* ukončující podmínka - funguje jako DO-WHILE cyklus podmínka se vyhodnocuje až na konci kroku a pokud platí, tak proběhne další krok */ ) SELECT c.hodnota, c.delitel, p.pocatecni_delitel, c.hodnota * p.pocatecni_delitel * 1000/* lze kalkulovat s oběma dotazy z WITH klauzule */ FROM cyklus c CROSS JOIN/* připojí dotaz s parametry */ parametry p;Dotaz postupně dělí počáteční hodnotu dělitelem, který se v každém kroku také mění. A dělá to tak dlouho, než výsledek dělení není dostatečně malý. Takže máme cyklus, máme průběžně měněné parametry výpočtu a umíme i cyklus ve vhodnou chvíli zastavit. Naplňme tedy tuto kostru algoritmem Vincentyho metody, prozatím jen úlohou pro výpočet azimutu / vzdálenosti dvou bodů:
WITH parametry AS (
SELECT id,
7 cnt,/* počet Vincentyho iterací - lambdaDiff by mělo být <1E-12 */
49.74375 x1,/* x1, y1 souřadnice bodu odkud počítám vzdálenost */
15.3386333 y1,/* zde použit střed ČR https://coord.info/GC4PCQB */
x x2,
y y2,
6378.137 a,
6356.752314245 b/* a, b - paramertry referenčního elipsoidu WGS84 v km */
FROM geocache
WHERE dtfound <> 0/* vybere nalezené keše */
),
cyklus AS (
SELECT id,/* úvodní výpočty - názvosloví z http://www.movable-type.co.uk/scripts/latlong-vincenty.html */
1 cnt,
a,
b,
( (a - b) / a) f,
x1,
y1,
x2,
y2,
(x1 * PI() / 180) fi1,
(y1 * PI() / 180) lambda1,
(x2 * PI() / 180) fi2,
(y2 * PI() / 180) lambda2,
( (y2 * PI() / 180) - (y1 * PI() / 180) ) L,
(1 - ( (a - b) / a) ) * tan( (x1 * PI() / 180) ) tanU1,
(1 - ( (a - b) / a) ) * tan( (x2 * PI() / 180) ) tanU2,
1 / (sqrt(1 + ( (1 - ( (a - b) / a) ) * tan( (x1 * PI() / 180) ) ) * ( (1 - ( (a - b) / a) ) * tan( (x1 * PI() / 180) ) ) ) ) cosU1,
1 / (sqrt(1 + ( (1 - ( (a - b) / a) ) * tan( (x2 * PI() / 180) ) ) * ( (1 - ( (a - b) / a) ) * tan( (x2 * PI() / 180) ) ) ) ) cosU2,
( (1 - ( (a - b) / a) ) * tan( (x1 * PI() / 180) ) ) * (1 / (sqrt(1 + ( (1 - ( (a - b) / a) ) * tan( (x1 * PI() / 180) ) ) * ( (1 - ( (a - b) / a) ) * tan( (x1 * PI() / 180) ) ) ) ) ) sinU1,
( (1 - ( (a - b) / a) ) * tan( (x2 * PI() / 180) ) ) * (1 / (sqrt(1 + ( (1 - ( (a - b) / a) ) * tan( (x2 * PI() / 180) ) ) * ( (1 - ( (a - b) / a) ) * tan( (x2 * PI() / 180) ) ) ) ) ) sinU2,
( (y2 * PI() / 180) - (y1 * PI() / 180) ) lambda,
0 sinSigma,
0 cosSigma,
0 sigma,
0 sinAlfa,
0 cosSqAlfa,
0 cos2SigmaM,
0 C,
( (y2 * PI() / 180) - (y1 * PI() / 180) ) lambdaIter,
0 lambdaDiff
FROM parametry
UNION ALL
SELECT id,/* iterační krok */
cnt + 1 cnt,
a,
b,
f,
x1,
y1,
x2,
y2,
fi1,
lambda1,
fi2,
lambda2,
L,
tanU1,
tanU2,
cosU1,
cosU2,
sinU1,
sinU2,
lambdaIter lambda,
(sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) sinSigma,
(sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) cosSigma,
(atan2( (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ), (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) ) ) sigma,
(cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ) sinAlfa,
(1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) cosSqAlfa,
(CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END) cos2SigmaM,
f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) C,
(L + (1 - (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) ) * f * (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ) * ( (atan2( (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ), (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) ) ) + (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) * (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) * ( (CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END) + (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) * (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) * ( -1 + 2 * power( (CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END), 2) ) ) ) ) lambdaIter,
(L + (1 - (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) ) * f * (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ) * ( (atan2( (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ), (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) ) ) + (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) * (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) * ( (CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END) + (f / 16 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) * (4 + f * (4 - 3 * (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) ) * (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) * ( -1 + 2 * power( (CASE/* osetreni pokud se pocita czdalenost na rovniku -> deleni nulou */ WHEN (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) = 0 THEN 0 ELSE ( (sinU1 * sinU2 + cosU1 * cosU2 * cos(lambdaIter) ) - 2 * ( (sinU1 * sinU2) / (1 - power( (cosU1 * cosU2 * sin(lambdaIter) / (sqrt(power( (cosU2 * sin(lambdaIter) ), 2) + power( (cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ), 2) ) ) ), 2) ) ) ) END), 2) ) ) ) ) - lambda lambdaDiff
FROM cyklus
WHERE cnt < (/* iteruje dokud není dosaženo potřebného počtu kroků */
SELECT DISTINCT cnt
FROM parametry
)
),
vincenty AS (
SELECT id,/* výpočty následující po dokončení cyklu */
x1,
y1,
x2,
y2,
distance,
bearing,
floor(bearing) bearing_int,
bearing_fin/* hodnota je ve směru DO cílového bodu - pro zpětný nutno připočíst 180st */,
lambdaDiff
FROM (
SELECT *,
(cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) uSq,
(1 + ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 16384) * (4096 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -768 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (320 - 175 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) bigA,
( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) bigB,
( ( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) * sinSigma * (cos2SigmaM + ( ( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) / 4) * (cosSigma * ( -1 + 2 * power(cos2SigmaM, 2) ) - ( ( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) / 6) * cos2SigmaM * ( -3 + 4 * power(sinSigma, 2) ) * ( -3 + 4 * power(cos2SigmaM, 2) ) ) ) ) deltaSigma,
b * (1 + ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 16384) * (4096 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -768 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (320 - 175 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) * (sigma - ( ( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) * sinSigma * (cos2SigmaM + ( ( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) / 4) * (cosSigma * ( -1 + 2 * power(cos2SigmaM, 2) ) - ( ( ( (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) / 1024) * (256 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * ( -128 + (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) * (74 - 47 * (cosSqAlfa * ( (power(a, 2) - power(b, 2) ) / power(b, 2) ) ) ) ) ) ) / 6) * cos2SigmaM * ( -3 + 4 * power(sinSigma, 2) ) * ( -3 + 4 * power(cos2SigmaM, 2) ) ) ) ) ) distance,
CASE WHEN (180 / PI() * atan2(cosU2 * sin(lambdaIter), cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ) ) < 0 THEN 360 + (180 / PI() * atan2(cosU2 * sin(lambdaIter), cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ) ) ELSE (180 / PI() * atan2(cosU2 * sin(lambdaIter), cosU1 * sinU2 - sinU1 * cosU2 * cos(lambdaIter) ) ) END bearing,
CASE WHEN (180 / PI() * atan2(cosU1 * sin(lambdaIter), -sinU1 * cosU2 + cosU1 * sinU2 * cos(lambdaIter) ) ) < 0 THEN 360 + (180 / PI() * atan2(cosU1 * sin(lambdaIter), -sinU1 * cosU2 + cosU1 * sinU2 * cos(lambdaIter) ) ) ELSE (180 / PI() * atan2(cosU1 * sin(lambdaIter), -sinU1 * cosU2 + cosU1 * sinU2 * cos(lambdaIter) ) ) END bearing_fin
FROM cyklus
WHERE cnt = (/* vybere pouze poslední krok iterace */
SELECT DISTINCT cnt
FROM parametry
)
)
)
SELECT v.*,
g.*
FROM vincenty v
LEFT JOIN
geocache g ON v.id = g.id
ORDER BY dtfound,
dtfoundtime;Tento dotaz vám vrátí nalezené keše a jejich azimut a vzdálenost od geometrického středu ČR. Ano, je to dlouhé jak Lovosice. Ale pokud abstrahujete ty goniometrické vzorečky, tak je ten dotaz velmi podobný předchozímu příkladu. Ty vzorečky jsou dlouhé proto, že SQL neumí na úrovni záznamu definovat alias, takže pokud by měly být kratší, tak by každý alias byl definován v samostatném pod-dotazu. Což by ale bylo ještě nepřehlednější.
Hodnoty jsem ověřoval v online kalkulačce na stánkách NOAA, což považuju za spolehlivou autoritu. Navíc jsem je ověřoval i vůči Statoru a Geogetu a lišíme se někde na 7 až osmém destinném místě (odpovídá odchylce desetinkám milimetrů na desítkách kilometrů vzdálenosti). Není problém tento dotaz použít i pro definování uživatelské SQL funkce v SQLite Studiu (dle návodu v jednom z předchozích dílů).
Závěrem bych chtěl poděkovat Pe_Bo, kterému se mě zželelo a vyrobil velmi rychle linkovatelnou dll knihovnu s implementací obou Vincentyho úloh, tedy i projekci bodu. Nutno podotknout, že použití knihovny (jak přilinkovat knihovnu viz předchozí díl) má nezanedbatelný vliv na výkon dotazu, nehledě na to, že se tím dotaz vracející stejné výsledky stává mnohem přehlednější:
SELECT id, 49.74375 x1,/* x1, y1 souřadnice bodu odkud počítám vzdálenost */ 15.3386333 y1,/* zde použit střed ČR https://coord.info/GC4PCQB */ x x2, y y2, distance(49.74375, 15.3386333, x, y) / 1000 distance,/* funkce v dll počítá s metry */ azimut_start(49.74375, 15.3386333, x, y) bearing, floor(azimut_start(49.74375, 15.3386333, x, y) ) bearing_int, azimut_end(49.74375, 15.3386333, x, y) bearing_fin, g.* FROM geocache g WHERE dtfound <> 0 ORDER BY dtfound, dtfoundtime;Snad jsem vás dneska tím goniometrickým bumbrlíčkem moc neodradil. Osobně jsem tento problém vyčerpal tzv. do hloubky. Vyřešil jsem tím úlohu, která mě lákala více jak 4 roky. Uvidím, jestli mě napadne v brzké době nějaká jiná, začnu pátrat tzv. do šířky nad jiným problémem. Pokud by jste nějaký pěkný nad daty z GeoGetu měli, dejte vědět, za inspiraci budu jedině rád.
EDIT 21.4.2021 - SQLite Studio začalo po takřka třech letech vydávat nové verze. Přidali extension manager, takže se dll knihovny nemusí při každém spuštění programu připojovat load_extension funkcí, ale jen se k nim nastaví cesta a program si je načte automaticky při startu sám. Zároveň ale program přešel kompletně na 64bitovou architekturu a 32bitová verze už není k dispozici. To přineslo problém s knihovnami od Pe_Bo, protože ty byly 32bitové a v nové verzi nepoužitelné. Petr byl ale naprosto skvělý a kompresní i Vincentiho knihovnu takřka obratem na pořádání vyrobil i v 64bitové verzi. Otestoval jsem je, fungují správně a ke stažení jsou na https://drive.google...iew?usp=sharing
GeoGet
URWIGO
CWG
Úžasné. Stále jsem u toho prvního rekurzivního dotazu a bádám, jak to funguje. Pořád mi vychází, že když se volá cyklus, tak ten tím úvodním řádkem SELECT 1 counter musí counter přece vždy znovu nastavit na 1 a tedy tu inkrementaci torpédovat.