Přejít na obsah






Doporučeno
Fotka
* * * * * 1 hlasy

SQL&GG04: Vincentyho metoda

Příspěvek od tarmara , 10 leden 2018 · 5 854 Zobrazení

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)
Pro sofistikovanější prohledávání textu existují ještě regulární výrazy, ale to už by bylo mimo rámec toho seriálu a dovolím si tvrdit (dle hesla podle sebe soudím tebe) i mimo chápání většiny auditoria.

 

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.



  • 6



Fotka
Pontiac_CZ
dub 04 2018 6:43

Úž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. :)

    • 0

Úž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. :)

No já jsem to měl původně podobně. Má to určitě i své technické vysvětlení dané implementací. Ale protože tohle funguje nejen v SQLite, ale třeba i v Oraclu, tak jsem tomu prostě začal věřit. Občas se stane, že pak některé věci nefungují (např. hvězdičková konvence - výpis a možnost referencování všech sloupců přes znak "*"), ale takový už je život...

    • 0

Červenec 2018

P Ú S Č P S N
      1
2345678
9101112131415
161718 19 202122
23242526272829
3031     

Poslední komentáře

Reklama